;$Id: f_comp.pro,v 1.1 2016/11/14 23:52:28 xiangyu Exp $
size=1.5
size2=1.5
if !d.name eq 'PS' then begin
  device,xsize=30,ysize=12,yoffset=3
  !p.charthick=3 & !p.thick=5 & !x.thick=3 & !y.thick=3
  !p.charsize=size
endif else begin
   !p.charsize=size
end
;
!p.multi=[0,2,1]
!x.margin=[7.5,2.5]
!y.margin=[3.5,.5]
;
;---------size spectra of the Eulerian and swarm model-------------------------

!x.title='!8r!6 [!7l!6m]'
!y.title='!8f!6 !8r!6!dini!n/!8n!6!d0!n'

xr=[4,2500]
yr=[1.e-10,3.e0]
microm=1e-6
thick1=6
dir1='$PENCIL_HOME/xiangyu/turb_coag/grav/0D/nocondens3457coag_a1em5_sig02_s001a_fast_a15em6_1p00543_MC3b/'
dir2='$PENCIL_HOME/xiangyu/turb_coag/grav/0D/nocondens3457coag_a1em5_sig02_s001a_fast_a15em6_1p00543/'
dir3='../SW32_coag_grav1_ap10um_shima_262kp_32procs'
dir4='../SW32_coag_grav1_ap10um_shima_1048kp_32procs'

;restore,file=dir1+'/data/pdf.sav' & radii1=radii & prob1=prob & time1=t
restore,file=dir1+'/data/ad.sav'
restore,file=dir1+'/data/nd_all.sav'
radii1_0=ad& prob1_0=nd_all[0,*]
radii1_1=ad& prob1_1=nd_all[1,*]
radii1_2=ad& prob1_2=nd_all[2,*]
radii1_3=ad& prob1_3=nd_all[3,*]
;radii2_4=md_all[4,*]& prob2_4=nd_all[4,*]

restore,file=dir2+'/data/histo/histo0.sav' & radii2_0=adm & prob2_0=nd & time2_0=t
restore,file=dir2+'/data/histo/histo10.sav'& radii2_1=adm & prob2_1=nd & time2_1=t
restore,file=dir2+'/data/histo/histo20.sav'& radii2_2=adm & prob2_2=nd & time2_2=t
restore,file=dir2+'/data/histo/histo30.sav'& radii2_3=adm & prob2_3=nd & time2_3=t
restore,file=dir2+'/data/histo/histo40.sav'& radii2_4=adm & prob2_4=nd & time2_4=t
restore,file=dir3+'/data/pdf.sav' & radii3=radii & prob3=prob & time3=t
restore,file=dir4+'/data/pdf.sav' & radii4=radii & prob4=prob & time4=t
;
;----normalization------------
a0=10.e-6 
n0_Sw=1.e10
facSw=a0/n0_Sw
facSw=facSw/microm  ; f=tild{f}/r, here divide "microm" to compensate, XY. 
; The initial particle number density were 10^11 for the Smoluchowski
; simulation while it was set to 10^10 for the Swarm simulation. We
; must therefore rescale.
;
rescale=0.1

