function index_water,t_in,w_in

;-------------------------------------------------------------------
; Purpose:
;   this code computes the refractive indicies for liquid water 
;   based on the work of Ray, P., Broadband complex refractive.....
;   Applied Optics, 8, 1836-1844, 1972.
;
; input: 
;    t_in.......temperature (K)
;    w_in.......wavelength (microns)
; output:
;    index......complex index
;
; Source: Mark Hervig
;-------------------------------------------------------------------

      pi  = 3.14159
      t   = t_in-273.16    ; t in C
      wav = w_in
      dxr = 0.
      dxi = 0.

      if (t gt 50.0  or  t lt -20.0) then begin
        print,'Temperature must be between -20 and 50 C'
        return,0
      endif

;---- parameters from eqn. 7a-c, epsilon sub infinity, alpha, lambda sub s

      cpf   = 5.27137 + 0.021647*t - 0.00131198 * t^2
      alpha = -16.8129/(t+273.) + 0.0609265
      wavs  = 0.00033836 * exp(2513.98/(t+273.))
      cps   = 78.54*(1.-4.579e-3 *(t-25.) $
              +1.19e-5 *(t-25.)^2 -2.8e-8 *(t-25.)^3)

;---- complex permitivity, cp = cpr - i*cpi, Eqns 2,3

      term1 = wavs/(wav*1.e-4)
      cpr   = cpf + (cps-cpf)/(1.+term1^2)
      cpi   = (cps-cpf)*term1/(1.+term1^2)

;---- complex permitivity, cp = cpr - i*cpi, Eqns 5,6

      term1 = (wavs/(wav*1.e-4))^(1-alpha)
      term2 = (cps-cpf)*(1.+term1 * sin(alpha*pi/2.))
      term3 = 1.+2.*term1*sin(alpha*pi/2.)+term1^2
      cpr   = cpf + term2/term3

      term2 = (cps-cpf)*term1 * cos(alpha*pi/2.)
      cpi   = term2/term3+12.5664e8*wav*1.e-4/18.8496e10

;---- complex refractive index for specific wavelength regions

      tav = 0.0001 * (t-25.0) * exp((wav/4.0)^0.25)

      term1 = 3352.27^2 - (10000.0/wav)^2
      term2 = 1639.00^2 - (10000.0/wav)^2
      term3 = 588.240^2 - (10000.0/wav)^2
      term4 = 99.9140e4*term1/(term1^2+15.1963e4*(10000./wav)^2) $
            + 50.4835e3*term2/(term2^2+9246.2700*(10000./wav)^2) $
            + 84.4697e4*term3/(term3^2+10.7615e5*(10000./wav)^2)
      dxra  = (1. +tav) * (1.79907+term4)^0.5

      term1 = 1639.0^2 - (10000.0/wav)^2
      term2 = 688.24^2 - (10000.0/wav)^2
      term3 = 161.29^2 - (10000.0/wav)^2
      term4 = 52340.4*term1/(term1^2 +10399.2*(10000./wav)^2) $
            + 345005.*term2/(term2^2 +259913.*(10000./wav)^2) $
            + 43319.7*term3/(term3^2 +27661.2*(10000./wav)^2)
      dxrb  = (1. +tav) * (1.83899+term4)^0.5

      if (wav le 6.0) then $
        dxr = dxra $
      else if (wav gt 6.0 and wav le 7.0) then $
        dxr = dxra*(7.0-wav) + dxrb*(wav-6.0) $
      else if (wav gt 7.0 and wav le 340.0) then $
        dxr = dxrb $
      else begin
        dxr = -999.9
        print,'Wavelength out of range for real component, index_water.pro'
        goto,jump10
      endelse

;     dxid = sqrt(dxr^2 - cpr)
      dxid = cpi/(2.*dxr)
      
      if (wav le 2.97) then $
        dxi = dxid +0.27 *exp(-1.0*abs(alog(wav/2.97)/0.025)^2.0) $
                   +0.01 *exp(-1.0*abs(alog(wav/4.95)/0.06)^1.0)  $
      else if (wav gt 2.97 and wav le 4.95) then $ 
        dxi = dxid +0.27 *exp(-1.0*abs(alog(wav/2.97)/0.04)^2.0) $
                   +0.01 *exp(-1.0*abs(alog(wav/4.95)/0.06)^1.0) $
      else if (wav gt 4.95 and wav le 6.1) then $
        dxi = dxid +0.01 *exp(-1.0*abs(alog(wav/4.95)/0.05)^1.0) $
                   +0.12 *exp(-1.0*abs(alog(wav/6.10)/0.08)^2.0) $
      else if (wav gt 6.1 and wav le 17.0) then $ 
        dxi = dxid +0.12 *exp(-1.0*abs(alog(wav/6.10)/0.042)^0.6) $
                   +0.39 *exp(-1.0*abs(alog(wav/17.0)/0.165)^2.4) $
                   +0.41 *exp(-1.0*abs(alog(wav/62.0)/0.220)^1.8) $
      else if (wav gt 17.0 and wav le 62.0) then $ 
        dxi = dxid +0.39 *exp(-1.0*abs(alog(wav/17.0)/0.45)^1.3) $
                   +0.41 *exp(-1.0*abs(alog(wav/62.0)/0.22)^1.8) $
                   +0.25 *exp(-1.0*abs(alog(wav/300.)/0.40)^2.0) $
      else begin
        dxi = -999.9
        print,'Wavelength out of range for imaginary component, index_water.pro'
        goto,jump10
      endelse

      jump10:
      
      index = complex(dxr,dxi)

      return,index
      end