gallx       #----print all for paper
		
   	  #standard data
		do z=2093 , 2142, 1{		
		clear
		graphden 			$z	-2.5 	 0.0				
		clear
		graphden 			$z	-2.5 	-0.5			}			
		dev postlandfile "~/standard_vsfacto.ps"
		graphanax 			2093 	2142	 0.0			
		dev x11		
		graphanax 			2093	2142	-0.5								
		dev postlandfile "~/avg_2183.ps"
		graphdenx 			2183 	-2.5 	 0.0
		dev x11		
	  #figure 5 data				
		do z=3293 , 3342, 1{
		clear
		graphden 			$z 	-2.5 	 0.0			
		clear
		graphden 			$z	-2.5 	-0.5			}	
		clear
		graphanax 			3293 	3342	 0.0		
		clear
		graphanax 			3293 	3342	-0.5								
	  #rg data
		do z=3493 , 3542, 1{
		clear
		graphden 			$z 	-2.5 	 0.0		
		clear
		graphden 			$z	-2.5 	-0.5			}		
		dev postlandfile "~/rg_vsfacto.ps"
		graphanax 			3493 	3542	 0.0		
		dev x11		
		graphanax 			3493 	3542	-0.5				
	  #BIG rg data
		do z=3593, 3642, 1{
		clear
		graphden 			$z 	-2.5   0.0			
		clear
		graphden 			$z	-2.5 	-0.5			}			
		dev postlandfile "~/brg_vsfacto.ps"
		graphanax 			3593 	3642	 0.0		
		dev x11		
		graphanax 			3593 	3642	-0.5				  		
	  #extra
		dev postlandfile "~/m_vs_a.ps"
		extra
		extra2
		dev x11
	  #density	
		dev postlandfile "~/density.ps"
		den1 2093
		den2 2101
		den3 2117
		den4 2139
		x11					  
	  #big rg
		dev postlandfile "~/rg035070den.ps"
		den1 3517
		den2 3617
		plw2	
		dev x11		
	  #velocity
		dev postlandfile "~/velocity.ps"
		vel
		dev x11
	  #powerlaw
		dev postlandfile "~/powerlaw.ps"
		plw
		dev x11

plw2		#----num  graphs average file for paper
		expand 0.7
		set rrmin=-2.5
		set rrmax=0.0
		
		define num 3517
		data "~/find_density/$!num/avg_$!num.dat"
		read {i 1 ii 2 iii 3 iv 4 v 5 vi 6 vii 7 ix 9 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 }
		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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 u[$(k[0])]=$(ii[$j])
                  set s[$(k[0])]=$(iii[$j])
		      set w[$(k[0])]=$(iv[$j])
                  set x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$j])
		      set vrms2[$(k[0])]=$(vrms[$j])
			set k=1+$(k[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)
		location 3500 15500 3500 16000
		xlabel log [r/r_{200}] 
            #ylabel log[\rho/\sigma^3 (M_{\odot}km^3/Mpc^3/s^3) ]
		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+.5+1.875*s)
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15		
		label \alpha_{fit}  = $(tempfmt[0])
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		set alpha=$((-1)*($a))
		set alpharms=$rms
		move 90 90
		label [C] 		
				
		define num 3617
		data "~/find_density/$!num/avg_$!num.dat"
		read {i 1 ii 2 iii 3 iv 4 v 5 vi 6 vii 7 ix 9 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}
		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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 u[$(k[0])]=$(ii[$j])
                  set s[$(k[0])]=$(iii[$j])
		      set w[$(k[0])]=$(iv[$j])
                  set x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$j])
		      set vrms2[$(k[0])]=$(vrms[$j])
			set k=1+$(k[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)
		location 18500 31000 3500 16000
		xlabel log [r/r_{200}] 
            #ylabel log[\rho/\sigma^3 (M_{\odot}km^3/Mpc^3/s^3) ]
		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+.5+1.875*s)
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \alpha_{fit}  = $(tempfmt[0])
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		set alpha=$((-1)*($a))
		set alpharms=$rms
		move 90 90
		label [D] 		
		
		DEL Mhalo alpha alpharms aout bout 
		DEL chiwq cout facto i ii 
		DEL iii iv ix chisq k 
		DEL q r200 rr200 rrmax rrmin 
		DEL s temp tempfmt u v 
		DEL vi vii vrms vrms2 w 
		DEL xaxis xi y yaxis yfit 
		DEL z x x x x
		echo --------
		list set
		echo --------			
				  		
