pro test_dge

; puts the iwc and re onto the time-height grid.

ciret2_fname='c:/mace/temp/sgpmmcrC1ciret2Mace.c1.20000309.163000.cdf'

print, ciret2_fname
print, n_elements(ciret2_fname)

;for jjjkkk=0,n_elements(ciret2_fname)-1 do begin

cdfid=ncdf_open(ciret2_fname)	; open the netcdf file

print, 'opened ',ciret2_fname

time_did=ncdf_dimid(cdfid,'Time')
height_did=ncdf_dimid(cdfid,'Height')

ncdf_diminq, cdfid, time_did, char_strng, num_times
ncdf_diminq, cdfid, height_did, char_strng, num_heights

; get the id's of the variables to be read

num_sec_id=ncdf_varid(cdfid,'base_time')
time_offset_id=ncdf_varid(cdfid,'time_offset')
N0s_id=ncdf_varid(cdfid,'N0s')
mlambda_id=ncdf_varid(cdfid,'mlambda')
iwc_id=ncdf_varid(cdfid,'Ice_Water_Content')
height_id=ncdf_varid(cdfid,'Height')



; get the data



ncdf_varget, cdfid, num_sec_id, num_sec
ncdf_varget, cdfid, time_offset_id, time_offset
ncdf_varget, cdfid, height_id, height

ijday_ihrfrac_fm_numsec, num_sec, time_offset, hrfrac, jday, yy, mm, dd, hh, mi, ss


ncdf_varget, cdfid, iwc_id, iwc_in
ncdf_varget, cdfid, N0s_id, N0s_in
ncdf_varget, cdfid, mlambda_id, mlambda_in

;if max(radrat) gt 0.15 then iwc_in[where(radrat gt 0.1)]=-9999.

ncdf_close, cdfid

		dge=fltarr(num_times,num_heights)
		dge[*,*]=-9999.
		dge_alt=dge

; following loops attempts to solve equation 2.3 of Fu 1996.
count_it=0
dflag=1


print, 'calculating dge'
for jj=0,num_times-1 do begin
	for k=0,num_heights-1 do begin


	if N0s_in[k,jj] ne 0 and N0s_in[k,jj] gt -10000 then begin

			N0s=10.^((float(N0s_in[k,jj]))/1000.) & mlambda=10.^((float(mlambda_in[k,jj]))/1000.)


	;	aa=0.15775 & ba=-0.5371	; coeffecients for aspect ratio with D in mm
	;	term1=(aa*(3.+(2.*ba)))*(gamma((3.+(2.*ba))))/((abs(mlambda_in[k,jj]))^(1.+ba))
	;	term2=(2.+ba)*gamma(2.+ba)
	;	term3=((sqrt(3.))/4.)*aa*2.*ba*(gamma(2.*ba))/((abs(mlambda_in[k,jj]))^ba)
;
;		dge_alt[jj,k]=term1/(term2+term3)
;

		;aa=0.15775 & ba=-0.5371	; coeffecients for aspect ratio with D in mm
	 ;	term1=(aa*gamma(4.+(2.*ba)))/(mlambda^(1.+(2.*ba)))
	 ;	term2=(gamma(3.+ba))/(mlambda^((ba)))
	 ;	term3=(sqrt(3.)/4.)*aa*(gamma(3.+(2.*ba)))/(mlambda^((2.*ba)))

	;	aa=0.25 & ba=-0.3
	;	term1=aa*(gamma(4.+(2.*ba)))/((abs(mlambda))^(1.+(ba)))
	;	term2=gamma(3.+ba)
	;	term3=((sqrt(3.))/4.)*aa*(gamma(3.+(2.*ba)))/((abs(mlambda))^ba)



