;-------Gamma distribution
function gamma_fit, shape, slope, n, r
  return, n*r^shape*exp(-slope*r)*slope^(shape+1.)/gamma(shape+1)
end

size=2.
size2=2.
if !d.name eq 'PS' then begin
  device,xsize=24,ysize=18,yoffset=3
  !p.charthick=3 & !p.thick=5 & !x.thick=3 & !y.thick=3
  !p.charsize=size
endif else begin
   !p.charsize=size
end
!x.margin=[7.5,.5]
!y.margin=[3.5,.5]
;
;---------size spectra of the Eulerian model and Gamma fitting-------------------------

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=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------------
!x.title='!8r!6 [!7l!6m]'
!y.title='!8f!6 !8r!6!dini!n/!8n!6!d0!n'

xr=[4,2500]  ;XY
yr=[1.e-10,3.e0]
microm=1e-6
thick1=6
n0_Sw=1.e10
a0=10.e-6 
facSw=a0/n0_Sw
microm=1.e-6
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
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
dr=xx0(3)-xx0(2)
yy0=rescale*prob2_0*facSw/(dlnad2*xx0)
plot_oo,xx0,yy0,xr=xr,yr=yr,/nodata  ;XY
oplot,xx0,yy0,col=55
;oplot,xx0,prob2_0,col=55
xx3=radii3(*,it)/microm
yy3=prob3(*,it)*facSw/(dlnad3*xx3)
xx4=radii4(*,it)/microm
yy4=prob4(*,it)*facSw/(dlnad4*xx4)
;;---calculate slope and shape------------------
a1=(total(double(xx0)^1*yy0)/total(yy0))^(1./1.)
a24=(total(double(xx0)^24*yy0)/total(yy0))^(1./24.)
;cutoff=where(xx0 ge a1 and xx0 le a24)
cutoff=where(xx0 eq xx0)
xx0_cutoff=xx0(cutoff)
save,file='xx0_cutoff.sav',xx0_cutoff
yy0_cutoff=yy0(cutoff)
m0=total(xx0_cutoff^0*prob2_0)
m1=total(xx0_cutoff^1*prob2_0)
m2=total(xx0_cutoff^2*prob2_0)
m3=total(xx0_cutoff^3*prob2_0)
m6=total(xx0_cutoff^6*prob2_0)
a1=m1/m0
a2=(m2/m0)^.5
shape=(a2^2-2*a1^2)/(a1^2-a2^2)
slope=(shape+1)/a1
n=m0
print,'t=0s, m6*m0/m3^2=', m6*m0/m3^2
print,'a1=',a1
print,'a2=',a2
print,'shape=',shape
print,'slope=',slope
print,'n=',n
;xx0_min=min(xx0_cutoff)
;xx0_max=max(xx0_cutoff)
;xx0_cutoff=grange(xx0_min,xx0_max,n_elements(xx0_cutoff))
gamma_p=gamma_fit(shape,slope,n,xx0_cutoff)
;gamma_p=rescale*gamma_p*facSw/(dlnad2*xx0)
gamma_p=rescale*gamma_p*facSw
oplot,xx0_cutoff,gamma_p,li=2,thick=6
;
it=11
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)
aa3=radii4(*,it)/microm
bb3=prob4(*,it)*facSw/(dlnad4*aa3)
;;---calculate slope and shape------------------
a1=(total(double(aa1)^1*bb1)/total(bb1))^(1./1.)
a24=(total(double(aa1)^24*bb1)/total(bb1))^(1./24.)
;cutoff=where(aa1 ge a1 and aa1 le a24)
cutoff=where(aa1 eq aa1)
aa1_cutoff=aa1(cutoff)
bb1_cutoff=bb1(cutoff)
m0=total(aa1_cutoff^0*prob2_1)
m1=total(aa1_cutoff^1*prob2_1)
m2=total(aa1_cutoff^2*prob2_1)
m3=total(aa1_cutoff^3*prob2_1)
m6=total(aa1_cutoff^6*prob2_1)
a1=m1/m0
a2=(m2/m0)^.5
shape=(a2^2-2*a1^2)/(a1^2-a2^2)
slope=(shape+1)/a1
n=m0
print,'t=1000s, m6*m0/m3^2=', m6*m0/m3^2
print,'a1=',a1
print,'a2=',a2
print,'shape=',shape
print,'slope=',slope
print,'n=',n
;aa1_min=min(aa1_cutoff)
;aa1_max=max(aa1_cutoff)
;aa1_cutoff=grange(aa1_min,aa1_max,n_elements(aa1_cutoff))
gamma_p=gamma_fit(shape,slope,n,aa1_cutoff)
;gamma_p=rescale*gamma_p*facSw/(dlnad2*aa1)
gamma_p=rescale*gamma_p*facSw
oplot,aa1,gamma_p,li=2,thick=6
;
it=21
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
cc3=radii3(*,it)/microm
dd3=prob3(*,it)*facSw/(dlnad3*cc3)
cc4=radii4(*,it)/microm
dd4=prob4(*,it)*facSw/(dlnad4*cc4)
;;---calculate slope and shape------------------
a1=(total(double(cc0)^1*dd0)/total(dd0))^(1./1.)
a24=(total(double(cc0)^24*dd0)/total(dd0))^(1./24.)
cutoff=where(cc0 eq cc0)
cc0_cutoff=cc0(cutoff)
dd0_cutoff=dd0(cutoff)
m0=total(cc0_cutoff^0*prob2_2)
m1=total(cc0_cutoff^1*prob2_2)
m2=total(cc0_cutoff^2*prob2_2)
m3=total(cc0_cutoff^3*prob2_2)
m6=total(cc0_cutoff^6*prob2_2)
a1=m1/m0
a2=(m2/m0)^.5
shape=(a2^2-2*a1^2)/(a1^2-a2^2)
slope=(shape+1)/a1
n=m0
print,'t=2000s, m6*m0/m3^2=', m6*m0/m3^2
print,'a1=',a1
print,'a2=',a2
print,'shape=',shape
print,'slope=',slope
print,'n=',n
;cc0_min=min(cc0_cutoff)
;cc0_max=max(cc0_cutoff)
;cc0_cutoff=grange(cc0_min,cc0_max,n_elements(cc0_cutoff))
gamma_p=gamma_fit(shape,slope,n,cc0_cutoff)
;gamma_p=rescale*gamma_p*facSw/(dlnad2*cc0_cutoff)
gamma_p=rescale*gamma_p*facSw
oplot,cc0_cutoff,gamma_p,li=2,thick=6
;
it=31
dlnad3=alog(radii3(3,it))-alog(radii3(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
ee3=radii3(*,it)/microm
ff3=prob3(*,it)*facSw/(dlnad3*ee3)
ee4=radii4(*,it)/microm
ff4=prob4(*,it)*facSw/(dlnad4*ee4)
;-----arrow at t=3000s-------
;----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
;
;---calculate slope and shape------------------
;cutoff=where(ee0 ge a1 and ee0 le a24)
cutoff=where(ee0 eq ee0)
ee0_cutoff=ee0(cutoff)
ff0_cutoff=ff0(cutoff)
save,file='../python/ee0.sav',ee0
save,file='../python/ff0.sav',ff0
m0=total(ee0_cutoff^0*prob2_3)
m1=total(ee0_cutoff^1*prob2_3)
m2=total(ee0_cutoff^2*prob2_3)
m3=total(ee0_cutoff^3*prob2_3)
m6=total(ee0_cutoff^6*prob2_3)
a1=m1/m0
a2=(m2/m0)^.5
shape=(a2^2-2*a1^2)/(a1^2-a2^2)
slope=(shape+1)/a1
n=m0
print,'t=2000s, m6*m0/m3^2=', m6*m0/m3^2
print,'a1=',a1
print,'a2=',a2
print,'shape=',shape
print,'slope=',slope
print,'n=',n
;ee0_min=min(ee0_cutoff)
;ee0_max=max(ee0_cutoff)
;ee0_cutoff=grange(ee0_min,ee0_max,n_elements(ee0_cutoff))
gamma_p=gamma_fit(shape,slope,n,ee0_cutoff)
;gamma_p=rescale*gamma_p*facSw/(dlnad2*ee0_cutoff)
gamma_p=rescale*gamma_p*facSw
oplot,ee0_cutoff,gamma_p,li=2,thick=6

;xyouts,100.,.5e-2,'!6Gamma'
;-----arrow-------
; Legends
xx=400
dx=20
dx=150
;yy=1.5e-4
yy=1.e-2
dy=0.2
dy=0.1
legend,xx,dx,yy*dy^0,2,'Gamma',siz=size2
legend,xx,dx,yy*dy^1,0,'Euler',siz=size2,col=55
;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


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