it=0
;dlnad=alog(radii1(3,it))-alog(radii1(2,it))
dlnad3=alog(radii3(3,it))-alog(radii3(2,it))
dlnad4=alog(radii4(3,it))-alog(radii4(2,it))
dlnad2=alog(radii2_0(3))-alog(radii2_0(2))
xx0=radii2_0
yy0=rescale*prob2_0*facSw/(dlnad2*xx0)
plot_oo,xx0,yy0,xr=xr,yr=yr,/nodata
oplot,xx0,yy0,col=55
;print,'np,t=',total(prob1(*,it)),time1(it)
;oplot,radii1(*,it)/microm,prob1(*,it)/dlnad,col=122
xx3=radii3(*,it)/microm
yy3=prob3(*,it)*facSw/(dlnad3*xx3)
;oplot,xx3,yy3,li=0,col=122,thick=thick1
xx4=radii4(*,it)/microm
yy4=prob4(*,it)*facSw/(dlnad4*xx4)
;oplot,xx4,yy4,li=2,col=155
oplot,xx4,yy4,li=0,col=122
;
xx1=radii1_0
yy1=rescale*prob1_0*facSw/(dlnad2*xx1)
oplot,xx1,yy1,li=1,col=55
;
;---------------------------------------------------
it=11
;dlnad=alog(radii1(3,it))-alog(radii1(2,it))
dlnad3=alog(radii3(3,it))-alog(radii3(2,it))
dlnad4=alog(radii4(3,it))-alog(radii4(2,it))
dlnad2=alog(radii2_1(3))-alog(radii2_1(2))
aa1=radii2_1
bb1=rescale*prob2_1*facSw/(dlnad2*aa1)
oplot,aa1,bb1,col=55
aa2=radii3(*,it)/microm
bb2=prob3(*,it)*facSw/(dlnad3*aa2)
;oplot,aa2,bb2,li=0,col=122,thick=thick1
aa3=radii4(*,it)/microm
bb3=prob4(*,it)*facSw/(dlnad4*aa3)
;oplot,aa3,bb3,li=2,col=155
oplot,aa3,bb3,li=0,col=122
;print,'np,t=',total(prob1(*,it)),time1(it)
;STOP
;
xx1=radii1_1
yy1=rescale*prob1_1*facSw/(dlnad2*xx1)
oplot,xx1,yy1,li=1,col=55
;
;---------------------------------------------------
it=21
;dlnad=alog(radii1(3,it))-alog(radii1(2,it))
dlnad3=alog(radii3(3,it))-alog(radii3(2,it))
dlnad4=alog(radii4(3,it))-alog(radii4(2,it))
dlnad2=alog(radii2_2(3))-alog(radii2_2(2))
cc0=radii2_2
dd0=rescale*prob2_2*facSw/(dlnad2*cc0)
oplot,cc0,dd0,col=55
;oplot,radii1(*,it)/microm,prob1(*,it)/dlnad,col=122
cc3=radii3(*,it)/microm
dd3=prob3(*,it)*facSw/(dlnad3*cc3)
;oplot,cc3,dd3,li=0,col=122,thick=thick1
cc4=radii4(*,it)/microm
dd4=prob4(*,it)*facSw/(dlnad4*cc4)
;oplot,cc4,dd4,li=2,col=155
oplot,cc4,dd4,li=0,col=122
;print,'np,t=',total(prob1(*,it)),time1(it)

