function index_h2so4,wl,temp,wts,iread

;-------------------------------------------------------------------
;
;      purpose
;
;       This routine returns the refractive
;     index of a sulphuric acid aerosol as a function of 
;     % H2SO4 composition, temperature, and wavelength.
;
;       The refractive indicies for 300 K and a range of
;     wavelengths and compositions were taken from
;     Palmer and Williams. Interpolations in wavelength 
;     and composition are done.
;
;       Lorentz-Lorenz temperature corrections are used.
;
;       A file containing the refractive indicies for the desired
;     wavelengths and compositions as given in Palmer and 
;     Williams must be supplied.
;
;     source:  Mark Hervig,  U of Wyoming
;
; input
;     wl........wavelength in microns
;     wts.......wt. fraction h2so4
;     temp......temperature in K
;     iread.....0=read the data file,  1=take from saved
; output
;     rindex....complex refractive index
;----------------------------------------------------------------------

      ws = [.25,.38,.50,.75,.845,.956]
      ns = n_elements(ws)
      xr = fltarr(ns)
      xi = fltarr(ns)

      header1 = ' '

      pat = '/users/mhervig/data/refractive_index/'
      
      file = pat + 'index_H2SO4-300K_palmer.dat'
      
;---- read in refractive index file,  if iread=0
;     otherwise the data is saved and then restored

;      if (iread eq 0) then begin
      
;       print,'open and read h2so4 indices'

        openr,lun,file,/get_lun
	
        readf,lun,format=('(a64)'),header1
        readf,lun, nwav
        readf,lun,format=('(a64)'),header1
        datr = fltarr(7,nwav)
        dati = fltarr(7,nwav)
        readf,lun,datr
        readf,lun,format=('(a64)'),header1
        readf,lun,dati

        dxr_25=datr(1,*) & dxr_38=datr(2,*) & dxr_50=datr(3,*) 
        dxr_75=datr(4,*) & dxr_85=datr(5,*) & dxr_95=datr(6,*)
        dxi_25=dati(1,*) & dxi_38=dati(2,*) & dxi_50=dati(3,*) 
        dxi_75=dati(4,*) & dxi_85=dati(5,*) & dxi_95=dati(6,*)

        wavel = 10000.0/datr(0,*)
        nw    = nwav
	
        close,lun 
	free_lun,lun

;       save,dxr_25,dxr_38,dxr_50,dxr_75,dxr_85,dxr_95, $
;            dxi_25,dxi_38,dxi_50,dxi_75,dxi_85,dxi_95, $
;            wavel,nwav

;      endif

;     if (iread ne 0) then restore

;---- interpolate index in wave # for each wt%

      dwdw = 0.0

      for i=0,nwav-2 do begin
        if (wavel(i) ge wl and wavel(i+1) lt wl) then begin
          dwdw   = (wl-wavel(i+1))/(wavel(i)-wavel(i+1))
          xr(5)  = dwdw *(dxr_95(i)-dxr_95(i+1)) +dxr_95(i+1)
          xi(5)  = dwdw *(dxi_95(i)-dxi_95(i+1)) +dxi_95(i+1)
          xr(4)  = dwdw *(dxr_85(i)-dxr_85(i+1)) +dxr_85(i+1)
          xi(4)  = dwdw *(dxi_85(i)-dxi_85(i+1)) +dxi_85(i+1)
          xr(3)  = dwdw *(dxr_75(i)-dxr_75(i+1)) +dxr_75(i+1)
          xi(3)  = dwdw *(dxi_75(i)-dxi_75(i+1)) +dxi_75(i+1)
          xr(2)  = dwdw *(dxr_50(i)-dxr_50(i+1)) +dxr_50(i+1)
          xi(2)  = dwdw *(dxi_50(i)-dxi_50(i+1)) +dxi_50(i+1)
          xr(1)  = dwdw *(dxr_38(i)-dxr_38(i+1)) +dxr_38(i+1)
          xi(1)  = dwdw *(dxi_38(i)-dxi_38(i+1)) +dxi_38(i+1)
          xr(0)  = dwdw *(dxr_25(i)-dxr_25(i+1)) +dxr_25(i+1)
          xi(0)  = dwdw *(dxi_25(i)-dxi_25(i+1)) +dxi_25(i+1)
          goto,jump1
        endif
      endfor

      jump1:
      if (dwdw eq 0.0) then begin
        print,'wavelength out of range in index_h2o_h2so4.pro'
        rindex = complex(0.,0.)
        result = rindex
        return,result
      endif 

;---- interpolate refractive index to wts 

      dxr = interpol(xr,ws,wts)
      dxi = interpol(xi,ws,wts)

      indexPW = complex(dxr,dxi)

;---- find the solution densities (g/cm3) for temp and wt% H2SO4

      wtn   = 1.e-9                            ; wt. fraction hno3
      rhoPW = lta_density(wts,wtn,300.0)  ; P+W measured @ 300K
      rho   = lta_density(wts,wtn,temp)

;---- apply Lorentz - Lorenz correction for temperature

      rindex = lorentz(rhoPW,rho,indexPW)

      result = rindex
      return,result
      end