# energy.py # # Plots the probability distribution for the dynamo to show # energy releases of certain strengths. # import pylab as plt import numpy as np import pencil as pc from matplotlib.image import NonUniformImage #################################################################################################### ### making own histogram #################################################################################################### #dirs = ['../A64kf2e_rep', '../A64kf2a_rep', '../A64kf2b_rep', '../A64kf2c_rep', '../A64kf2d_rep', #'../A64_kf1a_rep', '../A64_kf1b_rep', '../A64_kf1c_rep', '../A64_kf1d_rep', '../A64_kf1e_rep'] #CAnew = 2.2*1.54135/2.23080 #Calp = np.array([2.2, 2.2, 2.2, 2.2, 2.2, CAnew, CAnew, CAnew, CAnew, CAnew]) #CS = np.array([.24, .50, 1.1, 5.6, 12.4, .24, .50, 1.1, 5.6, 12.4]) #D = Calp*CS #jpt = np.array([]) #i = 0 #tTot = 0. #for d in dirs: #ts = pc.read_ts(datadir = d+"/data", quiet = True, plot_data = False) #idx0 = len(jpt) #jpt = np.resize(jpt, (idx0+3*len(ts.t),2)) #jpt[idx0+0*len(ts.t):idx0+1*len(ts.t),0] = ts.jxpt #jpt[idx0+1*len(ts.t):idx0+2*len(ts.t),0] = ts.jypt #jpt[idx0+2*len(ts.t):idx0+3*len(ts.t),0] = ts.jzpt #jpt[idx0+0*len(ts.t):idx0+3*len(ts.t),1] = D[i] #jpt[idx0+1*len(ts.t):idx0+3*len(ts.t),1] = D[i] #jpt[idx0+2*len(ts.t):idx0+3*len(ts.t),1] = D[i] #tTot += ts.t[-1]-ts.t[0] #i = i + 1 ## normalize the dissipation (epsilon/\bar{espilon}) #jpt[:,0] = jpt[:,0]/abs(jpt[:,0]).mean() #y_bins = 30 ##logE = np.linspace(np.log10(abs(jpt[:,0]).min()), np.log10(abs(jpt[:,0]).max()), y_bins) #logE = np.linspace(np.log(abs(jpt[:,0]).min()), np.log(abs(jpt[:,0]).max()), y_bins) #B = D #D = np.sort(D) ## make 2d histogram #x_edges = np.zeros(len(D)+1) ##x_edges[1:-1] = (np.log10(D[1:]) - np.log10(D[0:-1]))/2 + np.log10(D[0:-1]) ##x_edges[0] = -(np.log10(D[1]) - np.log10(D[0]))/2 + np.log10(D[0]) ##x_edges[-1] = (np.log10(D[-1]) - np.log10(D[-2]))/2 + np.log10(D[-1]) #x_edges[1:-1] = (np.log(D[1:]) - np.log(D[0:-1]))/2 + np.log(D[0:-1]) #x_edges[0] = -(np.log(D[1]) - np.log(D[0]))/2 + np.log(D[0]) #x_edges[-1] = (np.log(D[-1]) - np.log(D[-2]))/2 + np.log(D[-1]) #y_edges = np.zeros(len(logE)+1) #y_edges[1:-1] = (logE[1:] - logE[0:-1])/2 + logE[0:-1] #y_edges[0] = -(logE[1] - logE[0])/2 + logE[0] #y_edges[-1] = (logE[-1] - logE[-2])/2 + logE[-1] #H, xe, ye = np.histogram2d(jpt[:,1], abs(jpt[:,0]), normed = True, bins = [10**x_edges, 10**y_edges]) #H, xe, ye = np.histogram2d(np.log10(jpt[:,1]), np.log10(abs(jpt[:,0])), normed = True, bins = [x_edges, y_edges]) H, xe, ye = np.histogram2d(np.log(jpt[:,1]), np.log(abs(jpt[:,0])), normed = True, bins = [x_edges, y_edges]) # prepare the plot width = 8 height = 6 plt.rc("figure.subplot", left=(100/72.27)/width) #plt.rc("figure.subplot", right=(width-30/72.27)/width) plt.rc("figure.subplot", right=(width-50/72.27)/width) plt.rc("figure.subplot", bottom=(80/72.27)/height) plt.rc("figure.subplot", top=(height-20/72.27)/height) figure = plt.figure(figsize = (width, height)) axes = plt.subplot(111) plt.ion() # plot #plt.imshow(np.log10(zip(*H)), extent = [10**((xe[0]+xe[1])/2), 10**((xe[-2]+xe[-1])/2), 10**((ye[0]+ye[1])/2), 10**((ye[-2]+ye[-1])/2)], interpolation = 'nearest', origin = 'lower', aspect = 'auto', cmap = 'gnuplot2') #image = NonUniformImage(axes, interpolation = 'nearest', cmap = 'gnuplot2', extent = [10**((xe[0]+xe[1])/2), 10**((xe[-2]+xe[-1])/2), 10**((ye[0]+ye[1])/2), 10**((ye[-2]+ye[-1])/2)]) image = NonUniformImage(axes, interpolation = 'nearest', cmap = 'gnuplot2', extent = [(xe[0]+xe[1])/2, (xe[-2]+xe[-1])/2, (ye[0]+ye[1])/2, (ye[-2]+ye[-1])/2]) #image.set_data(np.log10(D), logE, np.log10(zip(*H))) image.set_data(np.log(D), logE, np.log10(zip(*H))) plt.plot([-2,4], [0,0], linewidth = 2, linestyle = '--', color = 'black') axes.images.append(image) #plt.xlim(xmin = np.min(np.log10(D)), xmax = np.max(np.log10(D))) #plt.ylim(ymin = np.min(logE), ymax = np.max(logE)) plt.xlim(xmin = x_edges[0], xmax = x_edges[-1]) plt.ylim(ymin = y_edges[0], ymax = y_edges[-1]) #image.set_data(D, 10**logE, np.log10(zip(*H))) #axes.images.append(image) ##axes.set_xscale('log', basex = 10) ##axes.set_yscale('log', basey = 10) #plt.xlim(xmin = 10**x_edges[0], xmax = 10**x_edges[-1]) #plt.ylim(ymin = 10**y_edges[0], ymax = 10**y_edges[-1]) #image.set_data(np.log10(D), logE, np.log10(zip(*H))) #axes.images.append(image) ##axes.set_xscale('log', basex = 10) ##axes.set_yscale('log', basey = 10) #plt.xlim(xmin = x_edges[0], xmax = x_edges[-1]) #plt.ylim(ymin = y_edges[0], ymax = y_edges[-1]) ##plt.xlim(1e-2,1e2) ##plt.ylim(ymin = 1e-11, ymax = 1e1) plt.draw() # make plot pretty #plt.xlabel(r'$D$', size = 25) #plt.ylabel(r'$\epsilon/\bar{\epsilon}$', size = 25) #plt.xlabel(r'$\log_{10}(D)$', size = 25) #plt.ylabel(r'$\log_{10}(\epsilon/\bar{\epsilon})$', size = 25) plt.xlabel(r'$\ln(D)$', size = 25) plt.ylabel(r'$\ln(\epsilon/\bar{\epsilon})$', size = 25) # colorbar cb = plt.colorbar(image) #cb.set_label(r'$\log_{10}\left(p(D,E)\right)$', fontsize = 25) #cb.set_label(r'$\log_{10}\left(p\left[\log_{10}(D),\log_{10}(\epsilon/\bar{\epsilon})\right]\right)$', fontsize = 25) cb.set_label(r'$\log_{10}(p)$', fontsize = 25) cbytick_obj = plt.getp(cb.ax.axes, 'yticklabels') plt.setp(cbytick_obj, fontsize = 15, family = 'serif') plt.xticks(fontsize=20, family = 'serif') plt.yticks(fontsize=20, family = 'serif') axes.tick_params(axis = 'both', which = 'major', length = 8) axes.tick_params(axis = 'both', which = 'minor', length = 4) # create an offset between the xylabels and the axes for label in axes.xaxis.get_ticklabels(): label.set_position((0,-0.03)) for label in axes.yaxis.get_ticklabels(): label.set_position((-0.03,0)) plt.show() ##################################################################################################### ### 1d pdf with the above data ##################################################################################################### # prepare the plot width = 8 height = 6 plt.rc("figure.subplot", left=(100/72.27)/width) #plt.rc("figure.subplot", right=(width-30/72.27)/width) plt.rc("figure.subplot", right=(width-50/72.27)/width) plt.rc("figure.subplot", bottom=(80/72.27)/height) plt.rc("figure.subplot", top=(height-20/72.27)/height) figure = plt.figure(figsize = (width, height)) axes = plt.subplot(111) plt.ion() #HLin, xe = np.histogram(abs(jpt[:,0]), density = True, bins = 10**y_edges) #plt.loglog(10**logE, HLin, linewidth = 2) #plt.xlim(xmin = 1e-7, xmax = 10) #plt.ylim(ymin = 1e-7, ymax = 100) #HLin, xe = np.histogram(np.log10(abs(jpt[:,0])), density = True, bins = y_edges) HLin, xe = np.histogram(np.log(abs(jpt[:,0])), density = True, bins = y_edges) plt.semilogy(logE, HLin, linewidth = 2) plt.plot([0,0], [1e-6,1], linewidth = 2, linestyle = '--', color = 'black') plt.xlim(xmin = -14, xmax = 4) plt.ylim(ymin = 1e-6, ymax = 1) # make plot pretty #plt.xlabel(r'$\epsilon/\bar{\epsilon}$', size = 25) #plt.ylabel(r'$p(\epsilon/\bar{\epsilon})$', size = 25) #plt.xlabel(r'$\log_{10}(\epsilon/\bar{\epsilon})$', size = 25) #plt.ylabel(r'$p\left[\log_{10}(\epsilon/\bar{\epsilon})\right]$', size = 25) plt.xlabel(r'$\ln(\epsilon/\bar{\epsilon})$', size = 25) plt.ylabel(r'$p\left[\ln(\epsilon/\bar{\epsilon})\right]$', size = 25) plt.xticks(fontsize=20, family = 'serif') plt.yticks(fontsize=20, family = 'serif') axes.tick_params(axis = 'both', which = 'major', length = 8) axes.tick_params(axis = 'both', which = 'minor', length = 4) # create an offset between the xylabels and the axes for label in axes.xaxis.get_ticklabels(): label.set_position((0,-0.03)) for label in axes.yaxis.get_ticklabels(): label.set_position((-0.03,0)) plt.show() ##################################################################################################### ### nuE-1 vs. E ##################################################################################################### # prepare the plot width = 8 height = 6 plt.rc("figure.subplot", left=(100/72.27)/width) #plt.rc("figure.subplot", right=(width-30/72.27)/width) plt.rc("figure.subplot", right=(width-50/72.27)/width) plt.rc("figure.subplot", bottom=(80/72.27)/height) plt.rc("figure.subplot", top=(height-20/72.27)/height) figure = plt.figure(figsize = (width, height)) axes = plt.subplot(111) plt.ion() #HLin2, xe = np.histogram(abs(jpt[:,0]), bins = 10**y_edges) HLin2, xe = np.histogram(abs(jpt[:,0]), bins = np.exp(y_edges)) #plt.plot(10**logE, HLin2/10**logE, linestyle = '-', linewidth = 2, markersize = 0, markeredgewidth = 1.5, color = 'black', drawstyle = 'steps-mid') plt.plot(np.exp(logE), HLin2/np.exp(logE), linestyle = '-', linewidth = 2, markersize = 0, markeredgewidth = 1.5, color = 'black', drawstyle = 'steps-mid') plt.gca().set_xscale('log') plt.gca().set_yscale('log') # make plot pretty plt.xlabel(r'$\epsilon/\bar{\epsilon}$', size = 25) plt.ylabel(r'$\nu \times (\epsilon/\bar{\epsilon})^{-1}$', fontsize = 25) #x = np.array([4e33, 1e36]) #y = 3e15**(1.8)*x**(-1.8)/10. #plt.plot(x, y, linestyle = '--', linewidth = 2, color = 'k') #plt.annotate(r'$E^{-1.8}$', xy=(4e34, 1e-35), xycoords='data', fontsize = 25) #plt.xlim(xmin = 6e31, xmax = 1e36) #plt.ylim(ymin = 3e-38, ymax = 5e-35) plt.xticks(fontsize=20, family = 'serif') plt.yticks(fontsize=20, family = 'serif') axes.tick_params(axis = 'both', which = 'major', length = 8) axes.tick_params(axis = 'both', which = 'minor', length = 4) # create an offset between the xylabels and the axes for label in axes.xaxis.get_ticklabels(): label.set_position((0,-0.03)) for label in axes.yaxis.get_ticklabels(): label.set_position((-0.03,0)) plt.show() ##################################################################################################### ### using Axel's data ##################################################################################################### ## read the data from Axel's file #data = np.loadtxt('../idl/ppdf_form.dat') ## prepare the plot #width = 8 #height = 6 #plt.rc("figure.subplot", left=(100/72.27)/width) ##plt.rc("figure.subplot", right=(width-30/72.27)/width) #plt.rc("figure.subplot", right=(width-50/72.27)/width) #plt.rc("figure.subplot", bottom=(80/72.27)/height) #plt.rc("figure.subplot", top=(height-20/72.27)/height) #figure = plt.figure(figsize = (width, height)) #axes = plt.subplot(111) #plt.ion() ## plot #plt.imshow(data, aspect = 'auto', interpolation = 'nearest', origin = 'lower', #extent = [0.528000, 27.28, 10**-4.95500, 10**3.95500], cmap = 'gnuplot2') #plt.gca().set_xscale('log') #plt.gca().set_yscale('log') ## make plot pretty #plt.xlabel(r'$D$', size = 25) #plt.ylabel(r'$E$', size = 25) ## colorbar #cb = plt.colorbar() #cb.set_label(r'$\log_{10}(p)$', fontsize = 25) #cbytick_obj = plt.getp(cb.ax.axes, 'yticklabels') #plt.setp(cbytick_obj, fontsize = 15, family = 'serif') #plt.xticks(fontsize=20, family = 'serif') #plt.yticks(fontsize=20, family = 'serif') #axes.tick_params(axis = 'both', which = 'major', length = 8) #axes.tick_params(axis = 'both', which = 'minor', length = 4) ## create an offset between the xylabels and the axes #for label in axes.xaxis.get_ticklabels(): #label.set_position((0,-0.03)) #for label in axes.yaxis.get_ticklabels(): #label.set_position((-0.03,0)) #plt.show()