function ismtau,w,NH=NH,fH2=fH2,He1=He1,HeII=HeII,Fano=Fano,ikev=ikev,$
	wam=wam,bam=bam,mam=mam, _extra=e
;+
;function	ismtau
;	returns the optical depth for a given photon energy at specified
;	column density.  i.e., the $\tau$ of $e^{-\tau}$.
;
;	ingredients:
;	* in the EUV regime, cannibalizes ISMEUV.PRO for
;	  -- HI photoionization (Spitzer 19??, "Physical Processes in the
;	     Interstellar Medium", p105), between 0.03 keV (413.3A) and 911.75A
;	  -- HeII photoionization (ibid), between 0.03 keV (413.3A) and 911.75A
;	  -- HeI with auto-ionization resonances (Marr & West 1976;
;	     Oza 1986, Phys.Rev.A 33 824; Fernley et al. 1987, J.Phys.B 20,
;	     6457; Rumph, Bowyer, & Vennes 1994, AJ 107, 2108) below 503.97A
;	  -- NOTE: these are extended down to 43.657 AA when the X-ray
;	     calculations are not being done.
;	* in the X-ray regime (0.03 - 10 keV), either
;	  -- Morrison & McCammon (1983, ApJ 270, 119) polynomial fit between
;	     0.03 keV and 10 keV, or
;	  -- Balucinska-Church & McCammon (1992, ApJ 400, 699) cross-sections,
;	     which allows abundance variations, or
;	  -- Wilms, Allen, & McCray (2000, ApJ, in press) calculations
;	     that use cross-sections from Verner & Yakovlev (1995, A&AS,
;	     109, 125), which also allows abundance variations
;	* in the hard X-ray regime and beyond (>10 keV),
;	  -- Tanaka & Bleeker (1977, ???) power-law extended past 10 keV
;	* Molecular H (Cruddace et al. 1977, ApJ 187, 497; Kashyap et
;	  al. 1992, ApJ 391, 684) between 12.3985A (1 keV) to 911.75A
;
;	to extract only the cross-sections, set NH=1
;	to obtain transmission factors, then use exp(-ismtau(...))
;
;syntax
;	tau=ismtau(w,NH=NH,fH2=fH2,He1=He1,HeII=HeII,Fano=Fano,/ikev,$
;	/wam,/bam,/mam, abund=abund,noHeH=noHeH,verbose=verbose)
;
;parameters
;	w	[INPUT; required] photon energy values at which to compute
;		the optical depth.
;		* default units are [Ang], unless IKEV is set, when they are
;		  assumed to be [keV]
;
;keywords
;	NH	[INPUT; default=1] atomic H column density [cm^-2]
;	fH2	[INPUT; default=0.26] N(H2)/N(HI)
;		* the default value comes from the asymptotic value of
;		  N(H2)/N(HI) computed with the 3 component model (cold
;		  HI, warm HI, H2) of Bloemen, J.B.G.M. (1987, ApJ 322, 694)
;		* explicitly set fH2 to 0 to avoid including this factor.
;	He1	[INPUT; default=0.1*NH] neutral He column density [cm^-2]
;	HeII	[INPUT; default=0.1*He1] ionized He column density [cm^-2]
;	Fano	[INPUT] if set, the 4 strongest auto-ionizing resonances
;		of He I are included; the shape of these resonances are
;		given by a Fano profile (cf. Rumph, Bowyer, & Vennes 1994,
;		AJ 107, 2108).
;		* input wavelength grid between 190 A and 210 A should be
;		  ~0.01 A for this to be useful
;	ikev	[INPUT] if set, assumes that W are in units of keV
;	wam	[INPUT] if set, uses Wilms, Allen, & McCray cross-sections
;		* NOT IMPLEMENTED YET -- if set, defaults to Balucinska-Church
;		  & McCammon.
;	bam	[INPUT] if set, uses Balucinska-Church and McCammon updates
;		to Morrison & McCammon cross-sections in 30 eV - 10 keV range
;	mam	[INPUT] if set, uses Morrison and McCammon photoionization
;		cross-sections in the 30 eV - 10 keV range
;		* This is the default
;		* WAM takes precedence over BAM takes precedence over MAM.
;		* if MAM is explicitly set to 0 and neither BAM nor WAM are
;		  set, then the original ISMEUV code is used from 43 AA longwards
;	_extra	[INPUT ONLY] pass defined variables to subroutines
;		BAMABS: ABUND, VERBOSE, NOHEH
;
;history
;	written by Vinay Kashyap (Feb97), based on
;	  ISMEUV.PRO (W.Landsman [Oct94], via ISM.C@cea-ftp.cea.berkeley.edu
;	  (Pat Jelinsky, Todd Rumph, et al.)) and ISMABS.PRO (VK; Dec91)
;	changed default on HeII from 0 to 0.1*He1 (VK; Jan MM)
;	changed default on NH from 1e20 to 1.0 (VK; MarMM)
;	changed keyword KEV to IKEV; also for some reason H and He
;	  photoionization from ISMEUV had been cut off at the C edge,
;	  now corrected; added keywords WAM, BAM, and MAM, and call to BAMABS
;	  (VK; OctMM)
;	bug correction: BAM was not filtering on energy; allowed setting
;	  MAM=0 to recover original ISMEUV calc (VK; Jun03)
;	follow-up bug correction (VK; Jul03)
;-

;	usage
nw=n_elements(w)
if nw eq 0 then begin
  print,'Usage: tau=ismtau(W,NH=NH,fH2=fH2,He1=He1,HeII=HeII,/Fano,/ikev,$'
  print,'       /wam,/bam,/mam, abund=abund,verbose=verbose)'
  print,'  returns optical depth at given photon energy for given column'
  return,0.
endif

nrg=[float(abs(w))]	;avoid touching the input, convert scalar to vector,
			;avoid -ve values in wvl/keV, ensure R*4

;	check keywords
if n_elements(NH) eq 0 then NH=1.			;N(HI)
if n_elements(fH2) eq 0 then fH2=0.26			;N(H2)/N(HI)
if n_elements(He1) eq 0 then He1=0.1*NH			;N(HeI)
if n_elements(HeII) eq 0 then HeII=0.1*He1		;N(HeII)
;	if input photon "energy" is in [Ang], convert to [keV]
if not keyword_set(ikeV) then begin
  wvl=nrg & oo=where(nrg gt 0,moo)
  if moo gt 0 then nrg(oo)=12.3985/nrg(oo) else return,0*w
endif else begin
  wvl=0*nrg & oo=where(nrg gt 0,moo)
  if moo gt 0 then wvl(oo)=12.3985/nrg(oo) else return,0*w
endelse
;
iMAM=1 & iBAM=0 & iWAM=0
if n_elements(MAM) gt 0 then begin
  if MAM(0) eq 0 then begin
    iMAM=0 & iBAM=0 & iWAM=0
  endif
endif
if keyword_set(BAM) then begin
  iMAM=0 & iBAM=1 & iWAM=0
endif
if keyword_set(WAM) then begin
  message,'Wilms, Allen, McCray cross-sections not implemented yet',/info
  message,'using Balucinska-Church & McCammon data',/info
  iMAM=0 & iBAM=1 & iWAM=0
endif

;	initialize
tauH=0*wvl & tauH2=tauH & tauHeII=tauH & tauHeI=tauH & tauX=tauH & tauTB=tauH

;	HI photoionization (from ISMEUV)
ev30 = 12.3985/(30.e-3)
	;Ang corresponding to 30 eV, boundary of ISMEUV applicability
	;when WAM, BAM, or MAM is set
if iMAM eq 0 and iBAM eq 0 and iWAM eq 0 then ev30=43.657
r=wvl/911.75
;DO NOT USE THIS: oo=where(r ge 43.657/911.75 and r lt 1,moo)
oo=where(r ge ev30/911.75 and r lt 1,moo)
if moo gt 0 then begin
  z=0*r & z(oo)=sqrt(r(oo)/(1.-r(oo)))
  tauH(oo) = NH(0) * 3.44e-16 * (r(oo)^4) *$
    exp(-4.0*z(oo)*atan(1./z(oo)))/(1.0-exp(-2*!pi*z(oo)))
endif

;	HeII photoionization (from ISMEUV); just like above, but charge=2
r=4*wvl/911.75
;WRONG-> oo=where(r gt 43.657/911.75 and r lt 1,moo)
oo=where(r ge ev30/911.75 and r lt 1,moo)
if moo gt 0 then begin
  z=0*r & z(oo)=sqrt(r(oo)/(1.-r(oo)))
  tauHeII(oo)=HeII(0)*3.44e-16*(r(oo)^4)*$
    exp(-4.0*z(oo)*atan(1./z(oo)))/((1.0-exp(-2*!pi*z(oo)))*4)
endif

;	HeI (from ISMEUV)
;	for wvl>46A (Marr & West 1976)
c1=[ -2.953607d+01, 7.083061d+00,  8.678646d-01, -1.221932d+00,$
      4.052997d-02, 1.317109d-01, -3.265795d-02,  2.500933d-03 ]
;	for wvl<46A (Marr & West 1976)
c2=[ -2.465188d+01, 4.354679d+00, -3.553024d+00, 5.573040d+00,$
     -5.872938d+00, 3.720797d+00, -1.226919d+00, 1.576657d-01 ]
;	parameters for auto-ionization resonances
;	from Oza (1986, Phys.Rev.A 33, 824)
q=[ 2.81d, 2.51d, 2.45d, 2.44d ]
;	from Fernley et al. (1987, J.Phys.B 20,6457)
nu=[ 1.610d, 2.795d, 3.817d, 4.824d ]
fano_gamma=[ 2.64061d-03, 6.20116d-04, 2.56061d-04, 1.320159d-04 ]
esubi=3.0d - 1.0d/nu^2 + 1.807317d
;
;INCORRECT-> oo=where(wvl gt 0 and wvl lt 503.97,moo)
oo=where(wvl gt ev30 and wvl lt 503.97,moo) & x=0*wvl & y=x & eps=x
if moo gt 0 then begin
  x(oo)=alog10(wvl(oo))
  o1=where(wvl ge 46.0 and wvl lt 503.97,mo1)
  o2=where(wvl gt 0 and wvl lt 46.0,mo2)
  if mo1 gt 0 then y(o1)=poly(x(o1),c1)
  if mo2 gt 0 then y(o2)=poly(x(o2),c2)
  if keyword_set(Fano) then begin
    eps(oo)=911.2671/wvl(oo)
    for i=0,3 do begin
      x(oo)=2.0*((eps(oo)-esubi(i))/fano_gamma(i))
      y(oo)=y(oo)+alog10((x(oo)-q(i))^2/(1.+x(oo)*x(oo)))
    endfor
  endif
  tauHeI(oo)=He1(0) * 10.^(y(oo))
endif

;	H2 (from ISMABS)
if not keyword_set(iwam) then begin
  ; Wilms, Allen & McCray have their own computation of H2 cross-section
  oo=where(nrg gt 0 and nrg lt 1,moo)
  if moo gt 0 then tauH2(oo)=NH(0)*fh2(0)*0.16e-22*nrg(oo)^(-3.54)
endif

;	X-Ray regime (from ISMABS)

if keyword_set(iMAM) then begin
  ;	Morrison & McCammon (1983, ApJ 270, 119) fit coefficients
  c0 = [ 17.3,      34.6,  78.1,  71.4,  95.5,  308.9, 120.6, 141.3, 202.7,$
       342.7,    352.2, 433.9, 629.0, 701.2 ]
  c1 = [ 608.1,    267.9,  18.8,  66.8, 145.8, -380.6, 169.3, 146.8, 104.7,$
       18.7,      18.7,  -2.4,  30.9,  25.2 ]
  c2 = [ -2150.0, -476.1,   4.3, -51.4, -61.1,  294.0, -47.7, -31.5, -17.0,$
       0.0,        0.0,  0.75,   0.0,   0.0 ]
  el = [ 0.030,    0.100, 0.284, 0.400, 0.532,  0.707, 0.867, 1.303, 1.840,$
       2.471,    3.210, 4.038, 7.111, 8.331 ]
  eh = [ 0.100,    0.284, 0.400, 0.532, 0.707,  0.867, 1.303, 1.840, 2.471,$
       3.210,    4.038, 7.111, 8.331, 10.00 ]
  nbin=n_elements(el) & sig0=1e-24
  for i=0,nbin-1 do begin
    oo=where(nrg gt el(i) and nrg le eh(i),moo)
    if moo gt 0 then tauX(oo)=NH(0)*sig0*$
      (c0(i)+c1(i)*nrg(oo)+c2(i)*nrg(oo)^2)/(nrg(oo))^3
  endfor
endif

if keyword_set(iBAM) then begin
  ;	Balucinska-Church & McCammon, 1992, ApJ 400, 699
  oo=where(nrg ge 0.03 and nrg le 10.,moo)
  if moo gt 0 then tauX(oo)=bamabs(nrg(oo),/ikeV,Fano=Fano, _extra=e)*NH(0)
endif

;	Hard X-ray regime (from ISMABS)
;	Tanaka & Bleeker coefficients
c0=2e-22 & c1=-2.5 & oo=where(nrg gt 10,moo)
if moo gt 0 then tauTB(oo)=NH(0)*c0*(nrg(oo))^(c1)

;	add up ever'thin
tau = tauH + tauH2 + tauHeII + tauHeI + tauX + tauTB

return,tau
end