den1 	1	#----num  graphs average file for paper
		define num $1
		expand 0.7
		set rrmin=-2.5
		set rrmax=0.0
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 iv 4 vi 6}
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {conx1 5 conx2 6 conx3 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 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 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 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])
	     #DENSITY/rho_o
		#divide by rho_0 in M_sol/cm^3
		set iv=iv+63-LG(2.3563)
		set w=w+63-LG(2.3563)
		location 3500 15500 18500 31000 
		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
		ltype 3
			set preverror=99999999
			set chi=.0338*LN(r200)-1.7922
			#echo $(1/$(chi[0]))
			#chi=1/concentration
			#   =r_s/r_200
			do const=3.8,4.3,.01{
			  #dentrial is trial values of log(den/den_o)
			  set dentrial=$const-iii-2*LG(1+((10**iii)/$(chi[0])))
			  set error=0
			  do j=0,dimen(iii)-1,1{
			    set error=$(error[0])+abs($(dentrial[$j])-$(iv[$j]))
			  }
			  if ($(error[0])<$(preverror[0])){
			    set preverror=$(error[0])
			    set holdconst=$const
			  }
			}
			if($(holdconst[0])==3.5){echo decrease limits>>>>>>>>>>>>>>>}
			if($(holdconst[0])==4.3){echo increase limits>>>>>>>>>>>>>>>}
			set dentrial=$(holdconst[0])-iii-2*LG(1+((10**iii)/$(chi[0])))
			con iii (dentrial+2*iii)
		ltype 0
		ptype 1 1
            points iii (iv+2*iii)
	 	con iii (iv+2*iii)
		lsq s w s q rms 
		con s (q+2*s+.1)
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \beta_{fit}  = $(tempfmt[0])
		set beta=$((-1)*($a))
		set betarms=$rms
		move 90 90
		label [A] 
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		
		DEL 1concen 2concen 3concen Mhalo beta 
		DEL betarms chi concen conx1 conx2
		DEL conx3 dentrial error facto holdconst 
		DEL iii iv k preverror q 
		DEL r200 x x s temp
		DEL tempfmt vi w y x
		echo ---
		list set
		echo ---
		
den2	1     define num $1
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 iv 4 vi 6}
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {conx1 5 conx2 6 conx3 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 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 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 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])
	     #DENSITY/rho_o
		#divide by rho_0 in M_sol/cm^3
		set iv=iv+63-LG(2.3563)
		set w=w+63-LG(2.3563)
		location 18500 31000 18500 31000 
		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
		ltype 3
			set preverror=99999999
			set chi=.0338*LN(r200)-1.7922
			#echo $(1/$(chi[0]))
			#chi=1/concentration
			#   =r_s/r_200
			do const=3.8,4.3,.01{
			  #dentrial is trial values of log(den/den_o)
			  set dentrial=$const-iii-2*LG(1+((10**iii)/$(chi[0])))
			  set error=0
			  do j=0,dimen(iii)-1,1{
			    set error=$(error[0])+abs($(dentrial[$j])-$(iv[$j]))
			  }
			  if ($(error[0])<$(preverror[0])){
			    set preverror=$(error[0])
			    set holdconst=$const
			  }
			}
			if($(holdconst[0])==3.5){echo decrease limits>>>>>>>>>>>>>>>}
			if($(holdconst[0])==4.3){echo increase limits>>>>>>>>>>>>>>>}
			set dentrial=$(holdconst[0])-iii-2*LG(1+((10**iii)/$(chi[0])))
			con iii (dentrial+2*iii)
		ltype 0
		ptype 1 1
            points iii (iv+2*iii)
	 	con iii (iv+2*iii)
		#ctype red
		ptype 4 1
		lsq s w s q rms 
		con s (q+2*s+.1)
		#ctype default
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \beta_{fit}  = $(tempfmt[0])
		set beta=$((-1)*($a))
		set betarms=$rms
		move 90 90
		label [B] 
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		
		DEL 1concen 2concen 3concen Mhalo beta 
		DEL betarms chi concen conx1 conx2
		DEL conx3 dentrial error facto holdconst 
		DEL iii iv k preverror q 
		DEL r200 x x s temp
		DEL tempfmt vi w y x
		echo ---
		list set
		echo ---
				
den3	1	define num $1
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 iv 4 vi 6}
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {conx1 5 conx2 6 conx3 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 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 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 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])
	     #DENSITY/rho_o
		#divide by rho_0 in M_sol/cm^3
		set iv=iv+63-LG(2.3563)
		set w=w+63-LG(2.3563)
		location 3500 15500 3500 16000
		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
		ltype 3
			set preverror=99999999
			set chi=.0338*LN(r200)-1.7922
			#echo $(1/$(chi[0]))
			#chi=1/concentration
			#   =r_s/r_200
			do const=3.8,4.3,.01{
			  #dentrial is trial values of log(den/den_o)
			  set dentrial=$const-iii-2*LG(1+((10**iii)/$(chi[0])))
			  set error=0
			  do j=0,dimen(iii)-1,1{
			    set error=$(error[0])+abs($(dentrial[$j])-$(iv[$j]))
			  }
			  if ($(error[0])<$(preverror[0])){
			    set preverror=$(error[0])
			    set holdconst=$const
			  }
			}
			if($(holdconst[0])==3.5){echo decrease limits>>>>>>>>>>>>>>>}
			if($(holdconst[0])==4.3){echo increase limits>>>>>>>>>>>>>>>}
			set dentrial=$(holdconst[0])-iii-2*LG(1+((10**iii)/$(chi[0])))
			con iii (dentrial+2*iii)
		ptype 1 1
		ltype 0
            points iii (iv+2*iii)
	 	con iii (iv+2*iii)
		ptype 4 1
		lsq s w s q rms 
		con s (q+2*s+.1)
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \beta_{fit}  = $(tempfmt[0])
		set beta=$((-1)*($a))
		set betarms=$rms
		move 90 90
		label [C]  
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		
		DEL 1concen 2concen 3concen Mhalo beta 
		DEL betarms chi concen conx1 conx2
		DEL conx3 dentrial error facto holdconst 
		DEL iii iv k preverror q 
		DEL r200 x x s temp
		DEL tempfmt vi w y x
		echo ---
		list set
		echo ---
		