; first integral from size of 0 to 30 microns has aspect ratio of 1.
	asp=1.0 & lh=0. & rh=29.
	term3_30=(exp_integral_soln(3., -mlambda, lh, rh))
	term2_30=(exp_integral_soln(2., -mlambda, lh, rh))
	denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_30
	numer=asp*term3_30
	fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
	dge_30=fraction*(numer/denom)

	asp=0.8 & lh=30. & rh=79.
	term3_80=(exp_integral_soln(3., -mlambda, lh, rh))
	term2_80=(exp_integral_soln(2., -mlambda, lh, rh))
	denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_80
	numer=asp*term3_80
	fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
	dge_80=fraction*(numer/denom)


	asp=0.5 & lh=80. & rh=199.
	term3_200=(exp_integral_soln(3., -mlambda, lh, rh))
	term2_200=(exp_integral_soln(2., -mlambda, lh, rh))
	denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_200
	numer=asp*term3_200
	fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
	dge_200=fraction*(numer/denom)

	asp=0.34 & lh=200. & rh=499.
	term3_500=(exp_integral_soln(3., -mlambda, lh, rh))
	term2_500=(exp_integral_soln(2., -mlambda, lh, rh))
	denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_500
	numer=asp*term3_500
	fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
	dge_500=fraction*(numer/denom)

	asp=0.22 & lh=500. & rh=5000.
	term3_5000=(exp_integral_soln(3., -mlambda, lh, rh))
	term2_5000=(exp_integral_soln(2., -mlambda, lh, rh))
	denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_5000
	numer=asp*term3_5000
	fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
	dge_5000=fraction*(numer/denom)

	dge_alt=dge_30+dge_80+dge_200+dge_500+dge_5000

;
				dge[jj,k]=0. & numer=0. & denom=0.
				factor=1.0
			for j=5,30,5 do begin
				inte=(((float(j)*factor)^2)*float(j))*(5.*N0s*exp(-mlambda*float(j-2)))
				numer=numer+inte
				inte=(5.*N0s*exp(-mlambda*float(j-2)))*(((float(j-2)*factor)*float(j-2))+( ((float(j-2)*factor)^2)*((sqrt(3))/4.)   ))
				denom=denom+inte
			endfor
			;print, 'first ', numer, denom
				factor=0.8
			for j=31,79,3 do begin
				inte=(((float(j-1)*factor)^2)*float(j-1))*(2.*N0s*exp(-mlambda*float(j-1)))
				numer=numer+inte
				inte=(2.*N0s*exp(-mlambda*float(j-1)))*(((float(j-1)*factor)*float(j-1))+( ((float(j-1)*factor)^2)*((sqrt(3))/4.)   ))
				denom=denom+inte
			endfor
			;print, '2nd ', numer, denom
				factor=0.5
			for j=80,199,2 do begin
				inte=(((float(j)*factor)^2)*float(j))*(2.*N0s*exp(-mlambda*float(j)))
				numer=numer+inte
				inte=(2.*N0s*exp(-mlambda*float(j)))*(((float(j)*factor)*float(j))+( ((float(j)*factor)^2)*((sqrt(3))/4.)   ))
				denom=denom+inte
			endfor
			;print, '3rd ', numer, denom
				factor=0.34
			for j=199,500,2 do begin
				inte=(((float(j)*factor)^2)*float(j))*(2.*N0s*exp(-mlambda*float(j)))
				numer=numer+inte
				inte=(2.*N0s*exp(-mlambda*float(j)))*(((float(j)*factor)*float(j))+( ((float(j)*factor)^2)*((sqrt(3))/4.)   ))
				denom=denom+inte
			endfor
			;print, '4th ', numer, denom
			factor=0.22
			for j=500,2500,5 do begin
				inte=(((float(j-2)*factor)^2)*float(j-2))*(5.*N0s*exp(-mlambda*float(j-2)))
				numer=numer+inte
				inte=(5.*N0s*exp(-mlambda*float(j-2)))*(((float(j)*factor)*float(j-2))+( ((float(j-1)*factor)^2)*((sqrt(3))/4.)   ))
				denom=denom+inte
			endfor
			;print, 'fifth', numer, denom

			dge[jj,k]=numer/denom
			print, jj,k,dge[jj,k], dge_alt,mlambda,N0s, numer, denom

	endif else begin

		if N0s_in[k,jj] ne 0 and N0s_in[k,jj] gt -10000 then begin
			dge[jj,k]=50.
		endif

	endelse

	endfor
	;print, float(jj)/float(num_times)
		if count_it eq 100 then begin
			print, jj, num_times
			count_it=0
		endif
		count_it=count_it+1

endfor
print, 'done calculating dge'

end