pro refractivity, t,lm,ph,x1,x2,er,ei
; Carey 10/05
;     This program calculates the complex dielectric function for ice
;     and water at precipitation and cloud radar wavelengths.
;     The equations are taken from Ray, P. S., 1972: Broadband complex
;     refractive indices of ice and water.  Applied Optics, 11(8), 1836-1844.

;     Input:  1. liquid water (ph=0) vs ice (ph=1)
;             2. Wavelength (cm)
;             3. Temperature (C):  strictly applicable for ice : -20C to 0C
;						       water: -20C to 50C

;     Output: Dielectric function (real,imaginary) and
;             Refractive Index(real, imaginary) 

;     first get pi

      pi = 4.D0*atan(1.D0)

;     set functions for water or ice respectively

      if(ph EQ 0)then begin
        epss=78.54*(1.0-4.579e-3*(t-25.0)+1.19e-5*(t-25.0)^2- $
            2.8e-8*(t-25.0)^3)

        sigm=12.5664e8

        lams=0.00033836*exp(2513.98/(t+273))

        alph=-16.8129/(t+273)+0.0609265

        einf=5.27137+0.0216474*t-0.00131198*t^2
      endif

      if(ph EQ 1)then begin
        epss=203.168+2.5*t+0.15*t^2

        sigm=1.26*exp(-12500.0/((t+273.0)*1.9869))

        lams=9.990288e-4*exp(13200.0/((t+273.0)*1.9869))

        alph=0.288+0.0052*t+0.00023*t^2

        einf=3.168
      endif

;     calculate the dielectric function eps = er + i*ei

;     real

      er=einf+((epss-einf)*(1+(lams/lm)^(1-alph)*sin(alph*pi/2))) $
        / (1+2*(lams/lm)^(1-alph)*sin(alph*pi/2)+(lams/lm)^ $
        (2*(1-alph)))


;     imaginary component

      ei=((epss-einf)*(lams/lm)^(1-alph)*cos(alph*pi/2))/ $
        (1+2*(lams/lm)^(1-alph)*sin(alph*pi/2)+(lams/lm)^ $
        (2*(1-alph)))+(sigm*lm)/18.8496e10 


; refractive index m=n+ik=(n,k) from dielectric e=er+iei=(er,ei)
; x1=n and x2=k

      x1=( ( er + (er^2+ei^2)^0.5 ) /2 )^0.5 
;      simple approximation for x1 only valid for n2>>k2 or er>>ei
;      x1=er^0.5

      x2=ei/(2*x1)

;      print, x1, x2, er, ei
       print, x1, x2
stop
end