den4 1	define num $1  
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 iv 4 vi 6}
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {conx1 5 conx2 6 conx3 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 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 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 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])
	     #DENSITY/rho_o
		#divide by rho_0 in M_sol/cm^3
		set iv=iv+63-LG(2.3563)
		set w=w+63-LG(2.3563)
		location 18500 31000 3500 16000
		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
		ltype 3
			set preverror=99999999
			set chi=.0338*LN(r200)-1.7922
			#echo $(1/$(chi[0]))
			#chi=1/concentration
			#   =r_s/r_200
			do const=3.8,4.3,.01{
			  #dentrial is trial values of log(den/den_o)
			  set dentrial=$const-iii-2*LG(1+((10**iii)/$(chi[0])))
			  set error=0
			  do j=0,dimen(iii)-1,1{
			    set error=$(error[0])+abs($(dentrial[$j])-$(iv[$j]))
			  }
			  if ($(error[0])<$(preverror[0])){
			    set preverror=$(error[0])
			    set holdconst=$const
			  }
			}
			if($(holdconst[0])==3.5){echo decrease limits>>>>>>>>>>>>>>>}
			if($(holdconst[0])==4.3){echo increase limits>>>>>>>>>>>>>>>}
			set dentrial=$(holdconst[0])-iii-2*LG(1+((10**iii)/$(chi[0])))
			con iii (dentrial+2*iii)
		ptype 1 1
		ltype 0
            points iii (iv+2*iii)
	 	con iii (iv+2*iii)
		ptype 4 1
		lsq s w s q rms 
		con s (q+2*s+.1)
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \beta_{fit}  = $(tempfmt[0])
		set beta=$((-1)*($a))
		set betarms=$rms
		move 90 90
		label [D] 	
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		
		DEL 1concen 2concen 3concen Mhalo beta 
		DEL betarms chi concen conx1 conx2
		DEL conx3 dentrial error facto holdconst 
		DEL iii iv k preverror q 
		DEL r200 rrmax rrmin s temp
		DEL tempfmt vi w y x
		echo ---
		list set
		echo ---
		
