5. pySALEplotのサンプル
pySALEPlotの解説&実習 20161109pySALEPlot.pdf の実践編に掲載されてあるサンプル集
- DenTmp.py
- profile_example1.py
- トレーサー粒子に関するプロットサンプル
DenTmp.py
密度と温度のプロット
import pySALEPlot as psp import matplotlib.pyplot as plt from numpy import arange,sqrt,ma # Need this for the colorbars we will make on the mirrored plot from mpl_toolkits.axes_grid1 import make_axes_locatable # This example plotting script designed to plot # material and temperature in the Chicxulub example # If viridis colormap is available, use it here try: plt.set_cmap('viridis') except: plt.set_cmap('YlGnBu_r') # distances between tracers def get_distances(s,line): x=s.xmark[line] y=s.ymark[line] return sqrt((x[:-1]-x[1:])**2+(y[:-1]-y[1:])**2) # Define the maximum separation allowed when plotting lines maxsep=3. def make_colorbars(ax, p, f, units): # Create axes either side of the plot to place the colorbars divider = make_axes_locatable(ax) cx1 = divider.append_axes("right", size="5%", pad=0.4) cx2 = divider.append_axes("left", size="5%", pad=0.8) cb1 = fig.colorbar(p[0], cax=cx1) cb1.set_label(psp.longFieldName(f[0])+units[0]) cb2 = fig.colorbar(p[1], cax=cx2) cb2.set_label(psp.longFieldName(f[1])+units[1]) # Need to set the labels on the left for this colorbar cx2.yaxis.tick_left() cx2.yaxis.set_label_position('left') # Make an output directory dirname='DenTmp' psp.mkdir_p(dirname) # Open the datafile model=psp.opendatfile('jdata.dat') # Set the distance units to km model.setScale('km') # Set up a pylab figure fig=plt.figure(figsize=(8,4)) ax=fig.add_subplot(111,aspect='equal') # Loop over timesteps for i in arange(0,model.nsteps,2): # Set the axis labels ax.set_xlabel('r [km]') ax.set_ylabel('z [km]') # Set the axis limits ax.set_xlim([-100,100]) ax.set_ylim([-50,45]) # Read the time step 'i' from the datafile: # read two or more fields by making a list of their abbreviations step=model.readStep(['Tmp', 'Den'],i) # -- the classic iSALEPlot mirrored setup. Plot the second field # -- using negative x values p1=ax.pcolormesh(model.x,model.y,step.Den, cmap='Oranges') #,vmin=1,vmax=model.nmat+1) p2=ax.pcolormesh(-model.x,model.y,step.data[0], vmin=0,vmax=1500) # Material boundaries [ax.contour(model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]] [ax.contour(-model.xc,model.yc,step.cmc[mat],1,colors='k',linewidths=0.5) for mat in [0,1,2]] # Tracer lines for u in range(1,model.tracer_numu): tru=model.tru[u] # Plot the tracers in horizontal lines, every 5 lines for l in arange(0,len(tru.xlines),5): # Get the distances between pairs of tracers in xlines dist=get_distances(step,tru.xlines[l]) # Mask the xmark values if separation too big... means the line won't be connected here ax.plot(ma.masked_array(step.xmark[tru.xlines[l]][:-1],mask=dist > maxsep*tru.d[0]), step.ymark[tru.xlines[l]][:-1], c='#808080',marker='None',linestyle='-',linewidth=0.5) # Add colorbars; only need to do this once if i == 0: make_colorbars(ax, [p1, p2], [step.plottype[1], step.plottype[0]], [r'[kg/m$^3$]', '[K]']) ax.set_title('{: 5.2f} s'.format(step.time)) # Save the figure fig.savefig('{}/DenTmp{:05d}.png'.format(dirname,i),dpi=300) # Remove the field, ready for the next timestep to be plotted ax.cla()
profile_example1.py
プロットするX座標の変更
import pySALEPlot as psp import matplotlib.pyplot as plt # This example plotting script designed to plot # a vertical pressure profile of the Chicxulub example # Open the datafile model=psp.opendatfile('jdata.dat') # Set the distance units to mm model.setScale('km') # Set up a pyplot figure fig=plt.figure() ax=fig.add_subplot(111) # Set the axis labels ax.set_xlabel('Pressure [GPa]') ax.set_ylabel('Depth [km]') # Set the axis limits ax.set_ylim(-200,0) #ax.set_ylim(0,100) # Read the time steps from the datafile step0=model.readStep('Pre',0) step1=model.readStep('Pre',1) # Plot the pressure field ax.plot(step1.Pre[0,:]*1e-9,model.yc[0,:],'k:',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[0,:][0])) ax.plot(step1.Pre[50,:]*1e-9,model.yc[0,:],'m--',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[50,:][0])) ax.plot(step1.Pre[100,:]*1e-9,model.yc[0,:],'b-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[100,:][0])) ax.plot(step1.Pre[200,:]*1e-9,model.yc[0,:],'r-',label='{:3.1f} s x = {:3.1f} km'.format(step1.time,model.xc[200,:][0])) # Add a legend ax.legend(loc=0) # Save the figure fig.savefig('profile1.png')
profile_example2.py
トレーサー粒子の最大圧力(TrP)をプロット
import pySALEPlot as psp import matplotlib.pyplot as plt # This example plotting script designed to plot # a vertical pressure profile of the Chicxulub example # Open the datafile model=psp.opendatfile('jdata.dat') # Set the distance units to mm model.setScale('km') # Set up a pyplot figure fig=plt.figure() ax=fig.add_subplot(111) # Set the axis labels ax.set_xlabel('r [km]') ax.set_ylabel('z [km]') # Set the axis limits ax.set_xlim(0,100) ax.set_ylim(-50,45) # Read the time steps from the datafile nstep=0 step0=model.readStep('TrP',0) step1=model.readStep('TrP',200) # Plot the pressure field for u in range(model.tracer_numu): #loop over tracer clouds tstart = model.tru[u].start tend = model.tru[u].end scat = ax.scatter(step.xmark[tstart:tend],step.ymark[tstart:tend], c=step1.TrP[angle_range]*1e-9,vmin=0,vmax=30, cmap='viridis',s=4,linewidths=0) cb=fig.colorbar(scat) cb.set_label('Tracer Prak Pressure [GPa]') # Save the figure fig.savefig('profile2.png')
profile_example3.py
ある角度でのトレーサー粒子の最大圧力(TrP)をプロット
import pySALEPlot as psp import matplotlib.pyplot as plt import numpy as np # This example plotting script designed to plot # a vertical pressure profile of the Chicxulub example # Open the datafile model=psp.opendatfile('jdata.dat') # Set the distance units to mm model.setScale('km') # Set up a pyplot figure fig=plt.figure() ax=fig.add_subplot(111) # Set the axis labels ax.set_xlabel('r [km]') ax.set_ylabel('z [km]') # Set the axis limits ax.set_xlim(0,100) ax.set_ylim(-50,45) # Read the time steps from the datafile nstep=0 step0=model.readStep('TrP',0) step1=model.readStep('TrP',200) # Calculate angles of tracers angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark)) angle = -45. angle_range = np.where((angles > angle) & (angles < angle + 1.)) # Plot the pressure field scat = ax.scatter(step0.xmark[angle_range],step0.ymark[angle_range], c=step1.TrP[angle_range]*1e-9,vmin=0,vmax=30, cmap='viridis',s=4,linewidths=0) cb=fig.colorbar(scat) cb.set_label('Tracer Prak Pressure [GPa]') # Save the figure fig.savefig('profile3.png')
profile_example4.py
ある角度でのトレーサー粒子の最大圧力(TrP)と距離をプロット
import pySALEPlot as psp import matplotlib.pyplot as plt # This example plotting script designed to plot # a vertical pressure profile of the Chicxulub example # Open the datafile model=psp.opendatfile('jdata.dat') # Set the distance units to mm model.setScale('km') # Set up a pyplot figure fig=plt.figure() ax=fig.add_subplot(111) # Set the axis labels ax.set_xlabel('Distance from impact point [km]') ax.set_ylabel('Peak Pressure [GPa]') # Set the axis limits ax.set_xlim(0,30) ax.set_ylim(0,60) # Read the time steps from the datafile nstep=0 step0=model.readStep('TrP',0) step1=model.readStep('TrP',200) # Calculate angles of tracers angles = np.rad2deg(np.arctan(step0.xmark/step0.ymark)) angle = -45. angle_range = np.where((angles > angle) & (angles < angle + 1.)) r = np.sqrt(step0.xmark[angle_range]**2+step0.ymark[angle_range]**2) trp = step1.TrP[angle_range]*1e-9 ax.plot(np.sort(r),trp[np.argsort(r)],'k-') # Save the figure fig.savefig('profile4.png')
Keyword(s):
References:[SideMenu]