;
xx1=radii1_2
yy1=rescale*prob1_2*facSw/(dlnad2*xx1)
oplot,xx1,yy1,li=1,col=55
;
;---------------------------------------------------
it=31
dlnad3=alog(radii3(3,it))-alog(radii3(2,it))
;dlnad=alog(radii1(3,it))-alog(radii1(2,it))
dlnad4=alog(radii4(3,it))-alog(radii4(2,it))
dlnad2=alog(radii2_3(3))-alog(radii2_3(2))
ee0=radii2_3
ff0=rescale*prob2_3*facSw/(dlnad2*ee0)
oplot,ee0,ff0,col=55
;oplot,radii1(*,it)/microm,prob1(*,it)/dlnad,col=122
ee3=radii3(*,it)/microm
ff3=prob3(*,it)*facSw/(dlnad3*ee3)
;oplot,ee3,ff3,li=0,col=122,thick=thick1
ee4=radii4(*,it)/microm
ff4=prob4(*,it)*facSw/(dlnad4*ee4)
;oplot,ee4,ff4,li=2,col=155
oplot,ee4,ff4,li=0,col=122
;
xx1=radii1_3
yy1=rescale*prob1_3*facSw/(dlnad2*xx1)
oplot,xx1,yy1,li=1,col=55
;
;---------------------------------------------------
;-----arrow-------
;----swarm----
a1=(total(ee4^1*ff4)/total(ff4))^(1./1.)
a3=(total(ee4^3*ff4)/total(ff4))^(1./3.)
a6=(total(ee4^6*ff4)/total(ff4))^(1./6.)
a12=(total(double(ee4)^12*ff4)/total(ff4))^(1./12.)
a24=(total(double(ee4)^24*ff4)/total(ff4))^(1./24.)
fact1=10. & fact2=1.
a=a1 & arrow,a,yr[0]*fact1,a,yr[0]*fact2,col=122,/data,thick=1
a=a3 & arrow,a,yr[0]*fact1,a,yr[0]*fact2,col=122,/data,thick=1
a=a6 & arrow,a,yr[0]*fact1,a,yr[0]*fact2,col=122,/data,thick=1
a=a12 & arrow,a,yr[0]*fact1,a,yr[0]*fact2,col=122,/data,thick=1
a=a24 & arrow,a,yr[0]*fact1,a,yr[0]*fact2,col=122,/data,thick=1
;----Euler----
a1=(total(ee0^1*ff0)/total(ff0))^(1./1.)
a3=(total(ee0^3*ff0)/total(ff0))^(1./3.)
a6=(total(ee0^6*ff0)/total(ff0))^(1./6.)
a12=(total(double(ee0)^12*ff0)/total(ff0))^(1./12.)
a24=(total(double(ee0)^24*ff0)/total(ff0))^(1./24.)
fact1=.1 & fact2=1. & j=1
a=a1 & arrow,a,yr[j]*fact1,a,yr[j]*fact2,col=55,/data,thick=1
a=a3 & arrow,a,yr[j]*fact1,a,yr[j]*fact2,col=55,/data,thick=1
a=a6 & arrow,a,yr[j]*fact1,a,yr[j]*fact2,col=55,/data,thick=1
a=a12 & arrow,a,yr[j]*fact1,a,yr[j]*fact2,col=55,/data,thick=1
a=a24 & arrow,a,yr[j]*fact1,a,yr[j]*fact2,col=55,/data,thick=1
;-----arrow-------
; Legends
xx=200
dx=20
dx=100
yy=1.e-4
dy=0.2
dy=0.1
legend,xx,dx,yy*dy^0,0,'Euler',siz=size2,col=55
legend,xx,dx,yy*dy^1,0,'!8N!6!dp!n/!8N!6!dgrid!n=32',col=122,siz=size2
xyouts,10.,3.e-9,'!8a!6!d1!n',col=122
xyouts,40.,1.e-1,'!8a!6!d3!n',col=55
xyouts,200.,1.e-1,'!8a!6!d6!n',col=55
xyouts,500.,1.e-1,'!8a!6!d12!n',col=55
xyouts,1100.,1.e-1,'!8a!6!d24!n',col=55
;---------------------------------------------------a24-------------------------------------------------------------
; Define directories
;
dir01g='$PENCIL_HOME/xiangyu/turb_coag/grav/0D/nocondens3457coag_a1em5_sig02_s001a_fast_a15em6_1p00543/'
;dir01='$PENCIL_HOME/xiangyu/turb_coag/grav/0D/nocondens3457coag_a1em5_sig02_s001a_fast_a15em6_1p00543_MC2/'
;dir2='../SW32_coag_grav1_ap10um_shima_65kp_16procs/'
;dir3='../SW32_coag_grav1_ap10um_shima_131kp_32procs/'
;dir4='../SW32_coag_grav1_ap10um_shima_262kp_32procs/'
dir4='../SW32_coag_grav1_ap10um_shima_1048kp_32procs/'
;standard swarm model-----------------
;dirx1='../SW32_coag_grav1_ap10um_shima_65kp_16procs_standard/'
;dirx2='../SW32_coag_grav1_ap10um_shima_131kp_32procs_standard/'
;dirx3='../SW32_coag_grav1_ap10um_shima_262kp_32procs_standard/'
;
; Read in data
;
pc_read_ts,obj=ts01g,datadir=dir01g+'data/',/double
;pc_read_ts,obj=ts01,datadir=dir01+'data/',/double
;pc_read_ts,obj=ts2,datadir=dir2+'data/',/double
;pc_read_ts,obj=ts3,datadir=dir3+'data/',/double
pc_read_ts,obj=ts4,datadir=dir4+'data/',/double

;pc_read_ts,obj=tsx1,datadir=dirx1+'data/',/double
;pc_read_ts,obj=tsx2,datadir=dirx2+'data/',/double
;pc_read_ts,obj=tsx3,datadir=dirx3+'data/',/double
;
; Find scaling parameter for the time
;
n0=1e8
n=1e10
tscale=n/n0
;
thick1=5
thick2=3
tilde='!9!s!aA!n!r!6'
bar='!20!s!A$!n!r!6'
;
mm=1.e-3  ; millimeter
!x.title=tilde+'!8t!6 [s]'
;---------------------------a24-----------------
a24_01g=(ts01g.admom24/ts01g.admom0)^(1./24)
a24_4=(ts4.admom24/ts4.admom0)^(1./24)