vel 		#----num  graphs average file for paper
		expand 0.7
		set rrmin=-2.5
		set rrmax=0.0
		
		define num 2093
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 v 5 vi 6 vii 7 ix 9 xi 11 }
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$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])
	     #V RMS
	      location 3500 15500 24750 31000 
		set vrms=vrms-lg(1E5)
		set vrms2=vrms2-lg(1E5)
		#change units to M_sol*km^3/(mpc^3*s^3)
		#xlabel log(r/r_{200})  
		ylabel log[V_{rms} (km/s)]
	 	limits -4.5 2.0 .1 2.5
            box 0 2 0 0
		ptype 1 1
            points iii vrms	
		con iii vrms
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		#ctype red
		#ptype 4 1
		lsq s vrms2 s q rms 
		#con s q
		#ctype default
		move 90 86
		label [A] 
		#beta graph
		location 3500 15500 18500 24750 
		set dimen(betax)=2
		set dimen(betay)=2
		xlabel log [r/r_{200}]  
            ylabel \beta_{vel.disp.}
		limits -4.5 2.0 -.5 1.1
		box 1 2 0 0
		ptype 1 1
		points iii ix
		#con iii ix
		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
		#label y=$a*x+$b	
		#$a
		#$b
		#$rms
		#label rms: $rms11		
		move 90 86
 	
		
		define num 2101
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 v 5 vi 6 vii 7 ix 9 xi 11 }
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 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(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$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])
	     #V RMS
	      location 18500 31000 24750 31000
		set vrms=vrms-lg(1E5)
		set vrms2=vrms2-lg(1E5)
		#change units to M_sol*km^3/(mpc^3*s^3)
		#xlabel log(r/r_{200})  
		ylabel log[V_{rms} (km/s)]
	 	limits -4.5 2.0 .1 2.5
            box 0 2 0 0
		ptype 1 1
            points iii vrms	
		con iii vrms
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		#ctype red
		#ptype 4 1
		lsq s vrms2 s q rms 
		#con s q
		#ctype default
		move 90 86
		label [B] 
		 #beta graph
		location 18500 31000 18500 24750 
		set dimen(betax)=2
		set dimen(betay)=2
		xlabel log [r/r_{200}]  
            ylabel \beta_{vel.disp.}
		limits -4.5 2.0 -.5 1.1
		box 1 2 0 0
		ptype 1 1
		points iii ix
		#con iii ix
		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
		#label y=$a*x+$b	
		#$a
		#$b
		#$rms
		#relocate -4 -5.5
		#label rms: $rms11		
		move 90 86
 
		
		define num 2117
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 v 5 vi 6 vii 7 ix 9 xi 11 }
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 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(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$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])
	     #V RMS
	      location 3500 15500 9750 16000
		set vrms=vrms-lg(1E5)
		set vrms2=vrms2-lg(1E5)
		#change units to M_sol*km^3/(mpc^3*s^3)
		#xlabel log(r/r_{200})  
		ylabel log[V_{rms} (km/s)]
	 	limits -4.5 2.0 .1 2.5
            box 0 2 0 0
		ptype 1 1
            points iii vrms	
		con iii vrms
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		#ctype red
		#ptype 4 1
		lsq s vrms2 s q rms 
		#con s q
		#ctype default
		move 90 86
		label [C] 
		#beta graph
		location 3500 15500 3500 9750
		set dimen(betax)=2
		set dimen(betay)=2
		xlabel log [r/r_{200}]  
            ylabel \beta_{vel.disp.}
		limits -4.5 2.0 -.5 1.1
		box 1 2 0 0
		ptype 1 1
		points iii ix
		#con iii ix
		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
		#label y=$a*x+$b	
		#$a
		#$b
		#$rms
		#relocate -4 -5.5
		#label rms: $rms11		
		move 90 86 
				
		define num 2139
		data "~/find_density/$!num/avg_$!num.dat"
		read {iii 3 v 5 vi 6 vii 7 ix 9 xi 11 }
		data "~/find_density/$!num/summary_$!num.dat"
		read {r200 5 Mhalo 6 facto 7}
		data "~/find_density/$!num/navarrofit_$!num.dat"
		read {aout 1 bout 2 cout 3 chisq 4 conx1 5 conx2 6 conx3 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(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$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])
	     #V RMS
	      location 18500 31000 9750 16000
		set vrms=vrms-lg(1E5)
		set vrms2=vrms2-lg(1E5)
		#change units to M_sol*km^3/(mpc^3*s^3)
		#xlabel log(r/r_{200})  
		ylabel log[V_{rms} (km/s)]
	 	limits -4.5 2.0 .1 2.5
            box 0 2 0 0
		ptype 1 1
            points iii vrms	
		con iii vrms
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		#ctype red
		#ptype 4 1
		lsq s vrms2 s q rms 
		#con s q
		#ctype default
		move 90 86
		label [D] 
		
		#beta graph
		location 18500 31000 3500 9750
		set dimen(betax)=2
		set dimen(betay)=2
		xlabel log [r/r_{200}]  
            ylabel \beta_{vel.disp.}
		limits -4.5 2.0 -.5 1.1
		box 1 2 0 0
		ptype 1 1
		points iii ix
		#con iii ix
		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
		#label y=$a*x+$b	
		#$a
		#$b
		#$rms
		#relocate -4 -5.5
		#label rms: $rms11		
		move 90 86
				
		DEL 1concen 2concen 3concen Mhalo aout 
		DEL betax betay bout chisq concen 
		DEL conx1 conx2 conx3 cout facto 
		DEL iii ix k q r200
		DEL rrmax rrmin s u v 
		DEL vi vii vrms vrms2 w 
		DEL x xi y z x
		
		echo -------
		list set
		echo -------
		
