@mie_log_nrm
pro mie_spec_snoe

;This program finds the scattering coefficient as a function of wavelength
;using the index of refraction values from Shettle. I have compared his index values
;with those given by Warren, and they agree very well (to three or four places)
;distributions
;This program then finds the Mie scattering coefficient for backscattering as a
;function of wavelength from 0.2 to 1.0 micrometers.
;The first two pairs of plots show the index of refraction vs wavelength (real +imaginary)
;The second plot shows the backscattering coefficient for the Jensen model C, and particle
;density eaual to 100 cm^-3. Overplotted is the Raleigh scattering coefficient for comparison.
;The program is easily changed to accomodate other scattering angles

          ;Input data from S. Warren of indices vs. wavelength
          ;up to .3 microns, and from .3 to 1.0 micron the Shettle values
          ;lm is the wavelength in microns
lm=[.19,.191,.193,.195,.21,.25,.3,.337,.4,.488,.515,.55,.633,.694,.8,.86, 1.06] ;lm is in microns
mr=[1.4101,1.4081,1.4042,1.4007,1.38,1.3509,1.334,1.326,1.32,1.313,1.312,1.311,1.308,$
    1.306,1.304,1.303,1.3]
mi=[1.9E-8,1.75E-8,1.64E-8,1.59E-8,1.325E-8,8.623E-9,5.50e-9,4.50e-9,2.71e-9,1.75e-9,$
    2.19e-9,3.11e-9,1.09e-8,2.62e-8,8.e-8,2.15e-7,1.96e-6 ]
;mr is the real part of the index of refraction for various wavelengths lm
;mi is the imaginary part of the index of refraction

;We now define a wavelength vector of higher resolution in order to
;plot a smooth curve and also interpolate to an arbitrary wavelength
lam=[.215, .2365]  ;0.2+.02*findgen(45)
y2a=nr_spline(lm,mr)    ;calculates interpolating cubic spline coefficients
mrfit=nr_splint(lm,mr,y2a,lam)  ;Cubic spline fit for the wavelength array lam
oplot,lam,mrfit
;Repeat for imaginary index of refraction
;plot the index of refraction versus wavelength.
;plot,/ylog,lm,mi,psym=4,/ynozero,xtitle='wavelength',ytitle='imaginary index',$
;title='IMAGINARY INDEX OF REFRACTION',xrange=[0.2,1.0]
;Now find the spline fit to the imaginary index of refraction
y2=nr_spline(lm,mi) ;calculates interpolating cubic spline coefficients
mifit=nr_splint(lm,mi,y2,lam)
;oplot,lam,mifit
;stop,''
;
;Now compute scattered-light albedo spectrum for a number of angles
;The angles calculated are specified by NANG, the input to the Mie
;scattering program developed by Mike Callan. The actual angles are
;given by the formula
;ANG=180.*FINDGEN(NANG)/[2*(NANG-1)].
;for example, every 5 degrees would require NANG=19, yielding a phase function
;having the dimensions 2*(NANG-1)
;every 90 degrees requires NANG=2, etc.
;ANG=180.*[0,1,2,3,4, ETC]/36=[0,5,10,ETC]
;The smallest number (to save computer time) is NANG=2, so that the angles are [0,90,180].
;To truncate the integration over particle size at an upper limit, you can add the
;optional keyword MAX_RAD at the end, for example, MAX_RAD=10 would truncate
;the integration at 10 microns.
;Jensen's four particle size models, A, B, C, and D:
;Model A; r0=22 nm, sigma=1.35
;Model B; r0=30 nm, sigma=1.42
;Model C: r0=45 nm, sigma=1.48
;Model D; r0=48 nm, sigma=1.65

nmod=50

;rmod=[.022,.03,.045,.048]
;sigma=[1.35,1.42,1.48,1.65]
rmod=findgen(nmod)*.002 + .01
sigma=fltarr(nmod)+1.4

NANG=46
ang=180.* findgen(2*nang-1)/(2*(nang-1))

kext=fltarr(2,nmod,2*nang-1)
ksca=kext
p=fltarr(3,4,2)    ;phase function. First index=angle,
;2nd index=particle size model, 3rd index= wavelength

;Begin calculation of spectrum

phase=fltarr(2,nmod,2*nang-1)

for i=0,nmod-1 do begin  ;outer loop, particle size model number
    for j=0,1 do begin    ;inner loop, wavelength
    mie_log_nrm,lam(j),complex(mrfit(j),mifit(j)),NANG,rmod(i),sigma(i),ke,ks,p1

    for k=0,2*nang-2 do begin
        phase(j,i,k)=p1(k)
        kext(j,i,k)=ke
        ksca(j,i,k)=ks
        endfor
;    p(0,i,j)=p1(0) ;phase function at zero degrees scattering angle
;    p(1,i,j)=p1(1) ;phase function at 90 degrees scattering angle
;    p(2,i,j)=p1(2) ;phase function at 180 degrees scattering angle
    endfor
endfor
;
convf=1.E-8 ;conversion factor from sq microns to sq cm
;
kext=kext*convf ;kext is the extinction cross-section in sq.cm
ksca=ksca*convf ;kscais the scattering cross-section in sq.cm

save,file='snoe_mie.save',ang,nmod,rmod,sigma,phase,kext,ksca

;plot scattered spectrum at 180 degrees angle for particle model C
;!p.multi=0
;assume number density of PMC particles in cm^-3
;npmc=100.
;plot,/ylog,lam,npmc*ksca(2,*)*p(2,2,*)/(4.*!pi),xtitle='Wavelength (micrometers)',$
;ytitle='Scattering Coefficient (cm^-1)'
;Now calculate the Rayleigh scattering coefficient from Eric Shettle's program, ray_scat
;The input wavelength is in nm. kray is the Rayleigh scattering coefficient, computed from
;the Rayleigh cross-section (in cm^2), and a model atmosphere for the cold arctic summertime.
; ray_coef is the scattering coefficient at the altitude (2nd input), in km^-1.
;nden is the number density cm^-3.
;evaluate scattering coefficient at 83 km
;th=180. ;backscattered angle for comparing with lidar
;prb=0.75*(1.+cos(th*!pi/180.)^2)    ;backscattering phase function for Rayleigh scattering
;ray_scat,lam*1.E3,83.,ray_coef,kray,nden
;oplot,lam,1.E-3*ray_coef*(prb/(4.*!pi)),line=1  ;convert to cm-1
stop,''
end