plot_io,ts01g.t*tscale*10.,(ts01g.admom24/ts01g.admom0)^(1./24)/mm,yr=[1e-2,3.0],xr=[0,40*tscale],$
        ytit='!8a!6!d24!n [mm]',thick=thick1,li=0,xticks=2,/nodata
oplot,ts01g.t*tscale*10.,(ts01g.admom24/ts01g.admom0)^(1./24)/mm,col=55,thick=6
oplot,ts4.t*tscale,(ts4.admom24/ts4.admom0)^(1./24)/mm,li=0,col=122,thick=thick1
;
;------------a12--------------
a12_01g=(ts01g.admom12/ts01g.admom0)^(1./12)
a12_4=(ts4.admom12/ts4.admom0)^(1./12)

oplot,ts01g.t*tscale*10.,(ts01g.admom12/ts01g.admom0)^(1./12)/mm,col=55,thick=6,li=1
oplot,ts4.t*tscale,(ts4.admom12/ts4.admom0)^(1./12)/mm,col=122,thick=thick1,li=1
;
;------------a6--------------
a6_01g=(ts01g.admom6/ts01g.admom0)^(1./6)
a6_4=(ts4.admom6/ts4.admom0)^(1./6)

oplot,ts01g.t*tscale*10.,(ts01g.admom6/ts01g.admom0)^(1./6)/mm,col=55,thick=6,li=2
oplot,ts4.t*tscale,(ts4.admom6/ts4.admom0)^(1./6)/mm,col=122,thick=thick1,li=2
;------------a3--------------
a3_01g=(ts01g.admom3/ts01g.admom0)^(1./3)
a3_4=(ts4.admom3/ts4.admom0)^(1./3)

oplot,ts01g.t*tscale*10.,(ts01g.admom3/ts01g.admom0)^(1./3)/mm,col=55,thick=6,li=3
oplot,ts4.t*tscale,(ts4.admom3/ts4.admom0)^(1./3)/mm,col=122,thick=thick1,li=3
;
;;------------a1--------------
;a1_01g=(ts01g.admom1/ts01g.admom0)^(1./1)
;a1_4=(ts4.admom1/ts4.admom0)^(1./1)
;
;oplot,ts01g.t*tscale*10.,(ts01g.admom1/ts01g.admom0)^(1./1)/mm,col=55,thick=6,li=4
;oplot,ts4.t*tscale,(ts4.admom1/ts4.admom0)^(1./1)/mm,col=122,thick=thick1,li=4
;
xyouts,1.3e3,2.e-2,'!8a!6!d3!n',col=55
xyouts,1.8e3,1.e-1,'!8a!6!d6!n',col=55
xyouts,1.8e3,2.e-1,'!8a!6!d12!n',col=55
xyouts,3.e3,2.,'!8a!6!d24!n',col=55

;-------linear fit for different moments------
result24_01g=linfit(ts01g.t,alog(a24_01g))
result24_4=linfit(ts4.t,alog(a24_4))

result12_01g=linfit(ts01g.t,alog(a12_01g))
result12_4=linfit(ts4.t,alog(a12_4))

result6_01g=linfit(ts01g.t,alog(a6_01g))
result6_4=linfit(ts4.t,alog(a6_4))