plw  	#----num  graphs average file for paper
		expand 0.7
		set rrmin=-2.5
		set rrmax=0.0
		
		define num 2093
		data "~/find_density/$!num/avg_$!num.dat"
		read {i 1 ii 2 iii 3 iv 4 v 5 vi 6 vii 7  ix 9 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 }
		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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 u[$(k[0])]=$(ii[$j])
                  set s[$(k[0])]=$(iii[$j])
		      set w[$(k[0])]=$(iv[$j])
                  set x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$j])
		      set vrms2[$(k[0])]=$(vrms[$j])
			set k=1+$(k[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)
		location 3500 15500 18500 31000 
		xlabel log [r/r_{200}] 
            #ylabel log[\rho/\sigma^3 (M_{\odot}km^3/Mpc^3/s^3) ]
		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+.5+1.875*s)
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \alpha_{fit}  = $(tempfmt[0])
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		set alpha=$((-1)*($a))
		set alpharms=$rms
		move 90 90
		label [A]  
		
		define num 2101
		data "~/find_density/$!num/avg_$!num.dat"
		read {i 1 ii 2 iii 3 iv 4 v 5 vi 6 vii 7 ix 9 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 }
		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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 u[$(k[0])]=$(ii[$j])
                  set s[$(k[0])]=$(iii[$j])
		      set w[$(k[0])]=$(iv[$j])
                  set x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$j])
		      set vrms2[$(k[0])]=$(vrms[$j])
			set k=1+$(k[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)
		location 18500 31000 18500 31000
		xlabel log [r/r_{200}] 
            #ylabel log[\rho/\sigma^3 (M_{\odot}km^3/Mpc^3/s^3) ]
		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+.5+1.875*s)
		set tempfmt=sprintf('%5.3f',$((-1)*($a)))
		move 5 15
		label \alpha_{fit}  = $(tempfmt[0])
		move 5 7
		label rms = $(sprintf('%5.3f',$rms))
		set alpha=$((-1)*($a))
		set alpharms=$rms
		move 90 90
		label [B]      

		define num 2117
		data "~/find_density/$!num/avg_$!num.dat"
		read {i 1 ii 2 iii 3 iv 4 v 5 vi 6 vii 7 ix 9 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 }
		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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 u[$(k[0])]=$(ii[$j])
                  set s[$(k[0])]=$(iii[$j])
		      set w[$(k[0])]=$(iv[$j])
                  set x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$j])
		      set vrms2[$(k[0])]=$(vrms[$j])
			set k=1+$(k[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)
		location 3500 15500 3500 16000
		xlabel log [r/r_{200}] 
            #ylabel log[\rho/\sigma^3 (M_{\odot}km^3/Mpc^3/s^3) ]
		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+.5+1.875*s)
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \alpha_{fit}  = $(tempfmt[0])
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		set alpha=$((-1)*($a))
		set alpharms=$rms
		move 90 90
		label [C] 		
				
		define num 2139
		data "~/find_density/$!num/avg_$!num.dat"
		read {i 1 ii 2 iii 3 iv 4 v 5 vi 6 vii 7 ix 9 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}
		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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(k[0])
		set dimen(z)=$(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 u[$(k[0])]=$(ii[$j])
                  set s[$(k[0])]=$(iii[$j])
		      set w[$(k[0])]=$(iv[$j])
                  set x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$j])
                  set z[$(k[0])]=$(vii[$j])
		      set vrms2[$(k[0])]=$(vrms[$j])
			set k=1+$(k[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)
		location 18500 31000 3500 16000
		xlabel log [r/r_{200}] 
            #ylabel log[\rho/\sigma^3 (M_{\odot}km^3/Mpc^3/s^3) ]
		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+.5+1.875*s)
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \alpha_{fit}  = $(tempfmt[0])
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		set alpha=$((-1)*($a))
		set alpharms=$rms
		move 90 90
		label [D] 		
		
		DEL Mhalo alpha alpharms aout bout 
		DEL chisq cout facto i ii 
		DEL iii iv ix k q 
		DEL r200 rr200 rrmax rrmin s 
		DEL temp tempfmt u v vi 
		DEL vii vrms vrms2 w x 
		DEL xaxis xi y yaxis yfit 
		DEL z 1concen 2concen 3concen betax 
		DEL betay concen conx1 conx2 conx3
		echo ------
		list set
		echo ------
				
				
