gall        #----graphs/prints all
		define start 2093
		define end 2142
		
		do z=$start,$end, 1{
         	  clear
		  #hpcolor
		  #dev postfile "~/find_density/avg_$!num.ps"
		  graphden $z -2.5 0.0
		  x11
		  #cursor 
		  
		}
		clear
		#hpcolor
		#dev postfile "~/find_density/facto_vs_$!num.ps"
		graphana $start $end 0.0
		#x11
		cursor
				 
		clear
		#hpcolor
		#dev postfile "~/find_density/mass_lost_$!num.ps"
		gmloss $start $end
		#x11
		cursor
		
gden 1	graphden $1 -2.5 0.0
gdenr 1	graphden $1 -2.5 -0.5
		
graphden 3  #----graphs density: folder -2.5 0.0
		black
		expand 0.50000
		define num $1
		define rrmin $2
		define rrmax $3
		  data "~/find_density/$!num/avg_$!num.dat"
              read {i 1 iii 3 iv 4 v 5 vi 6 vii 7 viii 8 ix 9 x 10 xi 11}
		  data "~/find_density/$!num/summary_$!num.dat"
		  read {r200 5 Mhalo 6 facto 7}
		  data "~/find_density/$!num/navavg_$!num.dat"
		  read {rr200 1 xaxis 2 yaxis 3 yfit 4}
		  data "~/find_density/$!num/navarrofit_$!num.dat"
		  read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 7}
		  data "~/find_density/$!num/energy_avg_$!num.dat"
		  read {err200 2 dm 4 ke 5 pe 6 ii 7}
		if ($(conx1[0])>$(iii[$(dimen(iii)-1)])){set conx1=1000}
		if ($(conx2[0])>$(iii[$(dimen(iii)-1)])){set conx2=1000}
		if ($(conx3[0])>$(iii[$(dimen(iii)-1)])){set conx3=1000}		
		
		set vrms=lg(((10**vi)**2+(10**v)**2)**0.5)
		set k=0
		do j=0, dimen(iii)-1,1{
              if((iii[$j]>$rrmin) && (iii[$j]<$rrmax)){
		    set k=1+$(k[0])}} 
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
 		set dimen(z)=$(k[0])
		set dimen(viii2)=$(k[0])
		set dimen(vrms2)=$(k[0])
		set k=0
		do j=0, dimen(iii)-1,1{
                 if( (iii[$j]>$rrmin) && (iii[$j]<$rrmax) ){
	             set s[$(k[0])]=$(iii[$j])
		       set w[$(k[0])]=$(iv[$j])
                   set z[$(k[0])]=$(vii[$j])
		       set viii2[$(k[0])]=$(viii[$j])
		       set vrms2[$(k[0])]=$(vrms[$j])
                   set k=1+$(k[0])}}
		set dimen(concen)=2
		set dimen(1concen)=2
		set dimen(2concen)=2
		set dimen(3concen)=2
		set concen[0]=-6.0
		set concen[1]=12.0
		set 1concen[0]=$(conx1[0])
		set 2concen[0]=$(conx2[0])
		set 3concen[0]=$(conx3[0])
		set 1concen[1]=$(conx1[0])
		set 2concen[1]=$(conx2[0])
		set 3concen[1]=$(conx3[0])
		set dimen(xvect)=2
		set dimen(y1vect)=2
		set dimen(y2vect)=2
		set dimen(y3vect)=2
		set xvect[0]=-4.5
		set xvect[1]=2.0
		set y1vect[0]=-1.0
		set y1vect[1]=-1.0
		set y2vect[0]=-2.0
		set y2vect[1]=-2.0
		set y3vect[0]=-3.0
		set y3vect[1]=-3.0
		  
		#DENSITY/rho_o
		set iv=iv+63-LG(2.3563)
		set xi=xi+63-LG(2.3563)
		set w=w+63-LG(2.3563)
		#divide by rho_0 in M_sol/cm^3
		hex1
		xlabel log(r/r_{200})  
            ylabel log [(\rho/\rho_o)*(r/r_{200})^2]
		limits -4.5 2.0 -3.5 3.5
		box
		ptype 1 1
            #con iii (xi+2*iii)
            points iii (iv+2*iii)
	 	#con iii (iv+2*iii)
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		ptype 4 1
		lsq s w s q rms 
		con s (q+2*s+.1)
		move 5 12
		label \frac{\rho}{\rho_0}=$(sprintf('%5.3f',$a))*log\frac{ r }{r_{200}}+$(sprintf('%5.3f',$b))
		move 5 5
		label rms: $(sprintf('%5.3f',$rms))
		set beta=$((-1)*($a))
		set betarms=$rms
		
		#V RMS
		#change units to M_sol*km^3/(mpc^3*s^3)
		set vrms=vrms-lg(1E5)
		set vrms2=vrms2-lg(1E5)
		hex2
		xlabel log(r/r_{200})  
		ylabel log(V_{rms} [\frac{km}{s}] )
	 	limits -4.5 2.0 1 3.0
            box 
		ptype 1 1
            points iii vrms	
		lsq s vrms2 s q rms 
		con s q
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		move 5 80
		label y=$(sprintf('%5.3f',$a))*x+$(sprintf('%5.3f',$b))
		move 5 75
		label rms: $(sprintf('%5.3f',$rms))
		set vela=$a
		set velrms=$rms
		
		#beta vs rad
		hex5
		xlabel log(r/r_{200})  
            ylabel \beta_{vel.disp.}
		limits -4.5 2.0 -.5 1.1
		box
		ptype 1 1
		points iii ix
		set dimen(betax)=2
		set dimen(betay)=2
		set betax[0]=-4.5
		set betax[1]=2.0
		set betay[0]=0
		set betay[1]=0
		ltype 2
		con betax betay		
		ltype 0
				
		#DEN/VRMS^3
		set vii=vii+LG(2.9377E38)+LG(1E20)
		set z=z+LG(2.9377E38)+LG(1E20)
		#change units to M_sol*km^3/(mpc^3*s^3)
		hex3
		ylabel log [(\rho/\sigma^3)*(r/r_{200})^{1.875}]
		limits -4.5 2.0 -30 -18
            box
		ptype 1 1
		points iii (vii+1.875*iii)
		#con iii (vii+1.875*iii)
		lsq s z s q rms 
		con s (q+1.875*s)
		move 5 15
		label log(\rho/\sigma^3)=($(sprintf('%5.3f',$a)))log(r/r_{200})+($(sprintf('%5.1f',$b)))
		move 5 10
		label rms: $(sprintf('%5.3f',$rms))		
		set alpha=$((-1)*($a))
		set alpharms=$rms	
		
		#navarro, fig4
		hex6		
		xlabel log (\sigma^3/\rho) (\frac{Mpc^3s^3}{M_{sol}km^3})
            ylabel log (dM/dlog(\sigma^3/\rho))
		limits 50 200 -30 18.0
		box
		ptype 1 1
            points xaxis yaxis
		#con xaxis yaxis
		ptype 1 1
		points xaxis yfit
		con xaxis yfit
            move 5 15
		label y=($(sprintf('%5.3f',$(aout[0]))))x^2+($(sprintf('%5.3f',$(bout[0]))))x+($(sprintf('%5.3f',$(cout[0]))))
		move 5 10
		label chisq=$(sprintf('%5.3f',$(chisq[0])))
		set gamm=$(aout[0])
		set gammarms=$(chisq[0])	
		
		##slope of den v. r200
		#location 3500 11000 9500 15500
		#ylabel \rho v. r/r_{200} slope
		#limits -4.5 2.0 -5 0
            #box 0 2
		#ltype 1
		#points iii x
		#con 1concen concen
		#con xvect y1vect
		#con 2concen concen
		#con xvect y2vect
		#con 3concen concen
		#con xvect y3vect
		#ltype 0
		#rperi/rapo vs. r/r200
		#location 3500 11000 3500 9500 
		#xlabel log(r/r_{200})  
            #ylabel rperi/rapo
		#limits -4.5 2.0 0.0 0.5
            #box
		#points iii viii
		#con s q
		#move 5 13
		#label y=$(sprintf('%5.3f',$a))*x+$(sprintf('%5.3f',$b))
		#move 5 5
		#label rms: $(sprintf('%5.3f',$rms))
		
		#energy graph
		hex4
		xlabel log(r/r_{200})  
		ylabel log(energy per mass [cm^2/s^2] )
		limits -4.5 2.0 0 22
		box
		ptype 1 1
	     ctype dred
	      points err200 (ke/dm)
		#con err200 ke
		move 10 90
		label red: ke/dm
	     ctype dblue
		#con err200 pe
		points err200 (pe/dm)
		move 10 85
		label blue: pe/dm
	     ctype grey
		#con err200 dm
		points err200 dm
		move 10 80
		label grey: dm
	     ctype cyan
		con err200 ii
		move 10 75
		label cyan: potential from code
	 		
	      ctype black
		toplabel FILE: $num                 FACTO: $(sprintf('%5.2f',$(facto[0])))			
		if ($rrmax==0.0){
		write "~/find_density/.temp_fits_$!num.dat" $(facto[0]) $(beta[0]) $(betarms[0]) $(alpha[0]) $(alpharms[0]) $(gamm[0])  $(gammarms[0])   $(vela[0])    $(velrms[0])
		}
		if ($rrmax==-0.5){
		write "~/find_density/.temp_fits_$!num.datr" $(facto[0]) $(beta[0]) $(betarms[0]) $(alpha[0]) $(alpharms[0]) $(gamm[0])  $(gammarms[0])   $(vela[0])    $(velrms[0])
		}
		
		DEL i ii iii iv v
		DEL r200 Mhalo facto shell integ 
		DEL aout baout cout chisq vrms 
		DEL vi vii viii ix ten
		DEL eccen rr200 xaxis yaxis yfit
		DEL u s w x y 
		DEL z vii2 ten2 vrms2 k
		DEL gamm gammarms vela velrms x
		DEL alpha alpharms beta betarms concen
		DEL conx1 conx2 conx3 1concen 2concen
		DEL 3concen q rrmax rrmin viii2
		DEL bout xi xvect y1vect y2vect 
		DEL y3vect dm err200 ke pe
		DEL betax betay x x x
		echo ----
		list set
		echo ----

