diff --git a/gscripts/ghelp b/gscripts/ghelp index 0039a63..f7f5557 100755 --- a/gscripts/ghelp +++ b/gscripts/ghelp @@ -1,214 +1,216 @@ #!/usr/bin/env python print """ ################################################ gscripts/ ################################################ ################################################ # run preparation galloc : determine the memory allocation for gadget gmktimesteps : set output timesteps for cosmological simulations gmktimesteps -p 10 > outputtimesteps.txt gread_header : return some info on the header (incomplete, actually return the atime) gread_header snap.dat gheader : return some info on the header gheader snap.dat addgas : add gas to a cosmological box ################################################ # global statitics on the run gcpu : plot cpu info gcpu cpu.txt gcpu -o budget --legend cpu.txt gcpu -o gravity --legend cpu.txt gcpu -o hydro --legend cpu.txt genergy : plot info relative to the energy genergy --relative energy.txt genergy -o EnergyKineticFeedback energy.txt gplot_energy : plot info relative to the energy budget gplot_energy --legend energy.txt gsteps : plot info relative to the steps gsteps info.txt ################################################ pscripts/ ################################################ ################################################ # global output analysis # star formation rate and stellar mass (sfr.txt) pSfrvsTime --param params -o SFR sfr.txt pSfrvsTime --param params -o MSTARS sfr.txt pSfrvsTime --param params -o MNEWSTARS sfr.txt pSfrvsTime --param params -o sfr sfr.txt pSfrvsTime --param params -o mnewstars --integrate --derive --rf 100 sfr.txt pSfrvsTime --param params -o mstars sfr.txt # supernovae (chimie.txt) gSNvsTime -o NSNII,NSNIa --rf 1000 chimie.txt gSNvsTime -o EgySNThermal,EgySNKinetic --rf 1000 chimie.txt gSNvsTime -o EgySNThermal,EgySNKinetic --integrate chimie.txt gSNvsTime -o Nflag --rf 1000 chimie.txt gSNvsTime -o Nflag --integrate chimie.txt ################################################ # snapshot analysis # logRho - logT pplot --param params --select=gas --x logrho --y logT --xmin -3 --xmax 0 --ymin 1 --ymax 5 snap.dat pplot --param params --select=gas --x logrho --y logT --xmin -3 --xmax 0 --ymin 1 --ymax 5 --map snap.dat pplot --param params --select=gas --colorbar --log z --x logrho --y logT --z Tcool --xmin -3 --xmax 0 --ymin 1 --ymax 5 --zmin 1 --zmax 1e16 snap.dat # rho - T pplot --param params --select=gas --x rho --y T --log xy --xmin 1e-3 --xmax 1e0 --ymin 1e1 --ymax 1e5 snap.dat # R - T pplot --param params --select=gas --x R --y logT snap.dat # R - Fe pplot --param params --select=gas --x R --y Fe snap.dat # MgFe - Fe pplot --param params --select=stars1 --x Fe --y MgFe --xmin=-4 --xmax=0.5 --ymin=-1 --ymax=1.5 snap.dat pplot --param params --select=stars1 --x Fe --y MgFe --xmin=-4 --xmax=0.5 --ymin=-1 --ymax=1.5 --map snap.dat # Fe - Age pplot --param params --select=stars1 --x Age --y Fe --xmin=0 --xmax=14 --ymin=-4 --ymax=0.5 snap.dat # Tcool pplot --param params --select=gas --x logrho --y Tcool --xmin=-5 --xmax=2 --ymin 0 --ymax 1e16 --log y snap.dat ################################ # metallicite ################################ # MgFevsFe pMgFevsFe --select=stars1 --xmin=-4 --xmax=0.5 --ymin=-1 --ymax=1.5 --map snap.dat # NvsAge pNvsX --param params --select stars1 --x Age --nx 14 --xmin=0 --xmax=14 snap.dat # NvsFe pNvsX --param params --select stars1 --x Fe --nx 40 --xmin=-4 --xmax=0.5 --ymin=0 --ymax=0.15 snap.dat ################################ # profiles ################################ # grille spherique # density pSphericalProfile -o density --xmax 10 --rmax 10 --nr 64 --center histocenter --log xy --param params snap.dat # masse pSphericalProfile -o mass --xmax 10 --rmax 10 --nr 64 --center histocenter --log xy --param params snap.dat # integrated mass pSphericalProfile -o imass --xmax 10 --rmax 10 --nr 64 --center histocenter --log xy --param params snap.dat # velocity disp. in z pSphericalProfile -o sigmaz --xmax 10 --rmax 10 --nr 64 --center histocenter --param params snap.dat # circular velocity pSphericalProfile -o vcirc --xmax 10 --rmax 10 --nr 64 --center histocenter --param params --eps 0.1 snap.dat # grille cylindrique # density pCylindricalProfile -o sdens --xmax 10 --rmax 10 --nr 64 --center histocenter --log xy --param params snap.dat # masse pCylindricalProfile -o mass --xmax 10 --rmax 10 --nr 64 --center histocenter --log xy --param params snap.dat # integrated mass pCylindricalProfile -o imass --xmax 10 --rmax 10 --nr 64 --center histocenter --log xy --param params snap.dat # velocity disp. in z pCylindricalProfile -o sigmaz --xmax 10 --rmax 10 --nr 64 --center histocenter --param params snap.dat # circular velocity pCylindricalProfile -o vcirc --xmax 10 --rmax 10 --nr 64 --center histocenter --param params --eps 0.1 snap.dat - + +# luminosity profile +pCylindricalProfile -o Lum --xmax 10 --rmax 10 --nr 64 --center histocenter --param params snap.dat""" ############################################################ # velocities and velocity dispertions in cylindrical coords ############################################################ # extract all gExtractCylindrical --Rmax 10 --components="('gas','stars1','halo1')" --disk="('gas','stars1')" --nR 32 --eps=0.01 --nmin 3 --param params snap.dat # number of points per bins pplotdmp --mode=Numbers gas.dmp # circular velocity per component pplotdmp --mode=Velocities --legend gas.dmp stars1.dmp halo1.dmp # velocities and velocity disp pplotdmp --mode=Dispersions --legend gas.dmp # stability Toomre paramter Q pplotdmp --mode=Q gas.dmp """ diff --git a/gscripts/gmktimesteps b/gscripts/gmktimesteps index ebd0041..d0feec4 100755 --- a/gscripts/gmktimesteps +++ b/gscripts/gmktimesteps @@ -1,74 +1,74 @@ -#!/usr/bin/python +#!/usr/bin/env python ''' Create a liste of cosmological timesteps in the form of parameter a. Yves Revaz ven avr 14 12:35:53 CEST 2006 ''' from optparse import OptionParser import sys def parse_options(): usage = "usage: %prog [options]" parser = OptionParser(usage=usage) parser.add_option("-p", action="store", dest="p", type='int', default=5, metavar=" INT", help="precision value") parser.add_option("-v", action="store_true", dest="verbose", default = 0, help="verbose mode") (options, args) = parser.parse_args() files = args return files,options ########################################################## # # MAIN # ######################################################### # get options files,options = parse_options() p = options.p verbose = options.verbose da = 1./(2.**p) a = 0 while 1: a = a + da z = 1./a - 1.0 if z < 0: break if verbose: print "a = %12.10f z = %12.8f"%(a,z) else: print "%12.10f"%(a) diff --git a/pscripts/pCylindricalProfile b/pscripts/pCylindricalProfile index e0a43ce..529626d 100755 --- a/pscripts/pCylindricalProfile +++ b/pscripts/pCylindricalProfile @@ -1,426 +1,463 @@ #!/usr/bin/env python from optparse import OptionParser import Ptools as pt from pNbody import * from pNbody import units from pNbody import ctes from pNbody import thermodyn import string from scipy import optimize def parse_options(): usage = "usage: %prog [options] file" parser = OptionParser(usage=usage) parser = pt.add_postscript_options(parser) parser = pt.add_ftype_options(parser) parser = pt.add_reduc_options(parser) parser = pt.add_center_options(parser) parser = pt.add_select_options(parser) parser = pt.add_cmd_options(parser) parser = pt.add_display_options(parser) parser = pt.add_info_options(parser) parser = pt.add_limits_options(parser) parser = pt.add_log_options(parser) parse = pt.add_units_options(parser) parser.add_option("--rmax", action="store", dest="rmax", type="float", default = 50., help="max radius of bins", metavar=" FLOAT") parser.add_option("--nr", action="store", dest="nr", type="int", default = 32, help="number of bins in r", metavar=" INT") parser.add_option("--forceComovingIntegrationOn", action="store_true", dest="forceComovingIntegrationOn", default = False, help="force the model to be in in comoving integration") parser.add_option("--curve", action="store_true", dest="curve", default = False, help="draw curve instead of histogram") parser.add_option("-o", action="store", type="string", dest="o", default = 'density', help="value to plot") parser.add_option("--AdaptativeSoftenning", action="store_true", dest="AdaptativeSoftenning", default = False, help="AdaptativeSoftenning") parser.add_option("--ErrTolTheta", action="store", dest="ErrTolTheta", type="float", default = 0.5, help="Error tolerance theta", metavar=" FLOAT") parser.add_option("--eps", action="store", dest="eps", type="float", default = 0.1, help="smoothing length", metavar=" FLOAT") parser.add_option("--nt", action="store", dest="nt", type="int", default = 32, help="number of division in theta", metavar=" INT") (options, args) = parser.parse_args() pt.check_files_number(args) files = args return files,options ####################################### # MakePlot ####################################### def MakePlot(files,opt): # some inits colors = pt.Colors(n=len(files)) datas = [] # read files for file in files: nb = Nbody(file,ftype=opt.ftype) ################ # units ################ # define local units unit_params = pt.do_units_options(opt) nb.set_local_system_of_units(params=unit_params) # define output units # nb.ToPhysicalUnits() if opt.forceComovingIntegrationOn: nb.setComovingIntegrationOn() ################ # apply options ################ nb = pt.do_reduc_options(nb,opt) nb = pt.do_select_options(nb,opt) nb = pt.do_center_options(nb,opt) nb = pt.do_cmd_options(nb,opt) nb = pt.do_info_options(nb,opt) nb = pt.do_display_options(nb,opt) ################ # some info ################ print "---------------------------------------------------------" nb.localsystem_of_units.info() nb.ComovingIntegrationInfo() print "---------------------------------------------------------" # grid division rc = 1 f = lambda r:log(r/rc+1.) fm = lambda r:rc*(exp(r)-1.) ############################### # compute physical quantities ############################### ########################################## # Cylindrical_2drt_Grid ########################################## if opt.o=='sdens': G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm) x,t = G.get_rt() y = G.get_SurfaceDensityMap(nb) # comoving conversion if nb.isComovingIntegrationOn(): print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam) x = x*nb.atime/nb.hubbleparam # length conversion y = y/nb.atime**2*nb.hubbleparam # surface density conversion # output units out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) out_units_y = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength) fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitSurfaceDensity) x = x*fx y = y*fy # set labels xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$' ylabel = r'$\rm{Surface\,\,Density}\,\left[ M_{\odot}/kpc^2 \right]$' x,y = pt.CleanVectorsForLogX(x,y) x,y = pt.CleanVectorsForLogY(x,y) datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points')) if opt.o=='mass': G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm) x,t = G.get_rt() y = G.get_MassMap(nb) y = sum(y,axis=1) # comoving conversion if nb.isComovingIntegrationOn(): print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam) x = x*nb.atime/nb.hubbleparam # length conversion y = y/nb.hubbleparam # mass conversion # output units out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) out_units_y = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength) fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitDensity) x = x*fx y = y*fy # set labels xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$' ylabel = r'$\rm{Mass}\,\left[ M_{\odot} \right]$' x,y = pt.CleanVectorsForLogX(x,y) x,y = pt.CleanVectorsForLogY(x,y) datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points')) if opt.o=='imass': G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm) x,t = G.get_rt() y = G.get_MassMap(nb) y = sum(y,axis=1) y = add.accumulate(y) # comoving conversion if nb.isComovingIntegrationOn(): print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam) x = x*nb.atime/nb.hubbleparam # length conversion y = y/nb.hubbleparam # mass conversion # output units out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) out_units_y = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength) fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitDensity) x = x*fx y = y*fy # set labels xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$' ylabel = r'$\rm{Mass}\,\left[ M_{\odot} \right]$' x,y = pt.CleanVectorsForLogX(x,y) x,y = pt.CleanVectorsForLogY(x,y) datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points')) if opt.o=='sigmaz': G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm) x,t = G.get_rt() y = G.get_SigmaValMap(nb,nb.Vz()) y = sum(y,axis=1) # comoving conversion if nb.isComovingIntegrationOn(): print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam) x = x*nb.atime/nb.hubbleparam # length conversion #y = y/nb.hubbleparam # velocity print "!!! this is not defined for comobiles coords, please, do it !!!" sys.exit() # output units out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) out_units_y = units.UnitSystem('local',[units.Unit_km,units.Unit_Ms,units.Unit_s,units.Unit_K]) fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength) fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitVelocity) x = x*fx y = y*fy # set labels xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$' ylabel = r'$\sigma_z\,\left[ \rm{km}/\rm{s} \right]$' x,y = pt.CleanVectorsForLogX(x,y) x,y = pt.CleanVectorsForLogY(x,y) datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points')) if opt.o=='vcirc': G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=opt.nt,z=0,g=f,gm=fm) x,t = G.get_rt() # compute tree nb.getTree(force_computation=True,ErrTolTheta=opt.ErrTolTheta) # compute accel Accx,Accy,Accz = G.get_AccelerationMap(nb,eps=opt.eps,UseTree=True,AdaptativeSoftenning=opt.AdaptativeSoftenning) Ar = sqrt(Accx**2+Accy**2) Ar = sum(Ar,1)/opt.nt dPhit = Ar y = libdisk.Vcirc(x,dPhit) # comoving conversion if nb.isComovingIntegrationOn(): print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam) x = x*nb.atime/nb.hubbleparam # length conversion #y = y/nb.hubbleparam # velocity print "!!! this is not defined for comobiles coords, please, do it !!!" sys.exit() # output units out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) out_units_y = units.UnitSystem('local',[units.Unit_km,units.Unit_Ms,units.Unit_s,units.Unit_K]) fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength) fy = nb.localsystem_of_units.convertionFactorTo(out_units_y.UnitVelocity) x = x*fx y = y*fy # set labels xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$' ylabel = r'$\rm{Circular\,Velocity}\,\left[ \rm{km}/\rm{s} \right]$' x,y = pt.CleanVectorsForLogX(x,y) x,y = pt.CleanVectorsForLogY(x,y) datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points')) + if opt.o=='Lum': + + G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm) + + x,t = G.get_rt() + y = reshape(G.get_SurfaceValMap(nb,nb.luminosity_spec()*1e10),opt.nr) # in Lsol/(unitlength) + + + # comoving conversion + if nb.isComovingIntegrationOn(): + print " converting to physical units (a=%5.3f h=%5.3f)"%(nb.atime,nb.hubbleparam) + x = x*nb.atime/nb.hubbleparam # length conversion + y = y/nb.atime**2*nb.hubbleparam # surface density conversion + + + # output units + out_units_x = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) + out_units_y = units.UnitSystem('local',[units.Unit_pc,units.Unit_Ms,units.Unit_Myr,units.Unit_K]) + + fx = nb.localsystem_of_units.convertionFactorTo(out_units_x.UnitLength) + fy = nb.localsystem_of_units.convertionFactorTo(1/out_units_y.UnitSurface) # do not count L + + x = x*fx + y = y*fy + + + # set labels + + xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$' + ylabel = r'$\rm{Luminosity}\left[ L_{\odot}/\rm{pc}^2 \right]$' + + if opt.log: + x,y = pt.CleanVectorsForLogX(x,y) + x,y = pt.CleanVectorsForLogY(x,y) + + datas.append(pt.DataPoints(x, y, color=colors.get(),label=None,tpe='points')) + ################## # plot all ################## for d in datas: pt.plot(d.x,d.y,color=d.color) # set limits and draw axis xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log) pt.SetAxis(xmin,xmax,ymin,ymax,opt.log) pt.xlabel(xlabel,fontsize=pt.labelfont) pt.ylabel(ylabel,fontsize=pt.labelfont) pt.grid(False) ######################################################################## # MAIN ######################################################################## if __name__ == '__main__': files,opt = parse_options() pt.InitPlot(files,opt) pt.pcolors MakePlot(files,opt) pt.EndPlot(files,opt)