graphdenx 3 #----num  graphs average file for paper
		expand 0.7
		define num $1
		define rrmin $2
		define rrmax $3
		
		data "~/find_density/$!num/avg_$!num.dat"
		read {i 1 ii 2 iii 3 iv 4 v 5 vi 6 vii 7 viii 8 ix 9 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}
		
		 
		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(u)=$(k[0])
		set dimen(s)=$(k[0])
		set dimen(w)=$(k[0])
		set dimen(x)=$(k[0])
		set dimen(y)=$(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 u[$(k[0])]=$(ii[$j])
                  set s[$(k[0])]=$(iii[$j])
		      set w[$(k[0])]=$(iv[$j])
                  set x[$(k[0])]=$(v[$j])
		      set y[$(k[0])]=$(vi[$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])
		
		#DENSITY/rho_o
		#divide by rho_0 in M_sol/cm^3
		set iv=iv+63-LG(2.3563)
		set xi=xi+63-LG(2.3563)
		set w=w+63-LG(2.3563)
		location 3500 15500 18500 31000 
		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
		ltype 3
	#		set preverror=99999999
	#		set chi=.0338*LN(r200)-1.7922
	#		#echo $(1/$(chi[0]))
	#		#chi=1/concentration
	#		#   =r_s/r_200
	#		do const=2,5,.1{
	#		#dentrial is trial values of log(den/den_o)
	#		set dentrial=$const-LG((10**iii/10**$(chi[0]))*(1+10**iii/10**$(chi[0]))**2) 
	#		set error=0
	#		do j=0,dimesn(iii)-1,1{
	#		set error=$(error[0])+abs($(dentrial[$j])-$(iv[$j]))}
	#		if ($(error[0])<$(preverror[0])){
	#		set preverror=$(error[0])
	#		set holdconst=$const}}
	#		if($(holdconst[0])==2){echo increase limits>>>>>>>>>>>>>>>}
	#		if($(holdconst[0])==5){echo increase limits>>>>>>>>>>>>>>>}
	#		set dentrial=$(holdconst[0])-LG((10**iii/10**$(chi[0]))*(1+10**iii/10**$(chi[0]))**2) 	
	#		con iii (dentrial+2*iii)
		ltype 0
		ptype 1 1
            #ctype magenta
		#con iii (xi+2*iii)
		#ctype black
		ptype 1 1
            points iii (iv+2*iii)
	 	con iii (iv+2*iii)
		#ctype red
		ptype 4 1
		lsq s w s q rms 
		con s (q+2*s+.1)
		#ctype default
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15
		label \beta_{fit}  = $(tempfmt[0])
		move 90 90
		label [A] 
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		set beta=$((-1)*($a))
		set betarms=$rms
			      
		#V RMS
	      location 18500 31000 18500 31000 
		set vrms=vrms-lg(1E5)
		set vrms2=vrms2-lg(1E5)
		#change units to M_sol*km^3/(mpc^3*s^3)
		xlabel log(r/r_{200})  
		ylabel log[V_{rms} (km/s)]
	 	limits -4.5 2.0 .1 3.0
		box
		ptype 1 1
            points iii vrms	
		con iii vrms
		ltype 1
		con 1concen concen
		ltype 2
		con 2concen concen
		ltype 3
		con 3concen concen
		ltype 0
		#ctype red
		#ptype 4 1
		lsq s vrms2 s q rms 
		#con s q
		#ctype default
		move 90 86
		label [B] 
		set vela=$a
		set velrms=$rms
						
		#beta graph
		location 18500 31000 3500 16000 
		set dimen(betax)=2
		set dimen(betay)=2
		xlabel log [r/r_{200}]  
            ylabel \beta_{vel.disp.}
		limits -4.5 2.0 -.5 1.1
		box 1 2 0 0
		ptype 1 1
		points iii ix
		#con iii ix
		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
		#label y=$a*x+$b	
		#$a
		#$b
		#$rms
		#relocate -4 -5.5
		#label rms: $rms11		
		move 90 86
		label [D]
						
		#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)
		location 3500 15500 3500 16000
		xlabel log [r/r_{200}] 
            #ylabel log[\rho/\sigma^3 (M_{\odot}km^3/Mpc^3/s^3) ]
		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+.5+1.875*s)
		set temp=$((-1)*($a))
		set tempfmt=sprintf('%5.3f',temp)
		move 5 15		
		label \alpha_{fit}  = $(tempfmt[0])
		set temp=$rms
		set tempfmt=sprintf('%5.3f',temp)
		move 5 7
		label rms = $(tempfmt[0])
		move 90 90
		label [C] 
		set alpha=$((-1)*($a))
		set alpharms=$rms
		
				
		##navarro, fig4
		#location 18500 31000 3500 16000
		#xlabel log [\sigma^3/\rho (Mpc^3s^3/M_{\odot}/km^3) ]
            #ylabel log [dM/dlog(\sigma^3/\rho)]
		#limits 60 100 -5 20
		#box
		#ptype 1 1
            #points xaxis yaxis
		#con xaxis yaxis
		#ptype 1 1
		#set yfit=xaxis*xaxis*aout+xaxis*bout+cout
		#con xaxis yfit
            #set temp=(-1)*$(aout[0])
		#set tempfmt=sprintf('%5.3f',temp)
		#relocate $($fx1+($fx2-$fx1)*.01*5) $($fy1+($fy2-$fy1)*.01*15)
		#label \gamma_{fit}=$(tempfmt[0])
		#set temp=$(chisq[0])
		#set tempfmt=sprintf('%8.3f',temp)
		#relocate $($fx1+($fx2-$fx1)*.01*5) $($fy1+($fy2-$fy1)*.01*7)
		#label chisq=$(tempfmt[0])
		#set gamm=$(aout[0])
		#set gammarms=$(chisq[0])
		#relocate $($fx1+($fx2-$fx1)*.01*90) $($fy1+($fy2-$fy1)*.01*90)
		#label [E] 	
		
		set facto=sprintf ('%5.2f',facto)
		#toplabel .                      FILE: $num  FACTO: $(facto[0])	
		DEL i ii iii iv v
		DEL r200 Mhalo facto denstiy dentrial 
		DEL aout bout cout chisq vrms 
		DEL vi vii viii ix holdconst
		DEL rr200 xaxis yaxis yfit holdchi
		DEL u s w x y 
		DEL z vrms2 k xi concen
		DEL conx1 conx2 conx3 1concen 2concen
		DEL 3concen q rrmax rrmin viii2
		DEL alpha betax betay temp tempfmt
		DEL alpharms beta betarms gamm gammarms
		DEL preverror radius rvect dentrial density
		DEL chi const error 2concen vela
		delete velrms
		echo ----
		list set
		echo ----
		


