@mie_log_nrm
pro mie_spec
;device,filename='index_of_refraction.ps',/landscape
;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
!p.font=-1       ;LOAD HERSHEY (VECTOR FONTS)
;!p.charsize=1.0       ;DEFAULT SIZE=1.0
!p.charsize=1.3
;loadct,39       ;LOAD COLOR TABLE 39 (RAINBOW +WHITE)
          ;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
!p.multi=[0,2,1]    ;2nd element=no of plots in horizontal
         ;3rd element=no of plots in vertical
;PLOT1- REAL INDEX OF REFRACTION
plot,lm,mr,psym=2,yrange=[1.3,1.4],xrange=[0.2,1.0],xtitle='wavelength',$
ytitle='real index',title='REAL INDEX';,background=255;,color=0
    ;coeffr=poly_fit(lm,mr,2,yfit,errfit)
;We now define a wavelength vector of higher resolution in order to
;plot a smooth curve and also interpolate to an arbitrary wavelength
lam=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
rmod=[.022,.03,.045,.048]
sigma=[1.35,1.42,1.48,1.65]
kext=fltarr(4,45)
ksca=kext
p=fltarr(3,4,45)    ;phase function. First index=angle,
;2nd index=particle size model, 3rd index= wavelength
NANG=2
;Begin calculation of spectrum
for i=0,3 do begin  ;outer loop, particle size model number
    for j=0,44 do begin    ;inner loop, wavelength
    mie_log_nrm,lam(j),complex(mrfit(j),mifit(j)),NANG,rmod(i),sigma(i),ke,ks,p1
    kext(i,j)=ke
    ksca(i,j)=ks
    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
;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