gana 2	graphana $1 $2 0.0
ganar 2	graphana $1 $2 -0.5

graphana 3	#----num num2, graphs analysis of files from num to num2
		define num $1
		define num2 $2
		define rrmax $3
		black
		expand 0.500000
		erase
		do temp=$num, $num2,1{
              if ($rrmax==0.0){
		    data "~/find_density/.temp_fits_$!temp.dat"
		    read {facto 1 beta 2 betarms 3 alpha 4 alpharms 5 gamm 6  gammarms 7 vela 8 velrms 9}
		    write + "~/find_density/.temp_ana_$!num.dat" $(facto[0]) $(beta[0]) $(betarms[0]) $(alpha[0]) $(alpharms[0]) $(gamm[0]) $(gammarms[0])    $(vela[0])     $(velrms[0])
		    data "~/find_density/$!temp/navarrofit_$!temp.dat"
		    read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 7}
		    write + "~/find_density/.temp_info_$!num.dat" $(aout[0]) $(bout[0]) $(cout[0]) $(chisq[0]) $(conx1[0]) $(conx2[0]) $(conx3[0])
		    }
		  if ($rrmax==-0.5){
		    data "~/find_density/.temp_fits_$!temp.datr"
		    read {facto 1 beta 2 betarms 3 alpha 4 alpharms 5 gamm 6  gammarms 7 vela 8 velrms 9}
		    write + "~/find_density/.temp_ana_$!num.datr" $(facto[0]) $(beta[0]) $(betarms[0]) $(alpha[0]) $(alpharms[0]) $(gamm[0]) $(gammarms[0])    $(vela[0])     $(velrms[0])
		    data "~/find_density/$!temp/navarrofit_$!temp.dat"
		    read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 7}
		    write + "~/find_density/.temp_info_$!num.datr" $(aout[0]) $(bout[0]) $(cout[0]) $(chisq[0]) $(conx1[0]) $(conx2[0]) $(conx3[0])
		    }
		     
		  DEL facto beta betarms alpha x
		  DEL aout alpharms gamm bout conx3 
		  DEL gammarms vela velb velrms aout 
		  DEL bout cout chisq conx1 conx2 
		  DEL x x vela velrms x 
		  echo ----
		  list set
		  echo ----
		}
		if ($rrmax==0.0){
		data "~/find_density/.temp_ana_$!num.dat"
		read {facto 1 beta 2 betarms 3 alpha 4 alpharms 5 gamm 6 gammarms 7 vela 8 velrms 9}
		data "~/find_density/.temp_info_$!num.dat"
		read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 7}}
		if ($rrmax==-0.5){
		data "~/find_density/.temp_ana_$!num.datr"
		read {facto 1 beta 2 betarms 3 alpha 4 alpharms 5 gamm 6 gammarms 7 vela 8 velrms 9}
		data "~/find_density/.temp_info_$!num.datr"
		read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 7}}
		
		#DENSITY/rho_0
		hex1
		xlabel Scaling Parameter
            ylabel (log(\rho/\rho_o) vs. log (r/r200)) slope
		limits facto beta
            box
		ptype 5 1
            points facto beta
		ptype 1 1
		con facto beta
		lsq facto beta facto q rms 
		
		#V RMS
		hex2
		xlabel Scaling Parameter  
		ylabel (log V rms vs. log r/r200) slope
	 	limits facto vela
            box
		ptype 5 1
            points facto vela	
		con facto vela			
            ptype 1 1
		lsq facto vela facto q rms 
		
		#DEN/VRMS^3
		hex3
		xlabel Scaling Parameter
		ylabel (log \rho/\sigma^3 vs. log r/r_{200}) slope
		limits facto alpha
            box
		set dimen(nfw)=2
		set dimen(nfwy)=2
		set nfw[0]=-1.875
		set nfw[1]=-1.875
		set nfwy[0]=0.6
		set nfwy[1]=1.5
		ptype 1 1
		points nfwy nfw
		ltype 3
		con nfwy nfw
		relocate 1.5 -1.875
		label NFW
		ltype 0
		ptype 5 1
            points facto alpha
		con facto alpha
		ptype 1 1
		lsq facto alpha facto q rms 
				
		#concentration factors
		hex4
		xlabel Scaling Parameter
		ylabel log(\frac{ r }{r_{200}}) of density-slope
		limits facto -5 5
		box
		ptype 4 3
		points facto conx1
		con facto conx1
		ptype 3 3
		points facto conx3
		con facto conx3
		ptype 5 2
		points facto conx2
		con facto conx2
		relocate .5 4.25
		label squares:-1 astrix:-2 triangles:-3
		
		#navarro
		hex5
		xlabel Scaling Parameter
            ylabel quadratic term of 'Navarro fig 4' fit
		limits facto gamm
            box
		ptype 5 1
            points facto gamm
		con facto gamm
            ptype 1 1
		lsq facto gamm facto q rms 
				
		#rms of den/rms^3 fits
		hex6
		xlabel Scaling Parameter
		ylabel \rho/\sigma^3 vs. log r/r_{200} RMS value
		limits facto alpharms
		box
		ptype 5 1
		black
		points facto alpharms
		con facto alpharms
		
		toplabel FACTO vs. averaged final quantities for folders $num to $num2
		DEL aout bout chisq conx1 conx2 
		DEL conx3 cout beta alpha betarms 
		DEL facto gammarms nfw alpharms gamm 
		DEL q vela velb velrms nfwy

		
		echo ----
		list set
		echo ----

