;if !d.name eq 'PS' then begin
;  device,xsize=18,ysize=12,yoffset=3
;  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
;end
;
;!p.charsize=1.7
;!x.margin=[8.8,.5]
;!y.margin=[3.2,.5]
;!x.title='!8r!6 [m]'
;!y.title='!7C'
function gamma2014,D,mu,lambda,N0
N0=9.92e10
   return, N0*D^mu*exp(-lambda*D) 
end

mu=1.62
mu=-1.62
lam=1e-2
N0=9.92e10  ;total number of cloud droplets:w
D=(1+findgen(1000)); diameter
mydata=gamma2014(D,mu,lam)
print,mydata

xr=minmax(D)
yr=[1e11,1e14]
yr=minmax(mydata)
plot_oo,D,mydata,xr=xr,yr=yr
n=1 & D1=(total(mydata*D^n)/total(mydata))^(1./float(n))
n=3 & D3=(total(mydata*D^n)/total(mydata))^(1./float(n))
n=12 & D12=(total(1d0*mydata*D^n)/total(mydata))^(1./float(n))
n=24 & D24=(total(1d0*mydata*D^n)/total(mydata))^(1./float(n))
;
oo=[1,1]
oplot,oo*D1,yr,li=1
oplot,oo*D3,yr,li=2
oplot,oo*D12,yr,li=3
oplot,oo*D24,yr,li=0
;
STOP
;mu=1.62
;;lambda=1.05e4
;lambda=1.0
D=findgen(1000); diameter
D[0]=1
;----Define Gamma funciton------
;2014_Mcfarquhar
;----Initial guess of the parameters------
mu=1.62
lambda=1.0
;-------define a vector of weights-----
Y=D
weights = 1.0/Y
;-------Fit Gamma fucntion--------
yfit=CURVEFIT(D,Y,weights,mu,lambda,function_name='gamma2014')
print,'fit',mu,lambda

END 
;plot_oo,D,fr,yr=[1.e3,5.e10]