result3_01g=linfit(ts01g.t,alog(a3_01g))
result3_4=linfit(ts4.t,alog(a3_4))
;----slope as the inset-----
!p.position=[.645,.7,.75,.96]
!p.charsize=.8
!x.title='!7z'
!y.title='!7g'
moments=[3,6,12,24]
slope_swarm=[result3_4(1),result6_4(1),result12_4(1),result24_4(1)]
slope_Euler=[result3_01g(1),result6_01g(1),result12_01g(1),result24_01g(1)]
plot,moments,slope_swarm,/nodata,/noerase,yticks=3
oplot,moments,slope_swarm,col=122
;oplot,moments,slope_Euler,col=55
!p.position=0
;
print,'mv idl.ps ~/tex/xiangyu/SwarmSmolu_numerics/fig/moments_shima_Nx32_a24.ps'
;
!p.multi=0
;
STOP
;----------------------------------radius--------------------------------------------------
;
!x.title=tilde+'!8t!6 [s]'
microm=1.e-6
yr=[1e-5,2e-4]/microm
;plot_io,ts01.t*tscale*10,ts01.admom1/ts01.admom0/microm,yr=yr,xr=[0,40*tscale],$
plot_io,ts01g.t*tscale*10,ts01g.admom1/ts01g.admom0/microm,yr=yr,xr=[0,40*tscale],$
        ytit=bar+'!8r !6[!7l!n!6m]',thick=thick1,li=0,xticks=2,/nodata
oplot,ts01g.t*tscale*10,ts01g.admom1/ts01g.admom0/microm,col=55,thick=6
;oplot,ts01.t*tscale*10,ts01.admom1/ts01.admom0/microm,col=55,li=1,thick=8
;oplot,ts2.t*tscale,ts2.admom1/ts2.admom0/microm,li=1,col=122,thick=thick1
;oplot,ts3.t*tscale,ts3.admom1/ts3.admom0/microm,li=2,col=122,thick=thick1
oplot,ts4.t*tscale,ts4.admom1/ts4.admom0/microm,li=0,col=122,thick=thick1

;standard swarm model-------------
;oplot,tsx1.t*tscale,tsx1.admom1/tsx1.admom0/microm,li=1,col=177,thick=thick1
;oplot,tsx2.t*tscale,tsx2.admom1/tsx2.admom0/microm,li=2,col=177,thick=thick1
oplot,tsx3.t*tscale,tsx3.admom1/tsx3.admom0/microm,li=0,col=177,thick=thick1

xx=300
dx=300
;yy=1.5e-4
yy=2.e-4/microm
;dy=0.82
dy=0.7
;legend,xx,dx,yy*dy^1,1,'!8N!6!dp!n/!8N!6!dgrid!n= 2',col=122,size=size2
;legend,xx,dx,yy*dy^2,2,'!8N!6!dp!n/!8N!6!dgrid!n= 4',col=122,size=size2
;legend,xx,dx,yy*dy^3,0,'!8N!6!dp!n/!8N!6!dgrid!n= 8',col=122,size=size2
;legend,xx,dx,yy*dy^4,0,'Euler',size=size2,col=55
xyouts,2000,120,'scheme I',col=177,size=size2
xyouts,1000,50,'scheme II',col=122,size=size2
xyouts,3000,17,'Euler',col=55,size=size2

;
;;---------a24 as inset--------------
;!p.charsize=size3
;;!p.position=[.52,.33,.628,.54]
;!p.position=[.53,.345,.628,.54]
;xr=[0.,4.e3]
;yr=[1e-2,3.0]
;!y.title='!8a!6!d24!n [mm]'
;plot_io,xr,yr,xtit=xtit,/nodata,/noerase,xticks=2,yticks=2
;oplot,ts2.t*tscale,(ts2.admom24/ts2.admom0)^(1./24)/mm,li=1,col=122,thick=thick1
;oplot,ts3.t*tscale,(ts3.admom24/ts3.admom0)^(1./24)/mm,li=2,col=122,thick=thick1
;oplot,ts4.t*tscale,(ts4.admom24/ts4.admom0)^(1./24)/mm,li=0,col=122,thick=thick1
;; standard swarm data--------
;oplot,tsx1.t*tscale,(tsx1.admom24/tsx1.admom0)^(1./24)/mm,li=1,col=177,thick=thick1
;oplot,tsx2.t*tscale,(tsx2.admom24/tsx2.admom0)^(1./24)/mm,li=2,col=177,thick=thick1
;oplot,tsx3.t*tscale,(tsx3.admom24/tsx3.admom0)^(1./24)/mm,li=0,col=177,thick=thick1
;!p.position=0
;---------------------------------


print,'mv idl.ps ~/tex/xiangyu/SwarmSmolu_numerics/fig/moments_shima_Nx32_a24.ps'
;
END