;$Id: gravity_test.pro,v 1.1 2016/06/29 13:08:16 xiangyu Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=6,yoffset=3
  device,xsize=18,ysize=10,yoffset=3
  !p.charthick=2 & !p.thick=2 & !x.thick=2 & !y.thick=2
endif else begin
   device,decompose=0
   loadct,5
end

;Swarm shima
;run1='SW256condens0_coag_grav0_turb_shima_f2em2_ngp_a24_2048CPU'
run1='SW16condens0_coag_grav_ap10um'
run2='SW16condens0_coag_grav2g_ap10um'

default,iread,0
if iread eq 0 then begin
  pc_read_ts,o=ts1,datadir='../'+run1+'/data',/double
  pc_read_ts,o=ts2,datadir='../'+run2+'/data',/double
  iread=1
endif

!p.charsize=1.1
!x.margin=[6.,2.]
!y.margin=[3.2,.5]
bar='!20!s!A$!n!r!6'
;
;  radius
;
s=1.
microm=1e-6
onethird=1./3.
rhow=1e3
a1=ts1.admom1/ts1.admom0
a2=ts2.admom1/ts2.admom0


; Find scaling parameter for the time
n0=1e8  ; 1e8 cloud droplets typically
tscale=n/n0
nevery=50
siz=0.8

!p.multi=[0,2,2]
xr=[0,6.e3]
yr=[10.,45]
tilde='!9!s!aA!n!r!6'
;!x.title=tilde+'!8t!6 [s]'
!x.title='!6'
!y.title='!8a!6!d1!n [!7l!6m]'
plot,every(s*ts1.t*tscale,nevery),every(a1/microm,nevery),xr=xr,yr=yr,/nodata
oplot,every(s*ts1.t*tscale,nevery),every(a1/microm,nevery),col=122,thick=5
oplot,every(s*ts2.t*tscale,nevery),every(a2/microm,nevery),li=2,col=55,thick=5
xyouts,2000,12.5,siz=siz,'g',col=122
xyouts,200,26,siz=siz,'2g',col=55

;----------third moment----------------
onethird=1./3.
a31=(ts1.admom3/ts1.admom0)^onethird
a32=(ts2.admom3/ts2.admom0)^onethird
;yr=[10.,35]
yr=[10.,400]
tilde='!9!s!aA!n!r!6'
;!x.title=tilde+'!8t!6 [s]'
!x.title='!6'
!y.title='!8a!6!d3!n [!7l!6m]'
plot,every(s*ts1.t*tscale,nevery),every(a31/microm,nevery),xr=xr,yr=yr,/nodata
oplot,every(s*ts1.t*tscale,nevery),every(a31/microm,nevery),col=122,thick=5
oplot,every(s*ts2.t*tscale,nevery),every(a32/microm,nevery),li=2,col=55,thick=5
;xyouts,2500,20.0,siz=siz,'!6Re=65',col=122
;xyouts,1000,27.0,siz=siz,'!6Re=660',col=55

;--------12th moments--------
a12_1=(ts1.admom12/ts1.admom0)^(1./12)
a12_2=(ts2.admom12/ts2.admom0)^(1./12)
mm=1.e-3  ; millimeter
print,minmax(a12_2/mm)
;yr=[10.,140]
;yr=[5.e-3,1.3e-1]
yr=[5.e-3,4.e0]
tilde='!9!s!aA!n!r!6'
!x.title=tilde+'!8t!6 [s]'
!y.title='!8a!6!d12!n [!6mm]'
plot,every(s*ts1.t*tscale,nevery),every(a12_1/mm,nevery),xr=xr,yr=yr,/nodata
oplot,every(s*ts1.t*tscale,nevery),every(a12_1/mm,nevery),col=122,thick=5
oplot,every(s*ts2.t*tscale,nevery),every(a12_2/mm,nevery),li=2,col=55,thick=5
;xyouts,2500,0.05,siz=siz,'!6Re=65',col=122
;xyouts,1000,0.07,siz=siz,'!6Re=660',col=55


;--------24th moments--------
a24_1=(ts1.admom24/ts1.admom0)^(1./24)
a24_2=(ts2.admom24/ts2.admom0)^(1./24)
mm=1.e-3  ; millimeter
;yr=[10.,140]
;yr=[5.e-3,2.e-1]
yr=[5.e-3,8.e0]
tilde='!9!s!aA!n!r!6'
!x.title=tilde+'!8t!6 [s]'
!y.title='!8a!6!d24!n [!6mm]'
plot,every(s*ts1.t*tscale,nevery),every(a24_1/mm,nevery),xr=xr,yr=yr,/nodata
oplot,every(s*ts1.t*tscale,nevery),every(a24_1/mm,nevery),col=122,thick=5
oplot,every(s*ts2.t*tscale,nevery),every(a24_2/mm,nevery),li=2,col=55,thick=5
;xyouts,2500,0.06,siz=siz,'!6Re=65',col=122
;xyouts,1000,0.12,siz=siz,'!6Re=660',col=55

!p.multi=0
print,'mv idl.ps ~/tex/nils/particles_coag/fig/gravity_test.ps'
;
END