; Interpolates the Mie scattering table for a given ssa and radius
; input -
;   ssa - scattering angle in degrees
;   rad - effective particle radius in nm. mean radius for gaussian, modal radius for lognormal
;   /lognorm - use lognormal distribution, shape parameter 1.4
;              Otherwise, use gaussian, 14 nm width
;   /max_r - return the maximum particle radius available in the table, instead of the phase function

function mie_phase,ssa,rad,lognorm=lognorm,max_r=max_r
  common mie_phase_common,mie_phase_f_gauss,mie_phase_f_lognorm,mie_phase_r,mie_phase_ssa
  common cips_mie_pfun,pfun_tag

  if keyword_set(lognorm) then begin
    if n_elements(mie_phase_f_lognorm) eq 0 then begin
      restore,!common_input_data+"/mie_lognorm_cips.sav"
      cips_pfun[0,*]=cips_pfun[1,*]
      mie_phase_f_lognorm=cips_pfun
    end
    if keyword_set(max_r) then return,(size(cips_pfun,/dim))[0]-1
    grid1=(n_elements(rad) ne n_elements(ssa))
    grid2=(n_elements(rad) ne 1)
    result=interpolate(mie_phase_f_lognorm,rad,ssa,grid=grid1)/interpolate(mie_phase_f_lognorm,rad,90,grid=grid2)
  end else begin
    if n_elements(rad) gt 1 and (n_elements(ssa) eq 1 or n_elements(ssa) eq n_elements(rad)) then begin
      if n_elements(ssa) gt 1 then begin
        result=ssa*!values.f_nan
        for i=0L,n_elements(result)-1 do result[i]=mie_phase(ssa[i],rad[i])
        return,result
      end else begin
        result=rad*!values.f_nan
        for i=0L,n_elements(result)-1 do result[i]=mie_phase(ssa,rad[i])
        return,result
      end 
    end
    if n_elements(mie_phase_f_gauss) eq 0 then begin
      if n_elements(pfun_tag) eq 0 then pfun_tag="cips_mie_pfun_14nm.sav"
      fn=!common_input_data+"/"+pfun_tag
      if ~file_test(fn) then begin
         print,"WARNING: "+fn+" was not found"
         fn=!common_input_data+"/cips_mie_pfun_14nm.sav"
         print,"WARNING: Using "+fn
      end
      print,fn
      restore,fn
      mie_phase_f_gauss=cips_pfun
      mie_phase_r=cips_rad
      mie_phase_ssa=cips_angle
    end
    if keyword_set(max_r) then return,(size(cips_pfun,/dim))[0]-1
;    grid1=(n_elements(rad) ne n_elements(ssa))
;    grid2=(n_elements(rad) ne 1)
;    result=interpolate(mie_phase_f_gauss,rad,ssa,grid=grid1)/interpolate(mie_phase_f_gauss,rad,90,grid=grid2)
     if ~finite(rad) then return,ssa*!values.f_nan
     i_rad=fix(rad)
     i_scat=fix(ssa)
   
     mie_s=size(mie_phase_f_gauss,/dimensions)
     if (i_rad ge mie_s[0]-1) || (max(i_scat) ge mie_s[1]-1) then begin
       result=ssa*!values.f_nan
     end else begin
       mp1 = mie_phase_f_gauss[i_rad  , i_scat]
       mp2 = mie_phase_f_gauss[i_rad+1, i_scat]
       w=where(mie_phase_ssa eq 90)
       w=w[0]
       mp90_1=mie_phase_f_gauss[i_rad,w]         ;Equivalent to normalization done in Scott's program (I hope)
       mp90_2=mie_phase_f_gauss[i_rad+1,w]
       mp1/=mp90_1
       mp2/=mp90_2
       result=mp2*0
       for iii=0L,n_elements(mp1)-1 do result[iii]=interpol([mp1[iii],mp2[iii]],[float(i_rad),float(i_rad+1)],rad)  
     end
    end
  if size(ssa,/n_dimensions) gt 0 then result=reform(result,size(ssa,/dimensions),/overwrite) else result=result[0]
  return,result
end