graphanax 3	#----num num2: graphs analysis for paper
		define num $1
		define num2 $2
		define rrmax $3
		expand 0.7
		erase
		#if($num==3493){
		#toplabel .                      RG DATA}
		#if($num==3293){
		#toplabel .                      FIG5 DATA}
		#if($num==2093){
		#toplabel .                      STANDARD DATA}
		
		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 alpha alpharms beta betarms conx1
		  	DEL conx2 conx3 facto gamm gammarms
			DEL aout bout chisq cout vela
			delete velrms
		  	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   BETA
		location 3500 15500 18500 31000
		xlabel Scaling Parameter
            ylabel \beta_{fit}
		limits 0.30 2.7 0.7 2.0
            box 1 2 0 0
		ptype 5 1
            points facto beta
		ptype 1 1
		lsq facto beta facto q rms 
		move 93 93
		label (a
		
		#DEN/VRMS^3   ALPHA
		location 18500 31000 18500 31000
		xlabel Scaling Parameter
		ylabel \alpha_{fit}
		limits 0.30 2.7 1.5 2.4
		box 1 2 0 0
		set dimen(nfw)=2
		set dimen(nfwy)=2
		set nfw[0]=1.875
		set nfw[1]=1.875
		set nfwy[0]=.35
		set nfwy[1]=2.0
		ptype 1 1
		points nfwy nfw
		ltype 3
		con nfwy nfw
		relocate $(nfwy[1]) $(nfw[0])
		label NFW
		ltype 0
		ptype 5 1
            points facto alpha
		ptype 1 1
		lsq facto alpha facto q rms 
		move 93 93
		label (b
		
		#concentration factors
		location 3500 15500 3500 16000
		xlabel Scaling Parameter
		ylabel r_{\beta=1,2,3}
		limits 0.30 2.7 -4.5 2
		box
		ptype 4 3
		points facto conx1
		ptype 3 3
		points facto conx3
		ptype 5 2
		points facto conx2
		ptype 1 1
		move 5 90
		label Squares: r_{\beta=1}  Stars: r_{\beta=2}  Triangles: r_{\beta=3}
		move 93 93
		label (d
						
		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))
		  	set sum1=0
		  	set sum2=0
		  	set sum3=0
		  	set sum4=0
		  	set sum5=0
		  	do j=0,$($(k[0])-1),1{
		  		set sum1=$(sum1[0])+$(masslost[$j])
		  		set sum2=$(sum2[0])+$(Mhalo[$j])
		  		set sum3=$(sum3[0])+$(shelllost[$j])
		  		set sum4=$(sum4[0])+$(r200[$j])
		  		set sum5=$(sum5[0])+$(mfraclost[$j])}
		  	set sum1=$(sum1[0])/$(k[0])
		  	set sum2=$(sum2[0])/$(k[0])
		  	set sum3=$(sum3[0])/$(k[0])
		  	set sum4=$(sum4[0])/$(k[0])
		  	set sum5=$(sum5[0])/$(k[0])
		  	if ($rrmax==0.0){
			write + "~/find_density/.temp_mass_$!num.dat" $(numm[0]) $(sum3[0]) $(sum4[0]) $(sum2[0]) $(facto[0]) $(sum1[0]) $(sum5[0]) 
			}
		  	if ($rrmax==-0.5){
			write + "~/find_density/.temp_mass_$!num.datr" $(numm[0]) $(sum3[0]) $(sum4[0]) $(sum2[0]) $(facto[0]) $(sum1[0]) $(sum5[0])
		  	}		
			DEL sum1 sum2 sum3 sum4 sum5
		}
		if ($rrmax==0.0){
		data "~/find_density/.temp_mass_$!num.dat"
		read {numm 1 shelllost 2 r200 3 Mhalo 4 facto 5 masslost 6 mfraclost 7}
		}
		if ($rrmax==0.0){
		data "~/find_density/.temp_mass_$!num.datr"
		read {numm 1 shelllost 2 r200 3 Mhalo 4 facto 5 masslost 6 mfraclost 7}
		}
				
		#Mhalo and masslost
		location 18500 31000 3500 16000
		set Mhalo=Mhalo/1e12
		set masslost=masslost/1e12
		xlabel Scaling Parameter
            ylabel Mass (10^{12} M_{\odot})
		set testlow=$(Mhalo[0])
		set testhigh=$(Mhalo[0])
		do z=0,dimen(Mhalo)-1,1{
			if ($(Mhalo[$z])>$(testhigh[0])){set testhigh=$(Mhalo[$z])}
		  	if ($(Mhalo[$z])<$(testlow[0])){set testlow=$(Mhalo[$z])}
		  	if ($(masslost[$z])>$(testhigh[0])){set testhigh=$(masslost[$z])}
		  	if ($(masslost[$z])<$(testlow[0])){set testlow=$(masslost[$z])}}
		limits 0.30 2.7 -.5 $($(testhigh[0])+$(testhigh[0])*.1)
		
		#limits 0.30 2.7 -2.5 32
		box
		ptype 4 0
		points facto Mhalo
		ptype 4 1 
		points facto masslost
		move 40 90
		label Squares: Final Halo Mass 
		move 40 83
		label Crosses: Mass Lost
		move 93 93
		label (c
				
		DEL Mhalo Mhalo facto facto halo 
		DEL halo k masslost masslost mfraclost
		DEL mfraclost numm numm r200 r200 
		DEL shelllost shelllost testhigh testlow x
		DEL a11 a21 a22 a31 aout 
  		DEL b11 b21 b22 b31 bout 
		DEL chisq conx1 conx2 conx3 cout
		DEL facto intcpt1 intcpt2 intcpt3 intcpt4 
		DEL intcpt5 intcpt6 nfw nfwy rms1 
		DEL rms11 rms2 rms21 rms22 rms3 
		DEL rms31 rms4 rms5 rms6 slope1 
		DEL slope2 slope3 slope4 slope5 slope6
		DEL a12 a32 b12 b32 q
		DEL alpha alpharms beta betarms gamm 
		DEL gammarms rms12 rms32 yfit2 sum5
		DEL vela velrms x x x
		echo ----
		list set
		echo ----
						
		
		
extra		#----graphs mass vs alpha (power law slope)
		#standard data
		expand 1.0001
		
	     #rrmax=0.0
		#standard data
		data "~/find_density/.temp_ana_2093.dat"
		read {aalpha 4}
		data "~/find_density/.temp_mass_2093.dat"
		read {aMhalo 4}
		#rg data
		data "~/find_density/.temp_ana_3493.dat"
		read {calpha 4}
		data "~/find_density/.temp_mass_3493.dat"
		read {cMhalo 4}
	      #big rg data
		data "~/find_density/.temp_ana_3593.dat"
		read {dalpha 4}
		data "~/find_density/.temp_mass_3593.dat"
		read {dMhalo 4}
	    		
		#mass vs. alpha
		set aMhalo=lg(aMhalo)
		set cMhalo=lg(cMhalo)
		set dMhalo=lg(dMhalo)
		location 3500 31000 3500 31000
		xlabel log (Mass M_{\odot})
		ylabel \alpha_{fit}
		limits 10 14 1.5 2.55
		box
		ltype 0
		
		ptype 3 3
		points aMhalo aalpha
		#con aMhalo aalpha
		
		ptype 4 3
		points cMhalo calpha
		#con cMhalo calpha
		
		ptype 5 3
		points dMhalo dalpha
		#con dMhalo dalpha
		
		ctype black
		expand .9
		#relocate $(aMhalo[0]) $($(aalpha[0])-0.05)
		#label Scaling Parameter=0.40
		#relocate $($(cMhalo[0])-1) $($(calpha[0])-0.05)
		#label Scaling Parameter=0.40
		move 60 96
		label Triangle:   Galactic Halos 
		move 60 92
		label Box:        Group Halos
		move 60 88
		label Pentagon:  Cluster Halos 
		move 5 31
		label NFW \alpha_{fit}=1.875
		expand 1.001
		set dimen(nfw)=2
		set dimen(nfwy)=2
		set nfw[0]=1.875
		set nfw[1]=1.875
		set nfwy[0]=-.5
		set nfwy[1]=15
		ltype 3
		con nfwy nfw
		
		DEL aMhalo bMhalo cMhalo aalpha balpha 
		DEL calpha dalpha nfw nfwy dMhalo
		delete dMhalo
		
		echo ------------
		list set
		echo ------------
		
		
extra2	#----graphs mass vs alpha (power law slope)
		#standard data
		expand 1.0001
		
	     #rrmax=-0.5
		#standard data
		data "~/find_density/.temp_ana_2093.datr"
		read {aalphar 4}
		data "~/find_density/.temp_mass_2093.datr"
		read {aMhalor 4}
		#rg data
		data "~/find_density/.temp_ana_3493.datr"
		read {calphar 4}
		data "~/find_density/.temp_mass_3493.datr"
		read {cMhalor 4}
	      #big rg data
		data "~/find_density/.temp_ana_3593.datr"
		read {dalphar 4}
		data "~/find_density/.temp_mass_3593.datr"
		read {dMhalor 4}	
		
		#mass vs. alpha
		set aMhalor=lg(aMhalor)
		set cMhalor=lg(cMhalor)
		set dMhalor=lg(dMhalor)
		location 3500 31000 3500 31000
		xlabel log (Mass M_{\odot})
		ylabel \alpha_{fit}
		limits 10 14 1.5 2.55
		box
		ltype 0
		
		ptype 3 0
		points aMhalor aalphar
		#con aMhalor aalphar
		
		ptype 4 0
		points cMhalor calphar
		#con cMhalor calphar
		
		ptype 5 0
		points dMhalor dalphar
		
		ctype black
		DEL aMhalor bMhalor cMhalor aalphar balphar 
		DEL calphar dalphar dMhalor dMhalor x
		
		echo ------------
		list set
		echo ------------