diff --git a/Gtools/pscripts/pSFRvsTime b/Gtools/pscripts/pSFRvsTime
new file mode 100755
index 0000000..85162ac
--- /dev/null
+++ b/Gtools/pscripts/pSFRvsTime
@@ -0,0 +1,298 @@
+#!/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("-o",
+ action="store",
+ dest="obs",
+ type="string",
+ default = 'SFR',
+ help="observable name",
+ metavar=" NAME")
+
+ parser.add_option("--x",
+ action="store",
+ dest="x",
+ type="string",
+ default = 'r',
+ help="x value to plot",
+ metavar=" STRING")
+
+ parser.add_option("--nx",
+ action="store",
+ type="int",
+ dest="nx",
+ default = 500,
+ help="number of bins")
+
+
+ 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")
+
+
+ (options, args) = parser.parse_args()
+
+
+ pt.check_files_number(args)
+
+ files = args
+
+ return files,options
+
+
+
+def get_value_and_label(nb,val):
+
+
+
+ if val == 'Age':
+ label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
+
+ out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
+
+ val = nb.StellarAge(units=out_units.UnitTime)
+
+ return val,label
+
+ elif val == 'CosmicTime':
+ label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
+
+ out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
+
+ val = nb.CosmicTime(units=out_units.UnitTime)
+
+ return val,label
+
+
+ elif val == 'Redshift':
+ label = r'$\rm{Redshift}$'
+
+ out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
+
+ val = nb.Redshift()
+ return val,label
+
+ elif val == 'ScalingFactor':
+ label = r'$\rm{a}$'
+
+ val = nb.tstar
+
+ return val,label
+
+
+
+ else:
+ label = r'%s'%val
+ print "val = %s"%val
+ exec("val = %s"%val)
+ return val, label
+
+
+
+
+#######################################
+# MakePlot
+#######################################
+
+
+def MakePlot(files,opt):
+
+
+ # some inits
+ colors = pt.Colors(n=len(files))
+
+
+ # 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 = nb.select(1)
+
+ nb = pt.do_reduc_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 "---------------------------------------------------------"
+
+
+
+ if opt.obs=='SFR':
+ opt.x = "CosmicTime"
+ opt.ylabel = r'$\rm{Sfr}\,\left[ M_\odot/yr \right]$'
+ opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
+
+ if opt.obs=='MSTARS':
+ opt.x = "CosmicTime"
+ opt.ylabel = r'$\rm{Stellar\,\,Mass}\,\left[ M_\odot \right]$'
+ opt.xlabel = r'$\rm{Time}\,\left[ Gyr \right]$'
+
+ if opt.obs=='MSTARSz':
+ opt.x = "Redshift"
+ opt.ylabel = r'$\rm{Stellar\,\,Mass}\,\left[ M_\odot \right]$'
+ opt.xlabel = r'$\rm{Redshift}$'
+
+ if opt.obs=='MSTARSa':
+ opt.x = "ScalingFactor"
+ opt.ylabel = r'$\rm{Stellar\,\,Mass}\,\left[ M_\odot \right]$'
+ opt.xlabel = r'$\rm{a}$'
+
+
+
+
+
+ ################
+ # get values
+ ################
+
+ x,xlabel = get_value_and_label(nb,opt.x)
+
+ # do the 1d histogram
+ G = libgrid.Generic_1d_Grid(opt.xmin,opt.xmax,opt.nx)
+
+ y = G.get_MassMap(x,nb.mass) # mass weighted histogram
+ x = G.get_r()
+
+
+ if opt.obs=='SFR':
+
+ y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
+ y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
+ dt = (x[1:]-x[:-1])*1e9 # Gyrs to yrs
+ dMdt = y[:-1] / dt
+
+ y = dMdt
+ x = x[1:]
+
+
+ if opt.obs=='MSTARS':
+
+ y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
+ y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
+
+ y = add.accumulate(y)
+
+
+ if opt.obs=='MSTARSz':
+ y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
+ y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
+
+ y = y[::-1]
+ y = add.accumulate(y)
+ y = y[::-1]
+
+ if opt.obs=='MSTARSa':
+ y_out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_yr,units.Unit_K])
+ y = y * nb.localsystem_of_units.convertionFactorTo(units.Unit_Ms)
+ y = add.accumulate(y)
+
+
+
+
+ pt.plot(x,y,color=colors.get())
+
+
+
+ if file == files[0]:
+ xmin,xmax,ymin,ymax = pt.SetLimits(opt.xmin,opt.xmax,opt.ymin,opt.ymax,x,y,opt.log)
+
+ # plot axis
+ pt.SetAxis(xmin,xmax,ymin,ymax,opt.log)
+ pt.xlabel(opt.xlabel,fontsize=pt.labelfont)
+ pt.ylabel(opt.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)
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Gtools/pscripts/pSFRvsTime.py b/Gtools/pscripts/pSFRvsTime.py
new file mode 120000
index 0000000..17b22c2
--- /dev/null
+++ b/Gtools/pscripts/pSFRvsTime.py
@@ -0,0 +1 @@
+pSFRvsTime
\ No newline at end of file
diff --git a/Gtools/pscripts/pget_Luminosity b/Gtools/pscripts/pget_Luminosity
new file mode 100755
index 0000000..5a92fce
--- /dev/null
+++ b/Gtools/pscripts/pget_Luminosity
@@ -0,0 +1,95 @@
+#!/usr/bin/python
+
+
+from optparse import OptionParser
+import Ptools as pt
+from pNbody import *
+from pNbody import units
+from pNbody import ctes
+import string
+
+from pNbody import units,io,ctes
+from scipy.optimize import bisect as bisection
+
+def parse_options():
+
+
+ usage = "usage: %prog [options] file"
+ parser = OptionParser(usage=usage)
+
+
+ parse = pt.add_units_options(parser)
+ parser = pt.add_ftype_options(parser)
+
+ (options, args) = parser.parse_args()
+
+
+ files = args
+
+ return files,options
+
+
+
+
+########################################################################
+# MAIN
+########################################################################
+
+if __name__ == '__main__':
+ files,opt = parse_options()
+
+
+
+ if len(files)==0:
+ sys.exit()
+
+ file = files[0]
+
+ nb = Nbody(file,ftype=opt.ftype)
+
+ unit_params = pt.do_units_options(opt)
+ HubbleParam = unit_params['HubbleParam']
+ nb.set_local_system_of_units(params=unit_params)
+
+ nb = nb.select(1)
+
+ Rmax = nb.rxyz().max()
+
+
+ ############################
+ # find tR using bissectrice
+ ############################
+ lmax=90
+
+ nb.L = nb.Luminosity()
+
+ Ltot = sum(nb.L)
+ rs = nb.rxyz()
+
+ def getL(r):
+ nb_s = nb.selectc(rs<r)
+ L = sum(nb_s.L)
+ return lmax - 100*L/Ltot
+
+ Rt = bisection(getL, 0.0, 50, args = (), xtol = 0.1, maxiter = 400)
+
+
+
+ nb_s = nb.selectc(rs<Rt)
+ L = sum(nb_s.L)
+
+
+ out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Msol,units.Unit_s,units.Unit_K])
+ fL = nb.localsystem_of_units.convertionFactorTo(out_units.UnitLength)
+
+ print
+ print "tidal radius : Rt = %g [kpc/h]"%(Rt*fL)
+ print "tidal radius : Rt = %g [kpc]"%(Rt*fL/HubbleParam)
+ print "Luminosity in tidal radius : L = %g [1e6 Lsol]"%(L/1e6)
+ print "Total luminosity : Ltot = %g [1e6 Lsol]"%(Ltot/1e6)
+
+
+
+
+
+
diff --git a/Gtools/pscripts/pget_VirialRadius b/Gtools/pscripts/pget_VirialRadius
new file mode 100755
index 0000000..f9009af
--- /dev/null
+++ b/Gtools/pscripts/pget_VirialRadius
@@ -0,0 +1,132 @@
+#!/usr/bin/python
+
+
+from optparse import OptionParser
+import Ptools as pt
+from pNbody import *
+from pNbody import units
+from pNbody import ctes
+import string
+
+from pNbody import units,io,ctes
+from scipy.optimize import bisect as bisection
+
+def parse_options():
+
+
+ usage = "usage: %prog [options] file"
+ parser = OptionParser(usage=usage)
+
+ parser.add_option("--dX",
+ action="store",
+ dest="dX",
+ type="float",
+ default = 200,
+ help="over density",
+ metavar=" STRING")
+
+ parse = pt.add_units_options(parser)
+ parser = pt.add_ftype_options(parser)
+
+ (options, args) = parser.parse_args()
+
+
+ files = args
+
+ return files,options
+
+
+
+
+########################################################################
+# MAIN
+########################################################################
+
+if __name__ == '__main__':
+ files,opt = parse_options()
+
+ X = opt.dX
+
+ # define local units
+ unit_params = pt.do_units_options(opt)
+ system_of_units = units.Set_SystemUnits_From_Params(unit_params)
+
+
+ G=ctes.GRAVITY.into(system_of_units)
+ H = ctes.HUBBLE.into(system_of_units)
+ HubbleParam = unit_params['HubbleParam']
+
+
+ rhoc = pow(H,2)*3/(8*pi*G)
+ rhoX = rhoc*X * unit_params['Omega0']
+
+ print "rhoX (code units, dX=%g)"%opt.dX,rhoX
+
+
+
+ # output system of units (the mass units is the hydrogen mass)
+ Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
+ Unit_atom.set_symbol('atom')
+ out_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_s,units.Unit_K])
+
+ funit = system_of_units.convertionFactorTo(out_units.UnitDensity)
+
+ atime = 1.0
+ print " converting to physical units (a=%5.3f h=%5.3f)"%(atime,HubbleParam)
+ rhoXu = rhoX/atime**3*HubbleParam**2 *funit
+ print "rhoX (atom/cm^3)",rhoXu
+ print "log10rhoX (atom/cm^3)",log10(rhoXu)
+
+
+ ###############################
+ # now, deal with files
+ ###############################
+
+ if len(files)==0:
+ sys.exit()
+
+ file = files[0]
+
+ nb = Nbody(file,ftype=opt.ftype)
+
+ Rmax = nb.rxyz().max()
+
+
+ ############################
+ # find rX using bissectrice
+ ############################
+ rs = nb.rxyz()
+
+ def getRes(r):
+
+
+ nb_s = nb.selectc(rs<r)
+ M = sum(nb_s.mass)
+ V = 4/3.*pi*r**3
+
+ return M/V - rhoX
+
+
+ rX = bisection(getRes, 0.1, Rmax, args = (), xtol = 0.001, maxiter = 400)
+
+ nb_s = nb.selectc(nb.rxyz()<rX)
+ MX = sum(nb_s.mass)
+ V = 4/3.*pi*rX**3
+
+
+ out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Msol,units.Unit_s,units.Unit_K])
+ fL = system_of_units.convertionFactorTo(out_units.UnitLength)
+ fM = system_of_units.convertionFactorTo(out_units.UnitMass)
+
+ print
+ print "Virial radius : r%d = %g [kpc/h]"%(X,rX*fL)
+ print "Virial mass : M%d = %g [Msol/h]"%(X,MX*fM)
+
+ print
+ print "Virial radius : r%d = %g [kpc]"%(X,rX*fL/HubbleParam)
+ print "Virial mass : M%d = %g [Msol]"%(X,MX*fM/HubbleParam)
+
+
+
+
+