PRO multifilt, wavelength, trans, elem, ele_thick, mesh_trans,$
               wave_max=wave_max,wave_min=wave_min, ang=ang,$
               file_name=file_name, densities=densities, $
               n_files=n_files, caldir=caldir

; =========================================================================
;+
;
;NAME: multifilt.pro
;PURPOSE: Calculates the transmission of a multilayer filter
;CALLING SEQUENCE: multifilt, wavelength, trans		
;OPTIONAL KEYWORD INPUT: 
;         elem - names of elements in layers (e.g. 'Al')
;         elem_thickness - thickness of layer in Angstroms
;         mesh_trans - transmission of mesh backing
;         wave_max - max wavelength in Angstroms
;         wave_min - min wavelength in Angstroms
;         ang - incident angle in radians
;  
;DEFAULT KEYWORDS: defaults are for DAG filter on SXT on Yokkoh
;OUTPUT:
;         wavelength - range of wavelengths in angstroms from
;                      wave_min to wave_max
;         trans - transmission as a function of wavelength
;         densities - material densities
;         n_files - material index filenames as read-in
;
;ASSUMPTIONS: medium on both sides of filter is
;             vacuum, index of refraction data is in columns of the
;             form wvaelength,delta,beta (where n=(1-delta) +i(beta)), 
;             files with index data are named n_elem.txt (e.g. n_Al.txt)
;
;
;REFERENCES: G.R Fowles, Intro to Modern Optics, pp96-101
;
;HISTORY:
;         written 09-July-99 by K.Reeves
;         modified 14-Mar-05 by M.Weber --- uses refraction files
;           which have a "wavelength [nm]" column instead
;           of "energy [eV]".
;         modified 14-May-07 by M.Weber --- added densities and n_files.
;  progver = 'v2008.Sep.05' ;--- (A.Davey) Added alternate path for 
;                                $XRT_CALIB, using CALDIR keyword.
   progver = 'v2008.Oct.01' ;--- (M.Weber) Reviewed version. Updated so
;                                that it can handle both the filters and
;                                the contamination layer.
;
;-
; =========================================================================


  
  XRT_CALIB = get_logenv('$SSW_XRT_CALIB')
  if keyword_set(caldir) then XRT_CALIB = caldir

  IF not keyword_set(mesh_trans) THEN mesh_trans = 1.0
  IF not keyword_set(wave_max) THEN wave_max = 400. ;Angstroms
  IF not keyword_set(wave_min) THEN wave_min = 1.   ;Angstroms
  IF not keyword_set(ang) THEN ang = 0.

;check that the same number of elements and thicknesses are entered
  IF n_elements(elem) ne n_elements(ele_thick) THEN BEGIN
    print, 'Number of thicknesses does not equal number of elements.'
    print, 'Try again.'
    RETURN
  ENDIF

  n_0 = 1.                      ;index of medium at entrance of filter (assumed vacuum)
  n_t = 1.                      ;index of medium at exit of filter (assumed vacuum)

  ss = where(elem ne '', nx)
  n_files = XRT_CALIB+'/n_'+elem[0:nx-1]+'.txt'

  FOR x=0,(nx-1) DO BEGIN ;for each layer
                                ;read in index data from file

    readcol, n_files[x], dens, format='X,A', numline=1, /silent
    pos = strpos(dens, '=')
    dens = strmid(dens, pos[0]+1, 10)
    if (x eq 0)                        $
      then densities = dens            $
      else densities = [densities, dens]

    readcol, n_files[x], wavelength, delta, beta, skip=2, /silent

    wavelength = wavelength * 10.  ;; Convert nm to Angstroms.

;; find range in index data that is between wave_max and wave_min
    
    id = where(wavelength ge wave_min and wavelength le wave_max, count)
    wavelength = wavelength(id)
    delta = delta(id)
    beta = beta(id)

    IF x eq 0 THEN BEGIN    
                                ;remember count, wave for
      count_temp = count        ;first set of data
      wave_temp  = wavelength 

    ENDIF ELSE BEGIN

      delta = interpol(delta, wavelength, wave_temp) ;interpolate so ranges
      beta  = interpol(beta, wavelength, wave_temp)  ;are the same
       
      wavelength = wave_temp         ;use wave, count vars from
      count      = count_temp        ;first set

    ENDELSE

    index = dcomplex((1-delta), (1.*beta)) ;n = (1-delta) + i(beta)
    i     = dcomplex(0,1)                  ;i = sqrt(-1)

    sina = n_0*sin(ang)/index              ;Snell's law 
    cosa = sqrt(1-(sina)^2)

    k  = 2.*!pi*index*cosa/wavelength      ;define wavevector
    kl = k*ele_thick(x)                    ;multiply be thickness

    realkl = float(kl)                     ;find real part
    imagkl = imaginary(kl) < 300.          ;find imaginary part

    kl = dcomplex(realkl, imagkl)

;; define transfer matrix

    M = [[[cos(kl)], [-i*sin(kl)/index]], [[-i*index*sin(kl)], [cos(kl)]]]
    
    IF x eq 0 THEN MM=M ELSE BEGIN  ;remember first matrix

;;  multiply current M matrix with previous
           
      FOR j=0, count-1 DO  MM(j, *, *) =reform(MM(j,*,*))#reform(M(j,*,*))  
    ENDELSE
  ENDFOR

;; calculate transmission coefficient and transmittance

  t = (2.*n_0) /                                                    $
      (MM(*,0,0)*n_0 + MM(*,0,1)*n_0*n_t + MM(*,1,0) + MM(*,1,1)*n_t)

  trans     = abs(t)^2 * mesh_trans

  layers    = elem
  thickness = ele_thick

  FOR y=0, (n_elements(elem)-1) DO BEGIN

    elem_temp  = elem(y)
    thick_temp = strcompress(ele_thick(y), /rem)

    IF y eq 0 THEN BEGIN

      elem_title  = elem_temp
      thick_title = thick_temp

    ENDIF ELSE BEGIN

      elem_title  = elem_title + ':' + elem_temp
      thick_title = thick_title + ',' + thick_temp

    ENDELSE
  ENDFOR

  if keyword_set(file_name) then begin
    elem_title = elem_title + ' filter'

    fileout = file_name
    title   = 'Transmittance for ' + elem_title

    save, file=fileout,file_name,wavelength,trans,title,elem,ele_thick,$
        mesh_trans, densities

  endif

  RETURN 
END