gmloss 2	#----num num2: graphs halo mass loss
		define num $1
		define num2 $2
		black
	      expand 0.5
		erase

		do temp=$num, $num2,1{
              data "~/find_density/$!temp/summary_$!temp.dat"
		  read {numm 1 halo 2 shelllost 4 r200 5 Mhalo 6 facto 7 masslost 8 mfraclost 9}
		  set k=dimen(numm)
		  do j=0,k-1,1{
		    write + "~/find_density/.temp_mass_$!num.dat" $(numm[$j]) $(halo[$j]) $(shelllost[$j]) $(r200[$j]) $(Mhalo[$j]) $(facto[$j]) $(masslost[$j]) $(mfraclost[$j])
		    }}
		data "~/find_density/.temp_mass_$!num.dat"
		read {numm3 1 halo3 2 shelllost3 3 r2003 4 Mhalo3 5 facto3 6 masslost3 7 mfraclost3 8}
		
		#R 200
		quad1
		xlabel Scaling Parameter
		ylabel R_{200}
		limits facto3 r2003
		box
		ptype 6 1
		points facto3 r2003
		
		#massfrac lost
		quad2
		xlabel Scaling Parameter
            ylabel % of Mass Lost
		limits facto3 mfraclost3
		box
		ptype 6 1
		points facto3 mfraclost3
		
		#Mhalo and masslost
		quad3
		xlabel Scaling Parameter
		#relocate (2000 2000)
            ylabel Mass (M_{sol})
		set testlow=$(Mhalo3[0])
		set testhigh=$(Mhalo3[0])
		do z=0,dimen(Mhalo3)-1,1{
		  if ($(Mhalo3[$z])>$(testhigh[0])){
		    set testhigh=$(Mhalo3[$z])
		  }
		  if ($(Mhalo3[$z])<$(testlow[0])){
		    set testlow=$(Mhalo3[$z])
		  }
		  if ($(masslost3[$z])>$(testhigh[0])){
		    set testhigh=$(masslost3[$z])
		  }
		  if ($(masslost3[$z])<$(testlow[0])){
		    set testlow=$(masslost3[$z])
		  }
		}
		limits facto3 $(testlow[0]) $(testhigh[0])
		box
		ptype 6 1
		points facto3 Mhalo3
		ptype 4 1 
		points facto3 masslost3
		relocate 1.2 1.2E12
		label  red:final Mhalo 
		relocate 1.2 1.1E12 
		label blue:mass lost via formation
		
		#shells lost
		quad4
		xlabel Scaling Parameter
            ylabel \# of shells lost
		limits facto3 shelllost3
		box
		ptype 6 1
		points facto3 shelllost3
		
		DEL Mhalo Mhalo3 facto facto3 halo 
		DEL halo3 k masslost masslost3 mfraclost
		DEL mfraclost3 numm numm3 r200 r2003 
		DEL shelllost shelllost3 testhigh testlow	x
		echo ----
		list set
		echo ----