c *****************************************
c Compton reflection from relativistic disk 
      subroutine relrefl(xr,refl,xm,gamma,sprf,  
     + par_refl,nrefl,jref,ioniz) 
c *****************************************
      implicit none
      integer ifl
      integer nrefl,jref,ioniz 
      real xr(jref),refl(jref) 
      real par_refl(nrefl)
      real gamma,ab_,ab0,temp,xi,xm,ec
c parameters number for the model used
      integer ndisk,nem
      PARAMETER(nem=5000)

      parameter(ndisk=6)
      integer nx
      real xmin,xmax
c internal integration parameters
      parameter(xmin=2.e-5,xmax=4.,nx=600)
      integer n_disk
      real disk_xmin, disk_xmax
      parameter(disk_xmin = 0.7, disk_xmax = 1.5, n_disk = 400)
     
      logical firstflag,fl1
      integer i,j,k,jlo,jhi
      real x(nx),x_disk(0:n_disk),dx
      real phot_ref(nx),phot_disk(n_disk),phot_rs(nx),spref(nx) 
      real logxmin, logxmax, logdmin, logdmax
      real pp,r,xlo,xhi,disx
      real in_disk(ndisk)
      real spint,eps,dln  
      external sprf 
      data firstflag/.true./,fl1/.true./ 
      data eps/1e-5/,ec/0./,ifl/0/   
      data in_disk/ndisk*9999./
      save firstflag,x,phot_ref,phot_disk,in_disk,
     >     phot_rs,x_disk,spref,dx

c     Driver for angle-dependent reflection from an exponentially-cutoff
c     power law and ionized medium.
c     It needs to be compiled together with xspexrav.f (neutral reflection).
c
c     See Magdziarz & Zdziarski, 1995, MNRAS.
c     See Zdziarski et al., 1995, ApJL January for description of 
c     calculation of ionization (based on Done et al. 1992).
c     The abundances are of Morrison & McCammon.
c
c     The output spectrum is the sum of the cutoff power law and the 
c     reflection component. 
c     The reflection component alone can be obtained 
c     for scale (see below) = rel_refl < 0. Then the actual 
c     reflection normalization is |scale|. Note that you need to 
c     change then the limits of rel_refl. The range of rel_refl in that case 
c     should exclude zero (as then the direct component appears).
c     
c     If E_c=0, there is no cutoff in the power law. 
c
c     This version allows for changes of the vector 'ear' between subsequent
c       calls.
c     and then calculates the broanening in a schwartzchild metric
c
c Fe abudance
      ab0=par_refl(1)
c  ab_       - abundance of elements heavier than He relative to
c              the solar abundances
      ab_=par_refl(2)
c  xi        - disk ionization parameter = L/nR^2
      xi=par_refl(3)
c  temp      - disk temperature (for ionization balance, Tmax=10^6 K)
      temp=par_refl(4)

      dln=alog(10.)

      if(firstflag)then
      logxmin = log10(xmin)
      logxmax = log10(xmax) 
      logdmin = log10(disk_xmin)
      logdmax = log10(disk_xmax) 
         dx=log10(xmax/xmin)/(nx-1)
         do i=1,nx
            x(i)=xmin*10**((i-1)*dx)
         enddo
         disx = log10(disk_xmax / disk_xmin) / n_disk
         x_disk(0) = disk_xmin
         do i = 1, n_disk
           x_disk(i) = disk_xmin * 10**(i * disx)
         enddo
      firstflag=.false. 
      dx=dx*dln
      endif

      if(ioniz.eq.1) 
     &  call greeion(x,spref,xm,ab_,ab0,sprf,gamma,ec,temp,xi,nx)
      if(ioniz.eq.0) 
     &  call greeneu(x,spref,xm,ab_,ab0,sprf,gamma,ec,nx)
c spref= x Fx
         do i=1,nx
         phot_ref(i) = spref(i)/x(i)**2  
         enddo 
       

c beta 
      if(in_disk(2).ne.par_refl(5))then
         in_disk(2)=par_refl(5)
         fl1=.true.
      endif
c Rin 
      if(in_disk(3).ne.par_refl(6))then
         in_disk(3)=par_refl(6)
         fl1=.true.
      endif
c Rout
      if(in_disk(4).ne.par_refl(7))then
         in_disk(4)=par_refl(7)
         fl1=.true.
      endif
c inclination 
      pp=180./3.1415*acos(xm) 
      if(in_disk(5).ne.pp)then
         in_disk(5)=pp 
         fl1=.true.
      endif
c redshift (we are in the source frame) 
      if(in_disk(6).ne.0.)then
         in_disk(6)=0. 
         fl1=.true.
      endif
      
      if(fl1) then
        in_disk(1) = 1.0
        call relline(x_disk, n_disk, in_disk, ifl, phot_disk)
         do i=1, n_disk
           if(phot_disk(i) .ge. 1.0e-8) then
             phot_disk(i) = phot_disk(i) / (x_disk(i) - x_disk(i-1))
           else
             phot_disk(i) = 0.0
           endif
         enddo
      endif

     
        do i=1,nx
           phot_rs(i)=1e-20 
        enddo
        do i = 1, nx
c f'(x)=int dx' f(x') line(x/x')
c x(i)=x 
          xlo = disk_xmin * x(i)
          xhi = disk_xmax * x(i)
          jlo = int(nx * (log10(xlo) - logxmin) / (logxmax-logxmin))
          jhi = int(nx * (log10(xhi) - logxmin) / (logxmax-logxmin))
          jlo = max(jlo, 1)
          jhi = min(jhi, nx)
          do j = jlo, jhi
c x(j)=x'
            r = x(i) / x(j)
            k = int(n_disk * (log10(r) - logdmin) / (logdmax-logdmin))
            if((k.ge.1) .and. (k.le.n_disk)) then
            phot_rs(i) = phot_rs(i) + phot_ref(j) * phot_disk(k)
            endif
          enddo
        enddo
      do i=1,nx
      phot_rs(i)=phot_rs(i)*dx 
      enddo

      fl1=.false.

      do i=1,jref
         pp=spint(xr(i),x,nx,phot_rs)
         refl(i)=pp*xr(i)**2 
      enddo

      return
      end

c ------------------------------------------------------------------- c
c interpolation 
      real function spint(xx,xt,jt,spt)
 
      implicit   none
      integer    jt
      real       xx,xt(jt),spt(jt)
 
      integer    il,ih
      real       sp
      save       ih
      data       ih/2/

      if((xx.gt.xt(1)).and.(xx.lt.xt(jt)))then
         if(xx.lt.xt(ih))then
            ih=2
         endif
 10      if(xx.gt.xt(ih))then
            ih=ih+1
            goto 10
         else
            il=ih-1
c            sp=spt(il)+
c     &         (spt(ih)-spt(il))*
c     &         (xx-xt(il))/(xt(ih)-xt(il))
        sp=spt(il)*exp(log(spt(ih)/spt(il))*
     *     log(xx/xt(il))/log(xt(ih)/xt(il))) 
         endif
      else
         sp=1e-20 
      endif
 
      spint=sp
 
      return
      end


      SUBROUTINE relline(earp,ne,param,ifl,photar)

      implicit none
      INTEGER ne,ifl
      REAL earp(0:ne),param(*),photar(ne)
      real ear(0:10000)

c  model to calculate the line shape for a rotating accretion
c  disk. does not include GR effects. note that if param(2) is
c  set to 10 then do the special case of an accretion disk
c  emissivity.
c  modified 11/30/95 by P.Zycki to include object's redshift
c  corrected 20/02/96 by P.Magdziarz 

c  parameters :
c	1        line energy
c       2        power law index for emissivity (10 for disk)
c       3        inner radius (GM/c**2)
c       4        outer radius (GM/c**2)
c       5        inclination  (degrees)
c       6        redshift

      INTEGER nr
      PARAMETER (nr=100)

      INTEGER n, k, ie
      REAL ri, ro, dlogr, dera, alp,
     &     tgcsi2, cosdel, fra, ra, dra,
     &     rd(nr), enel, spm, rilog10, zfac  
      REAL beta, betal, betah, zpo, zpol, zpoh
      REAL sincl, cincl, rafact, flux, total
      REAL enobl, enobh

      DATA dera /57.295779/

c convert input parameters - inclination set to radians.

      enel = param(1)
      alp = param(2)
      ri = param(3)
      rilog10 = alog10(ri)
      ro = param(4)
      zfac = 1+param(6)

      sincl = sin(param(5)/dera)
      cincl = cos(param(5)/dera)

      do n=1,ne
         ear(n) = zfac*earp(n)
      end do

c set spectrum to 0

      DO n = 1, ne
         photar(n) = 0.
      ENDDO

c trap case where inner radius is greater than outer

cjp
c      IF ( ri .GE. ro ) THEN
c         CALL xerror(
c     &     'Inner radius > outer radius  -  model is invalid', 5)
c         RETURN
c      ENDIF

c trap case of inner radius being less than the last stable orbit.

c      IF  (ri .LT. 6.) THEN
c         CALL xerror('Inner radius < 6 R_g  -  model is invalid', 5)
c         RETURN
c      ENDIF

c calculate log step in radius

      dlogr = (alog10(ro)-rilog10)/float(nr-1)

c calculate radii

      rd(1) = ri
      DO n = 2, nr
         rd(n) = 10.**(rilog10+float(n-1)*dlogr)
      ENDDO

c big loop for radii. note that the radius cannot reach 3 else
c the metric goes singular

      DO n = 1, nr - 1

         ra = (rd(n)+rd(n+1))/2.

         dra = rd(n+1) - rd(n)

         rafact = sqrt(1.-3./ra)

c if power-law index is less than ten use to calculate emissivity

         IF (alp.LT.9.9) THEN

            fra = ra**(alp)

c else use the accretion disk emissivity law.

         ELSE

            fra = (1.-sqrt(6./ra))/ra**3

         ENDIF

c loop over azimuthal angles

         DO k = 1, 179, 2

            beta = float(k)/dera
            betal = float(k-1)/dera
            betah = float(k+1)/dera

c calculate mean redshift (1+z = zpo) for the bin

            tgcsi2 = (sincl*sin(beta))
     &               **2/(1.-(sincl*sin(beta))**2)
            cosdel = sincl*cos(beta)
     &               /sqrt((sincl*cos(beta))**2+cincl**2)
            zpo = (1.+cosdel/sqrt(ra*(1.+tgcsi2)-2.))/rafact

c and the low and high redshifts for the bin

            tgcsi2 = (sincl*sin(betal))
     &               **2/(1.-(sincl*sin(betal))**2)
            cosdel = sincl*cos(betal)
     &               /sqrt((sincl*cos(betal))**2+cincl**2)
            zpol = (1.+cosdel/sqrt(ra*(1.+tgcsi2)-2.))/rafact

            tgcsi2 = (sincl*sin(betah))
     &               **2/(1.-(sincl*sin(betah))**2)
            cosdel = sincl*cos(betah)
     &               /sqrt((sincl*cos(betah))**2+cincl**2)
            zpoh = (1.+cosdel/sqrt(ra*(1.+tgcsi2)-2.))/rafact


c  enobl and enobh are the lower and upper observed energy from
c  this azimuthal and radial bin

            enobl = min(enel/zpol, enel/zpoh)
            enobh = max(enel/zpol, enel/zpoh)
            total = enobh - enobl

c  calculate total emission from this bin

            flux = fra*ra*dra*2./dera/zpo**3

c  find fractions of emission from this bin to place in each energy
c  range.

            IF (enobh .LT. ear(0) .OR. enobl .GT. ear(ne) ) GOTO 10

            DO ie = 1, ne
               IF ( ear(ie) .GE. enobh ) THEN
                  IF ( ear(ie-1) .LE. enobl ) THEN
                     photar(ie) = photar(ie) + flux
                     GOTO 10
                  ELSEIF ( ear(ie-1) .GE. enobh ) THEN
                     GOTO 10
                  ELSE
                     IF ( total .GT. 0. ) THEN
                        photar(ie) = photar(ie) 
     &                     + flux*(enobh-ear(ie-1))/total
                     ENDIF
                  ENDIF
               ELSEIF ( ear(ie) .GE. enobl ) THEN
                  IF ( ear(ie-1) .GE. enobl ) THEN
                     IF ( total .GT. 0. ) THEN
                        photar(ie) = photar(ie) 
     &                     + flux*(ear(ie)-ear(ie-1))/total
                     ENDIF
                  ELSE
                     IF ( total .GT. 0. ) THEN
                        photar(ie) = photar(ie) 
     &                     + flux*(ear(ie)-enobl)/total
                     ENDIF
                  ENDIF
               ENDIF
            ENDDO

 10         CONTINUE

         ENDDO
      ENDDO

c normalise values to total

      spm = 0.
      DO n = 1, ne
         spm = spm + photar(n)
      ENDDO

c write values

      IF ( spm .NE. 0 ) THEN
         DO n = 1, ne
            photar(n) = photar(n)/spm
         ENDDO
      ENDIF

      RETURN
      END

c**************************************************
C REFLECTION ROUTINES  FOR NEUTRAL AND IONIZED REFLECTION 
c**************************************************
      subroutine greeneu(x,spref,xm,ab_,ab0,spp,gamma,ec,jmax)
c     calculates the reflected spectrum
c     ab0 - iron log10 abundance relative to hydrogen
      implicit none 
      integer nem 
      PARAMETER(nem=5000) 
      REAL ab0,ab_,ddy,dy,dym,Ec,fjc,Gamma,GNR1,GRXYH1,
     &     GRXYL1,SIGMANEU,SPP,sr,srp,xjd,xjl,xkabs,Xm, 
     &     xmin,xmo
      REAL xnor , xrefmax , xtransh , xtransl , y , y0 , ymin
      INTEGER i,j,Jmax,jmin,jrefmax,jtransh,jtransl,nc,ncp 
      real x(nem),spref(nem)
      real pm1(19),pmx(26),ap(3)
      external spp
 
c  precision factor for Green's function integration
      ncp=100
c  ranges for different methods
      xmin=2.e-4
      xtransl=0.01957
      xtransh=0.02348
      ymin=.03333
c  angle dependent parameters
      xmo=max(.05,min(.95,xm))
      call pm1yn(xmo,pm1)
      call pmxyn(xmo,pmx)
      ap(2)=0.802-1.019*xmo+2.528*xmo**2-3.198*xmo**3+
     &      1.457*xmo**4+(.030/xmo)**4
      ap(3)=0.381
      xrefmax=1/(ymin+pmx(26))
 
      jmin=0
      do 50 j=1,jmax
         if(x(j).gt.xmin) goto 51
         jmin=j
 50   continue
 51   jtransl=0
      do 52 j=max(1,jmin),jmax
         if(x(j).gt.xtransl) goto 53
         jtransl=j
 52   continue
 53   jtransh=0
      do 54 j=max(1,jtransl),jmax
         if(x(j).gt.xtransh) goto 55
         jtransh=j
 54   continue
 55   jrefmax=0
      do 56 j=max(1,jtransh),jmax
         if(x(j).gt.xrefmax) goto 57
         jrefmax=j
 56   continue
 57   continue
      do 10 j=1,jmin
         spref(j)=0
 10   continue
c-----------------------
      do 20 j=jmin+1,jtransl
         xkabs=sigmaneu(j,x,jtransh,ab_,ab0,xnor) 
         spref(j)=spp(1/x(j),gamma,ec)*gnr1(xkabs,xmo)
 20   continue
c-----------------------
      if(jmin.ge.jtransl) then
         xkabs=sigmaneu(1,x,max(1,jtransh),ab_,ab0,xnor)
      endif
      ap(1)=xnor/1.21
      xjl=alog(xtransl)
      xjd=alog(xtransh)-xjl
      do 60 j=jtransl+1,jtransh
         fjc=.5*sin(3.14159*((alog(x(j))-xjl)/xjd-.5))
         xkabs=sigmaneu(j,x,jtransh,ab_,ab0,xnor)
         y=1/x(j)
         spref(j)=spp(y,gamma,ec)*gnr1(xkabs,xmo)
         dym=y-ymin
         dy=min(2,dym)
         ddy=(dy-pmx(26))/(ncp+1)
         y0=y-dy
         sr=spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
         dy=pmx(26)
         y0=y-dy
         sr=.5*(sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap))
         do 62 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
 62      continue
         sr=sr*ddy
         if(dym.gt.2.)then
            ddy=dym-2
            nc=int(ncp*ddy/(dym-pmx(26)))
            ddy=ddy/(nc+1)
            dy=dym
            y0=y-dy
            srp=spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
            dy=2
            y0=y-dy
            srp=.5*(srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap))
            do 64 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
 64         continue
            sr=sr+srp*ddy
         endif
         spref(j)=(.5-fjc)*spref(j)+(.5+fjc)*sr*x(j)
 60   continue
c-----------------------
      do 30 j=jtransh+1,jrefmax
         y=1/x(j)
         dym=y-ymin
         dy=min(2,dym)
         ddy=(dy-pmx(26))/(ncp+1)
         y0=y-dy
         sr=spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
         dy=pmx(26)
         y0=y-dy
         sr=.5*(sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap))
         do 32 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
 32      continue
         sr=sr*ddy
         if(dym.gt.2.)then
            ddy=dym-2
            nc=int(ncp*ddy/(dym-pmx(26)))
            ddy=ddy/(nc+1)
            dy=dym
            y0=y-dy
            srp=spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
            dy=2
            y0=y-dy
            srp=.5*(srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap))
            do 34 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
 34         continue
            sr=sr+srp*ddy
         endif
         spref(j)=sr*x(j)
 30   continue
c----------------------
      do 40 j=jrefmax+1,jmax
        spref(j)=0
 40   continue
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Reflection factor for nonrelativistic Compton reflection,
c  reflection orders >1 approximated with isotropic phase function,
c  and approximation of H-function by Basko ApJ 223,268,(1978),
c  valid up to 14keV;
c  sig - sigma_bf per hydrogen atom in sigma_Thomson units
c  xm  - cosine of the angle between normal to the slab and the
c        observation direction
 
      real function    gnr1(sig,xm)
      implicit         none
      real             sig,xm
      real             al0,al1,al2,xl,x2,x4,sq3
      parameter        (sq3=1.732051)
 
      al0=1/(1+sig/1.21)
      al1=sqrt(1-al0)
      al2=sq3*al1
      xl=log(1+1/xm)
      x2=xm*xm
      x4=x2*x2
 
      gnr1=.5*al0*xm*(.375*((3-2*x2+3*x4)*xl+(3*x2-1)*(.5-xm))+
     & ((1+(1/al1-1)*(1-log(1+al2)/al2))*(1+xm*sq3)/(1+xm*al2)-1)*xl)
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Green's functions for Compton reflection with asymtotic
c  absorption included; output is (dy+y0)*G(mu,dy,y0)*W(mu,dy,y0),
c  for .0316<y0<100.; in two parts according to dy=2.
c  grxyh1 - subroutine for dy>2
c  grxyl1 - subroutine for pmx(26)<dy<2
c  pm1   - vector of angle-dependent coefficients,
c          set-up by pm1yn subroutine
c  pmx   - vector of angle-dependent coefficients,
c          set-up by pmxyn subroutine
c  dy    - wavelength shift dy=1/x-1/x0
c  y0    - incident photon wavelength y0=1/x0
c  ap    - absorption parameters
 
      real function    grxyh1(pm1,pmx,dy,y0,ap)
      implicit         none
      real             pm1(19),pmx(26),dy,y0,ap(3)
      real             pyx(10),p1,p4,gr,aic1
 
      p1=log(y0)
      p4=1/y0
      pyx(5)=.230+.033*p1+.763*y0**(-.507)+.007*y0**(-1.553)
      pyx(6)=75.180*p4
      pyx(7)=-.100+.115*p1+.100*y0**1.623
      if(y0.lt.1.)then
         pyx(8)=pmx(19)+pmx(20)*p1
      else
         pyx(8)=pmx(21)+pmx(22)*p1+pmx(23)*(pmx(24)+p4)**pmx(25)
      endif
      p1=1/(y0+dy)
      gr=pyx(5)*(1+pyx(6)*p1)**(pyx(7)*(1+pyx(8)*p1))
      p1=1+dy
      gr=gr*pm1(18)*dy**(-1.5)*p1*
     &   (1+pm1(14)/p1)**(pm1(15)+pm1(16)*p1**pm1(17))
      gr=gr*exp(aic1(y0,(dy-.15),ap(1)))
 
      grxyh1=gr
 
      return
      end
 
      real function    grxyl1(pm1,pmx,dy,y0,ap)
      implicit         none
      real             pm1(19),pmx(26),dy,y0,ap(3)
      real             pyx(10),p1,p2,p3,p4,gr,aic1,xla0,xla
 
      p1=log(min(y0,10.))
      p2=p1*p1
      p3=p2*p1
      pyx(1)=p1*(pmx( 1)+pmx( 2)*p1+pmx( 3)*p2+pmx( 4)*p3)
      pyx(2)=p1*(pmx( 5)+pmx( 6)*p1+pmx( 7)*p2+pmx( 8)*p3)
      pyx(3)=p1*(pmx( 9)+pmx(10)*p1+pmx(11)*p2+pmx(12)*p3)
      pyx(4)=p1*(pmx(13)+pmx(14)*p1+pmx(15)*p2+pmx(16)*p3
     &                                             +pmx(17)*p3*p1)
      pyx(9)=1+pmx(18)*(.326*((y0*.1)**1.692-1)+p1-2.303)
      if(y0.lt.10.)then
         p1=y0+2
         p1=log(p1/(y0+dy))/log(p1/y0)
         p2=p1*p1
         p3=p2*p1
         gr=exp(pyx(1)+pyx(2)*p1+pyx(3)*p2+pyx(4)*p3)
      else
         p4=1/(10+dy)
         p1=log(12*p4)*5.4848
         p2=p1*p1
         p3=p2*p1
         p4=pyx(9)*(y0+dy)*p4
         gr=p4*exp(pyx(1)+pyx(2)*p1+pyx(3)*p2+pyx(4)*p3)
      endif
      p1=pm1(1)+pm1(2)*dy
      p2=(pm1(3)+pm1(4)*max(0.,(pm1(5)-dy))**pm1(6)+
     &   pm1(7)*(pm1(8)-dy)**pm1(9))**pm1(10)/
     &   (1+exp(pm1(11)*(pm1(12)-dy**pm1(13))))
      gr=gr*min(p1,p2)
      xla0=1/(1+ap(1)*y0**3)
      xla=1/(1+ap(1)*(y0+dy)**3)
      gr=gr*xla0*(ap(2)+(1-ap(2))*(1+(exp(ap(3)*dy)-2)*sqrt(1-xla))
     &     *exp(aic1(y0,dy,ap(1))))
 
      grxyl1=gr
 
      return
      end
 
      real function    aic1(y0,dy,apc)
      implicit         none
      real             y0,dy,apc
      real             yp,a0,a1,c0,c1,c2,d0l,d0h,d1l,d1h
      real             aih,ail
 
      yp=y0+dy
      a0=apc
      a1=apc**.3333
      c0=.5/a1
      c1=-3.4641*c0
      c2= 0.5774
      d0l=a0*y0**3
      d0h=a0*yp**3
      d1l=a1*y0
      d1h=a1*yp
      ail=y0*(3-log(1+d0l))+c0*log((1-d1l+d1l**2)/(1+d1l)**2)
     &                     +c1*atan(c2*(2*d1l-1))
      aih=yp*(3-log(1+d0h))+c0*log((1-d1h+d1h**2)/(1+d1h)**2)
     &                     +c1*atan(c2*(2*d1h-1))
 
      aic1=aih-ail
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Angle dependent coefficients for Green's function for Compton
c  reflection;
c  xm  - cosine of the angle between normal to the slab and the
c        observation direction
c  pmx - output coefficients vector
 
      subroutine       pmxyn(xm,pmx)
      implicit         none
      real             xm,pmx(26)
      real             p2,p3,p4
 
              p2=xm*xm
              p3=p2*xm
              p4=p3*xm
      pmx( 1)= .8827-.5885*xm+.7988*p2-.4117*p3
      pmx( 2)= .0517+.1076*xm-.1691*p2+.0678*p3
      pmx( 3)= .0014+.0043*xm
      pmx( 4)=-.0003-.0027*xm+.0021*p2
      pmx( 5)=-1.3259+1.6214*xm-3.7007*p2+2.1722*p3
     &        +exp(20.3*(xm-1.063))
      pmx( 6)=  .0790- .4029*xm+ .6619*p2+ .2210*p3
     &        -exp(17.7*(xm-1.041))
      pmx( 7)=  .0751- .1848*xm+ .4068*p2- .4126*p3
     &        +exp( 9.8*(xm-1.185))
      pmx( 8)=- .0020- .0394*xm+ .1004*p2- .0597*p3
      pmx( 9)= 3.2953-3.6996*xm+7.9837*p2-4.6855*p3
      pmx(10)=- .5278+ .9857*xm-1.7454*p2- .3464*p3
     &        +exp(27.6*(xm-0.959))
      pmx(11)=- .1919+ .5798*xm-1.2879*p2+1.4885*p3
     &        -exp(18.3*(xm-0.986))
      pmx(12)=  .0200+ .0832*xm- .0333*p2- .2370*p3
     &        +exp(16.6*(xm-1.086))
      pmx(13)=-2.2779+2.4031*xm-4.0733*p2+1.9499*p3
     &        -exp(19.6*(xm-0.968))
      pmx(14)=  .4790-1.0166*xm+3.1727*p2-4.0108*p3+3.0545*p4
     &        -exp(30.4*(xm-0.957))
      pmx(15)=  .1122- .3580*xm+ .4985*p2- .3750*p3- .5349*p4
     &        +exp(31.2*(xm-0.972))
      pmx(16)=- .0160- .0471*xm- .0536*p2+ .2006*p3+ .2929*p4
     &        -exp(30.6*(xm-1.009))
      pmx(17)=  .0005+ .0021*xm- .0447*p2+ .1749*p3- .2303*p4
     &        +exp(16.9*(xm-1.130))
      pmx(18)=.006+.089*xm-.102*p2+.056*p3
      pmx(19)=.618-.830*xm
      pmx(20)=.128-.132*xm
      pmx(21)=.632-.875*xm
      pmx(22)=-.672-.286*xm+.717*p2-.481*p3
      pmx(23)=.0126-.0160*xm+.0077*p2
      pmx(24)=.0111+.0030*xm-.0014*p2
      pmx(25)=-2.437-.328*xm-.260*p2+.279*p3
      pmx(26)=1-sqrt(1-(xm-.05)**2)
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Angle dependent coefficients for Green's function for Compton
c  reflection of incident photon at energy x0=1;
c  xm  - cosine of the angle between normal to the slab and the
c        observation direction
c  pm1 - output coefficients vector
 
      subroutine       pm1yn(xm,pm1)
      implicit         none
      real             xm,pm1(19)
      real             p2,p3,p4
 
              p2=.911-.549*(1-xm)**1.471
              p3=.254-.041*xm**(-3.798)
      pm1( 1)=p2-2*p3
      pm1( 2)=p3
      pm1( 3)=.161+.439*xm**.189-.791*xm**7.789
      pm1( 4)=max(0.,(-.871+1.740*xm**.176-1.088*xm**.907))
      pm1( 5)=.934+.054*xm**(-.666)
              p4=min(.524,xm)
      pm1( 6)=4.647-11.574*p4+11.046*p4*p4
      pm1( 7)=.012+.199*xm**.939+.861*xm**8.999
      pm1( 8)=2.150
              p2=xm*xm
              p3=p2*xm
      pm1( 9)=-1.281+2.374*xm-4.332*p2+2.630*p3
      pm1(10)=2.882+.035*xm**(-1.065)
      pm1(11)=30.28+69.29*xm
      pm1(12)=1.037*(1-sqrt(1-.874*p2))**.123
      pm1(13)=.123
              p4=p3*xm
      pm1(14)=56.50-27.54*xm+671.8*p2-1245.*p3+708.4*p4
      pm1(15)=-.897-.826*xm+2.273*p2-.984*p3
      pm1(16)=1.240-1.297*xm
      pm1(17)=-.490+1.181*xm-1.038*p2
      pm1(18)=1.858*(xm**2.033+0.745*xm**1.065)
      pm1(19)=1-sqrt(1-(xm-.05)**2)
 
      return
      end
 
c ------------------------------------------------------------------- c
 
      real function sigmaneu(in,x,n,ab_,ab0,xnor)
c     bound-free cross section for neutral matter with cosmic abundances
c     per hydrogen atom in units of the Thomson cross section
c     ab_ - log10 of abundance of elements heavier than He
c           relative to the Solar abundances
c     ab0 - iron log10 abundance relative to hydrogen
c     x - energy vector of the length n; simgabf is calculated for
c     element in from x
      implicit none  
      integer nem 
      PARAMETER(nem=5000) 
      integer i,in,n,nold  
      real x(nem),xold(nem),opac(nem),abo(17),ab(17)
      real ab0,ab_,totlxne,xnor
      real ab_in,ab0in
      logical fl1,fl2,firstflag
      save opac,abo,ab,xold,nold,fl1,fl2,firstflag,ab_in,ab0in
      data xold /nem*0./
      data nold /0/ 
      data ab_in,ab0in /2*999999./
      data firstflag,fl1,fl2 /3*.true./

      integer nu,nuo(17)
      data (nuo(i),i=1,17)
     &     /1,2,6,7,8,10,11,12,13,14,16,17,18,20,24,26,28/
cjp 
c for XSPEC 9 
c      include '../../../../../spectral/inc/sitedef.inc'
c for XSPEC 10 when local directory in $XANADU/local 
      include '../spectral/inc/sitedef.inc' 
      character*72 absfil
      integer lenact,ilun,ios
 
      if(firstflag)then
         firstflag=.false.
         call getlun(ilun)
c          ilun=14 
c this is Anders Grevese  
         absfil = XSPECDR(:lenact(XSPECDR))//'abun.dat'
c this is Anders Ebihara 
c PNS reflection function are for these abundances  
c         absfil = XSPECDR(:lenact(XSPECDR))//'abun_an_eb.dat'
c         absfil ='abun.dat'
c         absfil ='abun_an_eb.dat'
c        OPEN(ilun,file=absfil,form='formatted',status='old')
         call openwr(ilun, absfil, 'old', ' ', ' ', 0, 0, ios)
            read(ilun,*,iostat=ios)
            do i=1,17
               read(ilun,*) nu,abo(i)
               if(nu.ne.nuo(i))then
                  write(*,*) 'inadequate abun.dat file'
                  stop
               endif
            enddo
         close(unit=ilun)
         call frelun(ilun)
         ab(1)=abo(1)
         ab(2)=abo(2)
      endif
c     if the energy vector changed from the
c     previous call, recompute the opacity (in totlxs) and reset xold
      if((x(in).ne.xold(in)).or.(n.ne.nold)) then
         nold=n
         do 11 i=1,n
            xold(i)=x(i)
 11      continue
         fl2=.true.
      endif
      if(ab_.ne.ab_in)then
         do 12 i=3,17
            ab(i)=abo(i)+ab_
 12      continue
         ab(16)=ab(16)+ab0
         ab_in=ab_
         ab0in=ab0
         fl1=.true.
      endif
      if(ab0.ne.ab0in)then
         ab(16)=abo(16)+ab_+ab0
         ab0in=ab0
         fl1=.true.
      endif
c     If either the abundances or the energy vector changed,
c     recompute the opacity
      if(fl1.or.fl2) then
         do 10 i=1,n
            opac(i)=totlxne(ab,i,x,n,fl1,fl2,xnor)
 10      continue
      endif

      sigmaneu=opac(in)
 
      return
      end
 
C----------------------------------------------------------------------
C
      FUNCTION totlxne( AB ,in,x,n,fl1,fl2,xnor)
c     *** this version with fully ionized H and He!  AAZ, 18.05.94
c     Optimized version: calculetes cross section for abundances
c     contained in the vector AB at energy in 511.keV units given as
c     an element numbered in of the vector x; n is the length of x.
c     Flag fl1=.true. must indicate changes of AB, and fl2=.true.
c     changes of the energy vector x. For energy >10.keV the law
c     xnor*e^-3 is used, xnor is an output parameter.
c     Note: result in sigma_Thomson units.
C
C     Reference:
C      Monika Balucinska-Church and Dan McCammon
C      "Photoelectric Absorption Cross Sections with Variable Abundances"
C      Ap.J. 400, 699 (1992)
C
C
C     Description:
C     Calculates the effective absorption cross section in units of
c     sigma_Thomson per hydrogen atom at energy E in eV for the
c     abundances of the elements specified in vector AB
C
C     Method:
C     Calls seventeen functions that calculate the mass absorption coeffs
C     in cm**2/g for the following elements: H, He, C, N, O, Ne, Na, Mg,
C     Al, Si, S, Cl, Ar, Ca, Cr, Fe, Ni.
C
C     Deficiencies:
C     Works only in the range of energy from 30 eV to 10,000 eV.
C
C     Bugs:
C     None known -- please report any problems to authors
C
C     Authors:
C     Dan McCammon               (47413::MCCAMMON)
C     Monika Balucinska-Church   (19775::MBC)
C
C     History:
C     8.8.91  - original
C     8.1.92  - modified to remove VAX fortran77 extensions
C     12.3.93 - function HELIUM has been replaced by HELIUM_V2 with
C               Fano line profile
C     21.9.93 - Names of functions simplified: TOTLXS_V2 to TOTLX2,
C               HELIUM_V2 to HEL2.
C     24.1.94 - Names of functions simplified: TOTLX2 to TOTLXS,
C               HEL2 to HELIUM.
C
C     Parameters:
C     E - energy in eV
C     AB - vector of length 17 of Log10 abundances relative to
C          hydrogen = AB(1).  N.B. The order of the elements in AB
C          is as follows: H,He,C,N,O,Ne,Na,Mg,Al,Si,S,Cl,Ar,Ca,Cr,Fe,Ni.
C     TOTLX2 - effective cross section in cm**2/H atom for assumed
C             abundances of cosmic elements.  Note that this will be per
C             hydrogen atom for the relative abundances given in AB.
C
C     Type Definitions:
      IMPLICIT NONE
      integer i,nem 
      PARAMETER(nem=5000) 
      real xnor,xnors,xnor2,xnors2 
C
C     Local variables:
      INTEGER K
C          ( index through elements)
      REAL A(17,nem),x(nem),toth(nem),ap(17),ap1(17) 
      integer in,n
c          x  - energy points
c          in - index for x
c          n  - number of energy points
C          ( mass absorption cross sections in cm**2/g
C           for the cosmic elements -- in the same order as for AB)
C     Import:
      REAL E
C          ( energy in eV)
      REAL AB(17)
C          ( log10 abundances relative to hydrogen)
 
C     Export:
      REAL totlxne
C          ( effective cross section in cm**2/H atom)
C     External functions: ( mass absorption coefficients for the elements)
      REAL HYDRO1
      REAL HELIUM1
      REAL CARBON1
      REAL NITRO1
      REAL OXYGEN1
      REAL NEON1
      REAL SODIUM1
      REAL MAGNES1
      REAL ALUM1
      REAL SILICN1
      REAL SULFUR1
      REAL CHLORN1
      REAL ARGON1
      REAL CALC1
      REAL CHROM1
      REAL IRON1
      REAL NICKEL1
 
 
C     Local constants:
      REAL AW(17)
C          ( atomic weights of the elements)
      REAL AV
C          ( Avogadro's number) <- multiplied by sigma_Thomson
      INTEGER NUMB
C          ( number of elements)
      integer inh
c          ( limiting bin for asymptotic absorption )
 
      logical fl0,fl1,fl2,fl4 
      save a,aw,av,numb,ap,ap1,fl0,fl4,xnors,xnors2   
 
      DATA AW /1.00797, 4.0026, 12.01115, 14.0067, 15.9994, 20.183,
     &         22.9898, 24.312, 26.9815, 28.086, 32.064, 35.453,
     &         39.94, 40.08, 51.996, 55.847, 58.71/
 
c      DATA AV /6.022045E23/
      DATA AV /.400614/
 
      DATA NUMB /17/
 
      data fl0/.true./,fl4/.true./
 
C     Start:
 
c     normalizing factor for the asymptotic relation
      if(fl1)then
C        mass absorption coefficients for the elements at 10.keV
         if(fl0)then
            E=10000.
c            ap( 1)=HYDRO1(E)
c            ap( 2)=HELIUM1(E)
            ap( 1)=0
            ap( 2)=0
            ap( 3)=CARBON1(E)
            ap( 4)=NITRO1(E)
            ap( 5)=OXYGEN1(E)
            ap( 6)=NEON1(E)
            ap( 7)=SODIUM1(E)
            ap( 8)=MAGNES1(E)
            ap( 9)=ALUM1(E)
            ap(10)=SILICN1(E)
            ap(11)=SULFUR1(E)
            ap(12)=CHLORN1(E)
            ap(13)=ARGON1(E)
            ap(14)=CALC1(E)
            ap(15)=CHROM1(E)
            ap(16)=IRON1(E)
            ap(17)=NICKEL1(E)
            fl0=.false.
         endif
         xnors=0.
         do 30 k=1,numb
            xnors=xnors+AW(K)*ap(K)*(10.**(AB(K)-AB(1)))
 30      continue
         xnors=xnors/av*(0.01947)**3

c     normalizing factor for the asymptotic relation
C     mass absorption coefficients for the elements at 30 eV
         if(fl4)then
            E=30.
            ap1( 1)=0.
            ap1( 2)=0.
            ap1( 3)=CARBON1(E)
            ap1( 4)=NITRO1(E)
            ap1( 5)=OXYGEN1(E)
            ap1( 6)=NEON1(E)
            ap1( 7)=SODIUM1(E)
            ap1( 8)=MAGNES1(E)
            ap1( 9)=ALUM1(E)
            ap1(10)=SILICN1(E)
            ap1(11)=SULFUR1(E)
            ap1(12)=CHLORN1(E)
            ap1(13)=ARGON1(E)
            ap1(14)=CALC1(E)
            ap1(15)=CHROM1(E)
            ap1(16)=IRON1(E)
            ap1(17)=NICKEL1(E)  
            fl4=.false.
         endif

         xnors2=0.
         do 35 k=1,numb
            xnors2=xnors2+AW(K)*ap1(K)*(10.**(AB(K)-AB(1)))
 35      continue
         xnors2=xnors2/av*(5.87084e-5)**3

         fl1=.false.
      endif
      xnor=xnors
      xnor2=xnors2 
 
      if(fl2) then
         do 20 i=1,n
c           above 10 keV, use e**-3 law
            if(x(i).lt.0.01947) then
               E=x(i)*511000
c               A( 1,i)=HYDRO1(E)
c               A( 2,i)=HELIUM1(E)
               A( 1,i)=0
               A( 2,i)=0
               A( 3,i)=CARBON1(E)
               A( 4,i)=NITRO1(E)
               A( 5,i)=OXYGEN1(E)
               A( 6,i)=NEON1(E)
               A( 7,i)=SODIUM1(E)
               A( 8,i)=MAGNES1(E)
               A( 9,i)=ALUM1(E)
               A(10,i)=SILICN1(E)
               A(11,i)=SULFUR1(E)
               A(12,i)=CHLORN1(E)
               A(13,i)=ARGON1(E)
               A(14,i)=CALC1(E)
               A(15,i)=CHROM1(E)
               A(16,i)=IRON1(E)
               A(17,i)=NICKEL1(E)
               inh=i
            else
               toth(i)=x(i)**(-3)
            endif
 20      continue
         inh=inh+1
         fl2=.false.
      endif
 
      if(x(in).lt..01947.and.x(in).gt.5.87084e-5)then
         totlxne=0.
C        loop through elements
         DO 10 K=1,NUMB
            totlxne=totlxne+AW(K)*A(K,in)*(10.**(AB(K)-AB(1)))
10       CONTINUE
         totlxne=totlxne/av
      elseif(x(in).ge.0.01947) then
         totlxne=xnor*toth(in)
      else 
         totlxne=xnor2*toth(in)
      endif
 
      RETURN
      END
 
c ------------------------------------------------------------------ c
 
      subroutine greeion(x,spref,xm,ab_,ab0,spp,gamma,ec,temp,xi,jmax)
c     calculates the reflected spectrum
c     ab0 - iron log10 abundance relative to hydrogen
      implicit none 
      integer nem
      PARAMETER(nem=5000)
      REAL ab0,ab_,ddy,dy,dym,Ec,fjc,Gamma,GNRJ1,GRXYH1,
     &     GRXYL1,SIGMAION,SPP,sr,srp,Temp,Xi,xjd,xjl,
     &     xkabs,Xm
      REAL xmin,xmo,xnor,xrefmax,xtransh,xtransl,y,y0,ymin 
      INTEGER i,j,Jmax,jmin,jrefmax,jtransh,jtransl,nc,ncp 
      real x(nem),spref(nem)
      real pm1(19),pmx(26),ap(3)
      external spp
 
c  precision factor for Green's function integration
      ncp=100
c  ranges for different methods
      xmin=2.e-4
      xtransl=0.01957
      xtransh=0.02348
      ymin=.03333
c  angle dependent parameters
      xmo=max(.05,min(.95,xm))
      call pm1yn(xmo,pm1)
      call pmxyn(xmo,pmx)
      ap(2)=0.802-1.019*xmo+2.528*xmo**2-3.198*xmo**3+
     &      1.457*xmo**4+(.030/xmo)**4
      ap(3)=0.381
      xrefmax=1/(ymin+pmx(26))
 
      jmin=0
      do 50 j=1,jmax
         if(x(j).gt.xmin) goto 51
         jmin=j
 50   continue
 51   jtransl=0
      do 52 j=max(1,jmin),jmax
         if(x(j).gt.xtransl) goto 53
         jtransl=j
 52   continue
 53   jtransh=0
      do 54 j=max(1,jtransl),jmax
         if(x(j).gt.xtransh) goto 55
         jtransh=j
 54   continue
 55   jrefmax=0
      do 56 j=max(1,jtransh),jmax
         if(x(j).gt.xrefmax) goto 57
         jrefmax=j
 56   continue
 57   continue

      do 10 j=1,jmin
         spref(j)=0
 10   continue
c-----------------------
      do 20 j=jmin+1,jtransl
         xkabs=sigmaion(j,x,jtransh,gamma,temp,xi,ab_,ab0,xnor)
         spref(j)=spp(1/x(j),gamma,ec)*gnrj1(xkabs,xmo)
 20   continue
c-----------------------
      if(jmin.ge.jtransl) then
         xkabs=sigmaion(1,x,max(1,jtransh),gamma,temp,xi,ab_,ab0,xnor)
      endif
      ap(1)=xnor/1.21
      xjl=alog(xtransl)
      xjd=alog(xtransh)-xjl
      do 60 j=jtransl+1,jtransh
         fjc=.5*sin(3.14159*((alog(x(j))-xjl)/xjd-.5))
         xkabs=sigmaion(j,x,jtransh,gamma,temp,xi,ab_,ab0,xnor)
         y=1/x(j)
         spref(j)=spp(y,gamma,ec)*gnrj1(xkabs,xmo)
         dym=y-ymin
         dy=min(2,dym)
         ddy=(dy-pmx(26))/(ncp+1)
         y0=y-dy
         sr=spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
         dy=pmx(26)
         y0=y-dy
         sr=.5*(sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap))
         do 62 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
 62      continue
         sr=sr*ddy
         if(dym.gt.2.)then
            ddy=dym-2
            nc=int(ncp*ddy/(dym-pmx(26)))
            ddy=ddy/(nc+1)
            dy=dym
            y0=y-dy
            srp=spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
            dy=2
            y0=y-dy
            srp=.5*(srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap))
            do 64 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
 64         continue
            sr=sr+srp*ddy
         endif
         spref(j)=(.5-fjc)*spref(j)+(.5+fjc)*sr*x(j)
 60   continue
c-----------------------
      do 30 j=jtransh+1,jrefmax
         y=1/x(j)
         dym=y-ymin
         dy=min(2,dym)
         ddy=(dy-pmx(26))/(ncp+1)
         y0=y-dy
         sr=spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
         dy=pmx(26)
         y0=y-dy
         sr=.5*(sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap))
         do 32 i=1,ncp
            dy=dy+ddy
            y0=y-dy
            sr=sr+spp(y0,gamma,ec)*grxyl1(pm1,pmx,dy,y0,ap)
 32      continue
         sr=sr*ddy
         if(dym.gt.2.)then
            ddy=dym-2
            nc=int(ncp*ddy/(dym-pmx(26)))
            ddy=ddy/(nc+1)
            dy=dym
            y0=y-dy
            srp=spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
            dy=2
            y0=y-dy
            srp=.5*(srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap))
            do 34 i=1,nc
               dy=dy+ddy
               y0=y-dy
               srp=srp+spp(y0,gamma,ec)*grxyh1(pm1,pmx,dy,y0,ap)
 34         continue
            sr=sr+srp*ddy
         endif
         spref(j)=sr*x(j)
 30   continue
c----------------------
      do 40 j=jrefmax+1,jmax
        spref(j)=0
 40   continue
 
      return
      end
 
c ------------------------------------------------------------------ c
c  Reflection factor for nonrelativistic Compton reflection,
c  reflection orders >1 approximated with isotropic phase function,
c  and approximation of H-function by Basko ApJ 223,268,(1978),
c  valid up to 14keV;
c  sig - sigma_bf per hydrogen atom in sigma_Thomson units
c  xm  - cosine of the angle between normal to the slab and the
c        observation direction
 
      real function    gnrj1(sig,xm)
      implicit         none
      real             sig,xm
      real             al0,al1,al2,xl,x2,x4,sq3,xq3
      parameter        (sq3=1.732051)
 
      al0=1/(1+sig/1.21)
      al1=sqrt(1-al0)
      al2=sq3*al1
      xl=log(1+1/xm)
      x2=xm*xm
      x4=x2*x2

      if(al0.lt.0.999) then 
        gnrj1=.5*al0*xm*(.375*((3-2*x2+3*x4)*xl+(3*x2-1)*(.5-xm))+
     & ((1+(1/al1-1)*(1-log(1+al2)/al2))*(1+xm*sq3)/(1+xm*al2)-1)*xl)
      else
        xq3=1+xm*sq3
        gnrj1=.5*al0*xm*(.375*((3-2*x2+3*x4)*xl+(3*x2-1)*(.5-xm))+
     & ((1+sq3/2)*xq3*(1-xq3*sqrt(1-al0))-1)*xl)
      endif
      return
      end
 
c ------------------------------------------------------------------- c
 
      function sigmaion(in,x,n,gamma0,temp0,xi0,ab_0,ab0,asc0)
c     bound-free cross section for ionized matter with cosmic abundances
c     and variable iron abundance (in units of solar iron abundance)
c     per hydrogen atom in units of the Thomson cross section
c     e is in units of 511 keV
      implicit none  
      integer nem 
      PARAMETER(nem=5000) 
      REAL ab,ab0,ab_,ab_0,asc,asc0,gamma,gamma0,SIGMAION,temp,
     &     Temp0,WARMABSJ,xas,xi,Xi0
      INTEGER i,In,n,nold 
      real x(nem),opac(nem),xold(nem)
      logical fl1,firstflag
      save opac,gamma,temp,xi,ab_,ab,asc,xold,nold,fl1,firstflag
      data xold/nem*0./
      data nold/0/
      data gamma,temp,xi,ab_,ab/0.,0.,0.,99999.,99999./
      data firstflag,fl1/2*.true./

c In this version, changes of ear during one xspec session are allowed.
      if((gamma.ne.gamma0).or.(temp.ne.temp0).or.
     &   (xi.ne.xi0).or.(ab_.ne.ab_0).or.(ab.ne.ab0))then
         gamma=gamma0
         temp=temp0
         xi=xi0
         ab=ab0
         ab_=ab_0
         firstflag=.true.
         fl1=.true.
         xas=.01947
         asc=1.503e2*
     &       warmabsj(xas,gamma,temp,xi,ab_,ab,firstflag)*xas**3
      endif

      if((x(in).ne.xold(in)).or.(n.ne.nold))then
         nold=n
         do i=1,n
            xold(i)=x(i)
         enddo
         fl1=.true.
      endif

      if(fl1)then
         do 10 i=1,n
            opac(i)=1.503e2*
     &              warmabsj(x(i),gamma,temp,xi,ab_,ab,firstflag)
 10      continue
         fl1=.false.
      endif

      asc0=asc
      sigmaion=opac(in)
 
      return
      end

c ------------------------------------------------------------------- c

	function warmabsj(x0,gamma,temp,xi,ab_,ab0,firstflag)
c 	calculate cross-sections using Daltabuit and Cox, use cosmic 
c       abundances, and radiative recombination rates from 
c       Aldrovandi and Pequignot 1973 and explicit formula 
c       for recombination to hydrogenlike atoms
c       result in units of 10^-22
        implicit none
        real     warmabsj,x0,gamma,temp,xi,ab_,ab0 
        real     t4,z2,e1,e2,re,rp,ab_d,it4 
	real     aw(20),abo(20),ab(20),ion(20,30,10),y
        real     int(20,30),arec(20,30)
	real     denom,xp,sigmaeff,e(800),sigma(20,30,800)
	real*8   num(20,30),sum,ratio(20,30),mult,mul(30),xil
	real*8   dnum,dmul 
	integer  z(20),i,j,imax,n,ios,k,iab0
	logical  firstflag
c
c     ******new bit for xspec
cjp
c for XSPEC 9
c      include '../../../../../spectral/inc/sitedef.inc'
c for XSPEC 10 when local directory in $XANADU/local 
      include '../spectral/inc/sitedef.inc'
      character*72 absfil
      integer lenact,ilun10,ilun11  
c

c ------------------------------------------------------------------- c
c  set up constants and read in data if first time round

        save z,aw,abo,ab,ion,e,sigma,num,iab0
        data z(1)/0./

	if(z(1).eq.0) then
           iab0=0
        call getlun(ilun10)
        absfil = XSPECDR(:lenact(XSPECDR))//'iond_mm.dat'
cjp
        call openwr(ilun10, absfil, 'old', ' ', ' ', 0, 0, ios)
c          open(unit=10,file='iond_mm.dat',status='old',iostat=ios)
c          open(unit=11,file='mansig.dat',status='old')
c          ilun10=10
c          ilun11=11

	   i=1
	   read(ilun10,*,iostat=ios) z(i),aw(i),abo(i)
           if(z(i).eq.26) iab0=i
c  NB j=1 is ground state, j=z(i)=fully charged
	   do while (ios.eq.0)
	      do j=1,z(i),1
		 read(ilun10,*) ion(i,j,1),ion(i,j,2),ion(i,j,3),
     &               ion(i,j,4),ion(i,j,5),ion(i,j,6),ion(i,j,7),
     &               ion(i,j,8),ion(i,j,9),ion(i,j,10)
c  change units of coefficients
                 ion(i,j,2)=ion(i,j,2)*1.e+10
                 ion(i,j,4)=ion(i,j,4)*1.e+04
                 ion(i,j,5)=ion(i,j,5)*1.e-04
                 ion(i,j,7)=ion(i,j,7)*1.e-04
	      enddo
	      i=i+1
	      read(ilun10,*,iostat=ios) z(i),aw(i),abo(i)
              if(z(i).eq.26) iab0=i
	   enddo
	   imax=i-1

c       ******new bit for xspec
        close(unit=ilun10)
        call frelun(ilun10)

c  check not overflowing array bounds
           if(imax.gt.20.or.z(imax).gt.30) then 
              write(*,*) 'element array overflowing array bounds'
              write(*,*) 'number of elements= ',imax
              write(*,*) 'maximum atomic number= ',z(imax)
              stop
           end if
c  read in photoionization data from Reilman and Manson ApJ supp 40 1979
        call getlun(ilun11)
        absfil = XSPECDR(:lenact(XSPECDR))//'mansig.dat'
        call openwr(ilun11, absfil, 'old', ' ', ' ', 0, 0, ios)

           do i=1,imax,1
c  now go over all ion states - values come from Reilman
c  ApJ supp 1979 except hydrogenic which is analytic
              do j=1,z(i),1
c  ion identification header
                 read(ilun11,*) 
c  read in over all energies 
c  the table is valid for 5.eV-19907.7eV
                 do k=1,721,1
                    read(ilun11,*) e(k),sigma(i,j,k)
c  change units of sigma
                    sigma(i,j,k)=sigma(i,j,k)/6.6e-27
                 enddo
              enddo
           enddo
c       ******new bit for xspec
        close(unit=ilun11)
        call frelun(ilun11)

c           close(ilun10)
c           close(ilun11)

	endif

c  end of data set up
c ------------------------------------------------------------------- c

        if(firstflag) then
c  insert external iron abundance (ab0)
           ab_d=10.**ab_
           do i=1,imax
              if(z(i).gt.2) ab(i)=abo(i)*ab_d
           enddo
           if(iab0.gt.0) ab(iab0)=ab(iab0)*10.**ab0
           if(xi.gt.1.e-20)then
c  normalized temperature
              t4=temp*.0001
c  calculate recombination coefficients - radiative and dielectric
c  from Aldrovandi+Pequignot 1973 Astro, astrophys 25 137
              do i=1,imax,1
                 do j=1,z(i)-1,1
                    it4=ion(i,j,5)/t4
                    e1=0. 
                    if(it4.lt.40.) e1=exp(-it4) 
                    it4=ion(i,j,7)/t4
                    e2=0. 
                    if(it4.lt.40.) e2=exp(-it4) 
                    arec(i,j)=ion(i,j,2)*(t4**(-ion(i,j,3)))+
     &	               ion(i,j,4)*(t4**(-1.5))*e1*(1.+ion(i,j,6)*e2)
                 enddo
c  do radiative recomb to hydrogenic ion NB assumed gaunt factor=1
c  from Gould+Thakur 1970 ann phys 61 351.
                 z2=z(i)**2
                 y=15.8*z2/t4
                 arec(i,z(i))=1.033e-3*z2*(1.735+log(y)+1/(6.*y))
     &              /sqrt(t4)
              enddo
c  first normalise photon spectrum 
              denom=0.
              do k=1,720,1
                 denom=denom+(e(k)**(1-gamma))*(e(k+1)-e(k))
              enddo
              denom=1/denom
c  calculate the photoionization integrals and population balance
              do i=1,imax,1
                 do j=1,z(i),1
                    int(i,j)=0.0
c  do integral in log10 (eV) from lowest energy
c  to the maximum energy- use same energy grid
c  as the cross-sections are defined on
                    do k=1,720,1
                       int(i,j)=int(i,j)+
     &                    sigma(i,j,k)*(e(k)**(-gamma))*(e(k+1)-e(k))
                    enddo
                    int(i,j)=int(i,j)*denom
c  const=7.958d-2/2.43d+4 
		    ratio(i,j)=dlog(3.2749d-6*int(i,j)/arec(i,j))
	         enddo
	      enddo
c  calculate population numbers
c  set sum=1 as the sum=1+A21+A32A21+....
c  complicated because floating range limit 
              do i=1,imax,1
                 xil=dlog(dble(xi))
                 mult=0.d0
                 do j=z(i),1,-1
                    mul(j)=0.d0
                    do n=j,1,-1
                       mul(j)=mul(j)+ratio(i,n)
                    enddo
                    mul(j)=mul(j)+j*xil
                    if(mul(j).gt.mult) mult=mul(j)
                 enddo
                 do j=z(i),1,-1
                    mul(j)=mul(j)-mult
                 enddo
                 sum=0.d0
                 do j=z(i),1,-1
                    dmul=0.d0 
                    if(dabs(mul(j)).lt.7d1) dmul=dexp(mul(j)) 
                    sum=sum+dmul
                 enddo
                    dmul=0.d0 
                    if(dabs(mult).lt.7d1) dmul=dexp(-mult) 
                 sum=sum+dmul
                 num(i,1)=-mult-dlog(sum)
                 do j=2,z(i)+1,1
                    num(i,j)=num(i,j-1)+ratio(i,j-1)+xil
                 enddo
                 do j=1,z(i)+1,1
                    dnum=0d0 
                    if(dabs(num(i,j)).lt.7d1) dnum=dexp(num(i,j))
                    num(i,j)=dnum
                 enddo
	      enddo				
           else
              do i=1,imax,1
                 num(i,1)=1.
                 do j=2,z(i)+1,1
                    num(i,j)=0.
                 enddo
              enddo
           endif
           firstflag=.false.
        endif
c ------------------------------------------------------------------- c

        xp=x0*5.11e5
	sigmaeff=0.
        if(xp.lt.e(721)) then
c  pick closest energy bin
           k=1
           dowhile(e(k).lt.xp)
              k=k+1
           enddo
           re=(xp-e(k-1))/(e(k)-e(k-1))
           do i=1,imax,1
              do j=1,z(i),1
                 sigmaeff=sigmaeff+sngl(num(i,j))*ab(i)*
     &              (sigma(i,j,k-1)+re*(sigma(i,j,k)-sigma(i,j,k-1)))
              enddo
           enddo
        else
           re=(xp/e(721))**(-3)
           do i=1,imax,1
              rp=ab(i)*re
              do j=1,z(i),1
                 sigmaeff=sigmaeff+sngl(num(i,j))*sigma(i,j,721)*rp
              enddo
           enddo
        endif
        sigmaeff=sigmaeff*6.6e-5

        warmabsj=sigmaeff

	return
	end

c ------------------------------------------------------------------- c

C-------------------------------------------------------------------------
C
C  File XSCTNS.FOR
C
C        This set of subroutines calculates the photoelectric absorption
C   cross sections for the elements H, He, C, N, O, Ne, Na, Mg, Al, Si,
C   S, Cl, A, Ca, Cr, Fe, and Ni.  The result is in cm**2/g, given the
C   photon energy in eV.  These functions are valid only over the energy
C   range 30 - 10,000 eV, but do NOT check that the input energy is
C   within the valid range.  These functions are called by TOTLX2 to
C   calculate the total effective cross section, given a set of relative
C   abundances.  They can also be used by themselves.
C
C   history : 26 May, 1993 - original HELIUM function replaced by HEL2 ;
C              HELIUM not deleted from the program (19775::MBC)
C
C           : 21 Sep, 1993 - Function HELIUM now deleted, all remnants
C             of VAX Fortran 77 extensions removed, file renamed
C             XSCTN2 (19775::MBC)
C
C           : 22 Sep, 1993 - bug in FANO routine (called by HEL2) corrected
C             that resulted in incorrect line shape and energy (19775::MBC)
C
C           : 24 Jan, 1994 - name HEL2 changed to HELIUM
C   References:
C
C      Monika Balucinska-Church and Dan McCammon
C      "Photoelectric Absorption Cross Sections with Variable Abunances"
C      Ap.J. 400, 699 (1992)
C
C      All data are from:
C      B. L. Henke, P. Lee, T. J. Tanaka, R. L. Shimabukuro and B. K.
C      Fujikawa, 1982, Atomic Data and Nuclear Data Tables, vol 27, p 1.
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C  Real Function: ALUM
C  Source: Atomic & Nuclear Data Tables, January 1982
C
C  Description:
C      Calculates mass absorption coefficient (mu/rho) for aluminum.
C  History:  updated below L-edge (72.78 eV) - Aug 30, 1991 (MBC)
C
C  Usage:  FUNCTION ALUM(E)
C      E = Energy in eV.
C      ALUM returns mu/rho in cm**2/gm.
C
C  COMMON Blocks:
C      none
C  IMPLICIT
C      none
C  Subroutines called by ALUM:
C      none
C
C--------------------------------------------------------------------------
      FUNCTION ALUM1(E)
 
      REAL E, Elog, X, ALUM1
 
        Elog = Alog(E)
        If(E.LT.72.78) Then
            X = 26.90487 + (3. - 9.135221)*Elog + 1.175546*Elog*Elog
 
        Else If(E.LT.1559.9) Then
                X = -38.1232 + 29.5161*Elog - 4.45416*Elog*Elog +
     +                        0.226204*Elog*Elog*Elog
            Else
                X = 14.6897 + 4.22743*Elog - 0.344185*Elog*Elog +
     +                        8.18542E-3*Elog*Elog*Elog
        Endif
      X=X-3.*Elog
      ALUM1 = 0.
      IF(abs(X).lt.70.) ALUM1 = Exp(X) 
 
      RETURN
      END
C------------------------------------------------------------------------
C------------------------------------------------------------------------
C------------------------------------------------------------------------
C
C  REAL FUNCTION: ARGON
C  Source: Atomic & Nuclear Data Tables, January 1982
C  Description: ARGON calculates the mass absorbtion coefficient of argon.
C
C ****          works well for energies above 40 eV !!!
C
C  History: updated below L-edge (245 eV) - Aug 30, 1991 (MBC)
C
C  Usage: FUNCTION ARGON(E)
C      E = Energy in eV
C      ARGON returns the mass absorption cross section in cm**2/g
C
C  COMMON BLOCKS:
C      NONE
C
C  IMPLICIT
C      NONE
C  SUBROUTINES CALLED BY ARGON:
C      NONE
C
C---------------------------------------------------------------------------
      FUNCTION ARGON1(E)
 
      REAL E, Elog, X, ARGON1
 
        Elog = Alog(E)
        If(E.LT.245.0) Then
           X = -330.3509 + (267.7433 + 3.)*Elog - 78.90498*Elog*Elog
     &         + 10.35983*(Elog**3) - 0.5140201*(Elog**4)
 
            Else If(E.LT.3202.9) Then
               X = -5.71870 + (8.85812*Elog) + (-0.307357*Elog*Elog) +
     +             (0.00169351*(Elog**3)) + (-0.0138134*(Elog**4)) +
     +             (0.00120451*(Elog**5))
            Else
               X = 19.1905 + (2.74276*Elog) + (-0.164603*Elog*Elog) +
     +             (0.00165895*Elog*Elog*Elog)
        Endif
      X=X-3.*Elog
      ARGON1 = 0.
      IF(abs(X).lt.70.) ARGON1 = Exp(X) 
 
      Return
      End
C--------------------------------------------------------------------------
C------------------------------------------------------------------------------
C------------------------------------------------------------------------------
C     Real Function: CALC
C     Source: Atomic & Nuclear Data Tables, January 1982
C
C     Description:
C     Calculates mass absorption coefficient (mu/rho) for calcium
C
C     History: original Aug 6, 1991 (MBC)
C
C     Usage:   FUNCTION CALC(E)
C              E = Energy in eV
C              CALC returns mu/rho in cm**2/g
C
C     COMMON Blocks:
C            none
C
C     IMPLICIT
C            none
C-------------------------------------------------------------------------
C
      FUNCTION CALC1(E)
C
C
      REAL E,ELOG,X,CALC1
C
C
      Elog = Alog(E)
 
 
      IF(E.LT.349.31) THEN
      X=-873.972 + (865.5231 + 3.)*Elog - 339.678*Elog*Elog +
     &   66.83369*(Elog**3) - 6.590398*(Elog**4) +
     &    0.2601044*(Elog**5)
 
            ELSE IF(E.LT.4038.1) THEN
            X=-3449.707 + (2433.409 + 3.)*Elog
     &        - 682.0668*Elog*Elog + 95.3563*(Elog**3)
     &        - 6.655018*(Elog**4) + 0.1854492*(Elog**5)
 
                ELSE
                X=18.89376 + (3. - 0.2903538)*Elog
     &            - 0.1377201*Elog*Elog
 
      ENDIF
      X=X-3.*Elog
      CALC1 = 0.
      IF(abs(X).lt.70.) CALC1 = Exp(X)
 
 
      RETURN
      END
C
C--------------------------------------------------------------------------
C--------------------------------------------------------------------------
C
C  REAL FUNCTION: CARBON
C  Source: Atomic & Nuclear Data Tables, Jan. 1982
C
C  Description: CARBON calculates the mass absorbtion cross-section of carbon.
C
C  USAGE: FUNCTION CARBON(E)
C      E = Energy in EV
C      CARBON returns the mass absorption cross-section in cm**2/g
C
C  COMMON BLOCKS:
C      NONE
C
C  IMPLICIT
C      NONE
C
C  SUBROUTINES CALLED BY CARBON:
C      NONE
C
C-------------------------------------------------------------------------
      FUNCTION CARBON1(E)
C
C
      REAL E,ELOG,X,CARBON1
C
C
        Elog = Alog(E)
        If(E.LT.284.0) Then
            X = 8.74161 + (7.13348*Elog) + (-1.14604*Elog*Elog) +
     +          (0.0677044*Elog*Elog*Elog)
        Else
            X = 3.81334 + (8.93626*Elog) + (-1.06905*Elog*Elog) +
     +          (0.0422195*Elog*Elog*Elog)
        Endif
      X=X-3.*Elog
      CARBON1 = 0.
      IF(abs(X).lt.70.) CARBON1 = Exp(X) 
      Return
      End
C
C---------------------------------------------------------------------------
C---------------------------------------------------------------------------
C     Real Function: CHLORN
C     Source: Atomic & Nuclear Data Tables, January 1982
C
C     Description:
C     Calculates mass absorption coefficient (mu/rho) for chlorine
C
C     History: original Aug 6, 1991 (MBC)
C
C     Usage:   FUNCTION CHLORN(E)
C              E = Energy in eV
C              CHLORN returns mu/rho in cm**2/g
C
C     COMMON Blocks:
C            none
C
C     IMPLICIT
C            none
C-------------------------------------------------------------------------
C
      FUNCTION CHLORN1(E)
C
C
      REAL E,ELOG,X,CHLORN1
C
C
      Elog = Alog(E)
 
 
      IF(E.LT.202.0) THEN
      X=6253.247 +(3. - 8225.248)*Elog + 4491.675*Elog*Elog
     &  - 1302.145*(Elog**3) + 211.4881*(Elog**4)
     &  - 18.25547*(Elog**5) + 0.6545154*(Elog**6)
 
            ELSE IF(E.LT.2819.6) THEN
            X=-233.0502 + (143.9776 + 3.)*Elog
     &        - 31.12463*Elog*Elog +
     &        2.938618*(Elog**3) - 0.104096*(Elog**4)
 
                ELSE
                X=-23.74675 + (14.50997 + 3.)*Elog
     &            - 1.857953*Elog*Elog + 6.6208832E-2*(Elog**3)
 
      ENDIF
      X=X-3.*Elog
      CHLORN1 = 0.
      IF(abs(X).lt.70.) CHLORN1 = Exp(X)
 
 
      RETURN
      END
C
C---------------------------------------------------------------------------
C---------------------------------------------------------------------------
C     Real Function: CHROM
C     Source: Atomic & Nuclear Data Tables, January 1982
C
C     Description:
C     Calculates mass absorption coefficient (mu/rho) for chromium
C
C     History: original Aug 5, 1991 (MBC)
C
C     Usage:   FUNCTION CHROM(E)
C              E = Energy in eV
C              CHROM returns mu/rho in cm**2/g
C
C     COMMON Blocks:
C            none
C
C     IMPLICIT
C            none
C-------------------------------------------------------------------------
C
      FUNCTION CHROM1(E)
C
C
      REAL E,ELOG,X,CHROM1
C
C
      Elog = Alog(E)
 
 
      IF(E.LT.598.0) THEN
      X=-0.4919405 + (12.66939 + 3.)*Elog - 5.199775*Elog*Elog +
     &  1.086566*(Elog**3) - 0.1196001*(Elog**4) +
     &  5.2152011E-3*(Elog**5)
 
            ELSE IF(E.LT.691.0) THEN
            X=27.29282 +(3. - 2.703336)*Elog
 
                ELSE IF(E.LT.5988.8) THEN
                X=-15.2525 + (13.23729 + 3.)*Elog
     &            - 1.966778*Elog*Elog + 8.062207E-2*(Elog**3)
 
                     ELSE
                     X=8.307041 + (2.008987 + 3.)*Elog
     &                - 0.2580816*Elog*Elog
 
      ENDIF
      X=X-3.*Elog
      CHROM1 = 0.
      IF(abs(X).lt.70.) CHROM1 = Exp(X)
 
 
      RETURN
      END
C------------------------------------------------------------------------------
C-----------------------------------------------------------------------------
C
C     Real Funcion : HELIUM
C     Source : Marr, G. V., and West, J. B., Atomic and Nuclear Data Tables,
C                (1976) 18, 497.
C             Oza, D. H., (1986), Phys. Rev. A, 33,  824.
C             Fernley, J. A., Taylor, K. T., and Seaton, M. J., (1987),
C                J. Phys. B., 20, 6457.
C
C     Description :
C     calculates mass absorption coefficient (mu/rho) in cm2/g for neutral
C     helium for the given energy in eV.
C     Cross sections come from experimental data compiled by Marr and
C     West (Atomic Data and Nuclear Data Tables (1976) 18, 497).
C     The four strongest autoionization resonances are taken into account;
C     numbers come from Oza (Phys Rev A (1986), 33, 824), and Fernley et al.
C     (J. Phys B (1987) 20, 6457).
C
C     Deficiencies :
C     works in the energy range from 30 eV to 10,000 eV
 
C     Bugs :
C     if any are found please report to the authors
C
C     History :
C     this subroutine replaces the previous version of HELIUM which
C     calculated mass absoprtion coefficients based on Henke's data
C     (Henke, B. L., et al., (1982), Atomic and Nuclear Data Tables, 27, 1).
C     This version of HELIUM returns mass  absorption coefficients which
C     are in better agreement with the best experiments as well as
C     theoretical models (see Chen, W. F., Cooper, G., and Brion, C. E.,
C     (1991), Phys. Rev. A, 44, 186).  This fortran-77 version of the
C     subroutine is based on Jelinsky's program written in C
C     (EUVE Archive)
C
C     History :
C     Jan 4, 1993 : original (19775::MBC)
C     Feb 23, 1993 : comments added and modified to remove VAX
C                    fortran 77 extensions
C
C     Sep 21, 1993 : further remnants of VAX fortran 77 extensions
C                    have been removed (19775::MBC)
C
C     Sep 22, 1993 : bug in the FANO routine has been removed (19775::MBC)
C
C     Usage : FUNCTION HELIUM(E)
C            E = Energy in eV
C
C     Common Blocks :
C           none
C
C     Implicit :
C           none
C
C     Functions called by HELIUM
C           FANO
C
C------------------------------------------------------------------------------
      FUNCTION HELIUM1(E)
 
C    Type definitions :
C     IMPLICIT NONE
C    Global variables :
C    Structure definitions :
C    Function declarations :
      REAL FANO1
C    Local constants :
      INTEGER IP
C         ( index through loop )
      PARAMETER (IP=8)
      INTEGER IF
C         ( index through loop )
      PARAMETER (IF=4)
      REAL AV
C         ( Avogadro's number )
      PARAMETER (AV=6.022045E23)
      REAL AW
C         ( atomic weight of HYDROgen )
      PARAMETER (AW=4.0026E0)
C    Local variables :
       REAL LAMBDA
C          ( wavelength in Angstroms)
       REAL X
       REAL Y
       REAL SIGMA
C          ( cross section in cm2/atom)
       INTEGER I
C          ( index trough loop)
C     Import :
       REAL E
C          ( energy in eV)
C     Export :
       REAL HELIUM1
C          ( cross section in cm**2/g)
C    Local data :
       REAL C1(IP)
       REAL C2(IP)
       REAL Q(IF)
       REAL NU(IF)
       REAL GAMMA(IF)
 
C          ( polynomial coefficients for Marr and West data)
       DATA C1 /-2.953607E1, 7.083061E0, 8.678646E-1, -1.221932E0,
     +         4.052997E-2, 1.317109E-1, -3.265795E-2, 2.500933E-3/
 
C          ( polynomial coefficients for Marr and West data )
       DATA C2 /-2.465188E1, 4.354679E0, -3.553024E0, 5.573040E0,
     +          -5.872938E0, 3.720797E0, -1.226919E0, 1.576657E-1/
 
C          ( parameters Q for resonances (Fernley et al. 1987) )
       DATA Q /2.81E0, 2.51E0, 2.45E0, 2.44E0/
 
C          ( parameters NU for resonances (Oza 1986) )
       DATA NU /1.610E0, 2.795E0, 3.817E0, 4.824E0/
 
C          ( parameters GAMMA for resonances (Oza 1986) )
       DATA GAMMA /2.64061E-3, 6.20116E-4, 2.56061E-4,
     +             1.320159E-4/
 
 
C     Start :
 
       LAMBDA=12398.54E0/E
       X=ALOG10(LAMBDA)
 
       IF(LAMBDA.GT.503.97E0) THEN
         HELIUM1=0.E0
 
       ELSEIF(LAMBDA.LT.46.E0) THEN
         Y=0.E0
         DO 1 I=1,IP
 1         Y=Y + C2(I)*(X**(I-1))
 
       ELSE
         Y=0.E0
         DO 2 I=1,IP
 2         Y=Y + C1(I)*(X**(I-1))
 
         DO 3 I=1,IF
 3          Y=Y + ALOG10(FANO1(Q(I),NU(I),GAMMA(I),LAMBDA))
 
       ENDIF
 
       SIGMA=0.
       IF(abs(Y).lt.30.) SIGMA=10.E0**Y
       HELIUM1=SIGMA*AV/AW
 
       END
C
C----------------------------------------------------------------------------
C----------------------------------------------------------------------------
C    Real Function FANO
C
C    Source :
C             Fernley, J. A., Taylor, K. T., and Seaton, M. J., (1987),
C                J. Phys. B., 20, 6457.
C             Oza, D. H., (1986), Phys. Rev. A, 33,  824.
C
C    Description :
C     returns a Fano line profile for use by HELIUM. The form of
C     this line shape is taken from Fernley, Taylor and Seaton (above).
C
C    Deficiencies :
C     works in the energy range from 30 eV to 10,000 eV
C
C    Bugs :
C    if any are found please report to the authors
C
C    22 Sep, 1993 - bug fixed in calculations of EPSI - devided by 1.807317
C    not added (19775::MBC)
C
C    History :
C     Jan 4, 1993 : original, based on Jelinsky's program written in
C     C (EUVE Archive) - (19775::MBC)
C
C     Feb 23, 1993 : comments added and modified to remove VAX
C     fortran 77 extensions (19775::MC)
C
C     Sep 21, 1993 : further remnants of VAX fortran 77 extensions
C     removed (19775::MBC)
C
C
C    Common Blocks :
C           none
C    Implicit :
C           none
C
C    Functions called by FANO
C           none
C
C----------------------------------------------------------------------------
C
       FUNCTION FANO1(A,B,C,LAMBDA)
 
C    Type definitions :
C     IMPLICIT NONE
C    Global variables :
C    Structure definitions :
C    Function declarations :
C    Local constants :
C    Local variables :
       REAL EPS
C          ( energy in Rydbergs )
       REAL EPSI
       REAL X
C          ( log_10 of wavelength in Angstroms )
C     Import :
       REAL A
C          ( Q coefficient (Fernley et al. 1987) )
       REAL B
C          ( NU coefficient (Oza 1986) )
       REAL C
C          ( GAMMA coefficient (Oza 1986) )
       REAL LAMBDA
C          ( wavelength in Angstroms )
C     Export :
       REAL FANO1
C    Start :
 
       EPS=911.2671E0/LAMBDA
       EPSI=3.0E0 - 1.E0/(B*B) + 1.807317
       X=2.0*(EPS - EPSI)/C
       FANO1=(X-A)*(X-A)/(1.0E0 + X*X)
 
       END
C-----------------------------------------------------------------------
C       REAL FUNCTION HYDRO(E)
C
C       DATE:3/6/84
C       AUTHOR: ANGIE BETKER
C       SOURCE: ATOMIC AND NUCLEAR DATA TABLES, JANUARY 1982
C
C        History: modified: 6/5/87 - J. Bloch - Created F77 Vax/VMS version.
C                 updated Aug 30, 1991 (MBC)
C
C
C       DESCRIPTION: HYDRO CALCULATES MU/RHO FOR HYDROGEN IN CM**2/GM
C
C       COMMON BLOCKS:
C               NONE
C
C       SUBROUTINES CALLED:
C               NONE
C
C       IMPLICIT
C               NONE
C
C---------------------------------------------------------------------------
C
        FUNCTION HYDRO1(E)
C
C
        REAL E, Elog, X, HYDRO1
C
C
C        IF (E.LT.180) HYDRO1 = (10 ** 10.11) * (E ** -3.03)
C        IF ((E.GE.180).AND.(E.LT.750)) HYDRO1=(10**10.54)*(E**-3.23)
C        IF ((E.GE.750).AND.(E.LT.6000)) HYDRO1=(10**10.94)*(E**-3.37)
C        IF (E.GE.6000) HYDRO1 = (10 ** 10.42) * (E ** -3.24)
 
        Elog=Alog(E)
 
        X=21.46941 + (3. - 2.060152)*Elog - 0.1492932*Elog*Elog
     &    + 5.4634294e-3*(Elog**3)

      X=X-3.*Elog
      HYDRO1 = 0.
      IF(abs(X).lt.70.) HYDRO1 = Exp(X) 
 
        RETURN
        END
C
C----------------------------------------------------------------------
C------------------------------------------------------------------------------
C     Real Function: IRON
C     Source: Atomic & Nuclear Data Tables, January 1982
C
C     Description:
C     Calculates mass absorption coefficient (mu/rho) for iron
C     History: original Aug 6, 1991 (MBC)
C
C     Usage:   FUNCTION IRON(E)
C              E = Energy in eV
C              IRON returns mu/rho in cm**2/g
C
C     COMMON Blocks:
C            none
C
C     IMPLICIT:
C            none
C
C-------------------------------------------------------------------------
C
      FUNCTION IRON1(E)
C
C
      REAL E,ELOG,X,IRON1
C
C
      Elog = Alog(E)
 
      IF(E.LT.707.4) THEN
      X=-15.07332 + (18.94335 + 3.)*Elog - 4.862457*Elog*Elog +
     &    0.5573765*(Elog**3) - 3.0065542E-2*(Elog**4) +
     &    4.9834867E-4*(Elog**5)
 
          ELSE IF(E.LT.7111.2) THEN
          X=-253.0979 + (135.4238 + 3.)*Elog - 25.47119*Elog*Elog +
     &      2.08867*(Elog**3) - 6.4264648E-2*(Elog**4)
 
             ELSE
             X=-1.037655 + (4.022304 + 3.)*Elog
     &       - 0.3638919*Elog*Elog
 
      ENDIF
      X=X-3.*Elog
      IRON1 = 0.
      IF(abs(X).lt.70.) IRON1 = Exp(X)
 
 
      RETURN
      END
C-------------------------------------------------------------------------
C-------------------------------------------------------------------------
C     Real Function: MAGNES
C     Source: Atomic & Nuclear Data Tables, January 1982
C
C     Description:
C     Calculates mass absorption coefficient (mu/rho) for magnesium
C     History: original Aug 6, 1991 (MBC)
C
C     Usage:   FUNCTION MAGNES(E)
C              E = Energy in eV
C              MAGNES returns mu/rho in cm**2/g
C
C     COMMON Blocks:
C            none
C
C     IMPLICIT:
C            none
C-------------------------------------------------------------------------
C
      FUNCTION MAGNES1(E)
C
C
      REAL E,ELOG,X,MAGNES1
C
C
      Elog = Alog(E)
 
 
      IF(E.LT.49.45) THEN
      X=7.107172 + (0.7359418 + 3.)*Elog
 
            ELSE IF(E.LT.1303.4) THEN
            X=-81.32915 + (62.2775 + 3.)*Elog
     &         - 15.00826*Elog*Elog +
     &         1.558686*(Elog**3) - 6.1339621E-2*(Elog**4)
 
                ELSE
                X=-9.161526 + (10.07448 + 3.)*Elog
     &            - 1.435878*Elog*Elog + 5.2728362E-2*(Elog**3)
 
      ENDIF
      X=X-3.*Elog
      MAGNES1 = 0.
      IF(abs(X).lt.70.) MAGNES1 = Exp(X)
 
 
      RETURN
      END
C
C------------------------------------------------------------------------
C----------------------------------------------------------------------
C  REAL FUNCTION: NEON
C  Source: Atomic and Nuclear Data Tables, Jan. 1982
C
C  Description:  NEON calculates the mass absorption coefficient
C      for NEON1 gas.
C
C  Usage:  REAL FUNCTION NEON(E)
C      E = Energy in eV
C      NEON returns the mass absorption cross section in cm**2/g.
C
C  COMMON BLOCKS:
C      NONE
C
C  IMPLICIT
C      NONE
C
C  SUBROUTINES CALLED BY NEON:
C      NONE
C
C------------------------------------------------------------------------------
      FUNCTION NEON1(E)
C
C
      REAL E,ELOG,X,NEON1
C
C
        Elog = Alog(E)
        If(E.LT.867.) Then
            X = -3.04041 + (13.0071*Elog) + (-1.93205*Elog*Elog) +
     +          (0.0977639*Elog*Elog*Elog)
        Else
            X = 17.6007 + (3.29278*Elog) + (-0.263065*Elog*Elog) +
     +          (5.68290E-3*Elog*Elog*Elog)
        Endif
      X=X-3.*Elog
      NEON1 = 0.
      IF(abs(X).lt.70.) NEON1 = Exp(X)
       RETURN
      END
C
C---------------------------------------------------------------------------
C---------------------------------------------------------------------------
C
C  REAL FUNCTION: NICKEL
C  Source: Atomic & Nuclear Data Tables, January 1982
C
C  Description: NICKEL calculates the mass absorbtion coefficient for nickel.
C  History: updated below L-edge (853.6 eV) - Aug 30, 1991 (MBC)
C
C  Usage: REAL FUNCTION NICKEL(E)
C      E = Energy in eV
C      NICKEL returns the mass absorption cross section in cm**2/g.
C
C  COMMON BLOCKS:
C      NONE
C  IMPLICIT
C      NONE
C
C  SUBROUTINES CALLED BY NICKEL:
C      NONE
C
C
C---------------------------------------------------------------------------
      FUNCTION NICKEL1(E)
 
      REAL E, ELOG, X, NICKEL1
 
        Elog = Alog(E)
 
        If(E.LT.853.6) Then
        X = -7.919931 + (11.06475 + 3.)*Elog - 1.935318*Elog*Elog
     &      + 9.3929626e-2*(Elog**3)
 
            Else if(E.LT.8331.6) Then
            X = 3.71129 + (8.45098*Elog) + (-0.896656*Elog*Elog) +
     +      (0.0324889*Elog*Elog*Elog)
                Else
                X = 28.4989 + (0.485797*Elog)
        Endif
      X=X-3.*Elog
      NICKEL1 = 0.
      IF(abs(X).lt.70.) NICKEL1 = Exp(X) 
 
      Return
      End
C--------------------------------------------------------------------------
C--------------------------------------------------------------------------
C--------------------------------------------------------------------------
C
C  REAL FUNCTION: NITRO
C  Source: Atomic & Nuclear Data Tables, January 1982
C
C  Description:  NITRO calculates the mass absorption cross-section of
C        nitrogen as a function of energy.
C
C  Usage: REAL FUNCTION NITRO(E)
C      E = Energy in eV
C      NITRO returns the mass absorption cross-section in cm**2/g.
C
C  COMMON BLOCKS:
C      NONE
C
C  IMPLICIT
C      NONE
C
C  SUBROUTINES CALLED BY NITRO:
C      NONE
C
C--------------------------------------------------------------------------
      FUNCTION NITRO1(E)
C
C
      REAL E,ELOG,X,NITRO1
C
C
        Elog = Alog(E)
        If(E.LT.401.) Then
            X = 9.24058 + (7.02985*Elog) + (-1.08849*Elog*Elog) +
     +      (0.0611007*Elog*Elog*Elog)
        Else
            X = -13.0353 + (15.4851*Elog) + (-1.89502*Elog*Elog) +
     +          (0.0769412*Elog*Elog*Elog)
        Endif
      X=X-3.*Elog
      NITRO1 = 0.
      IF(abs(X).lt.70.) NITRO1 = Exp(X)
      RETURN
      END
C
C------------------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C  REAL FUNCTION: OXYGEN
C
C  Source:        Atomic & Nuclear Data Tables, January 1982
C
C  Description:  OXY Calculates the mass absorption cross-section of oxygen
C      as a function of energy.
C  History: updated above K-edge (531.7 eV) - Aug 30, 1991 (MBC)
C
C  Usage: REAL FUNCTION OXYGEN(E)
C      E = Energy in eV
C      OXYGEN returns the mass absorption cross-section in cm**2/g.
C
C  COMMON BLOCKS:
C      NONE
C
C  IMPLICIT
C      NONE
C  SUBROUTINES CALLED BY OXYGEN:
C      NONE
C
C------------------------------------------------------------------------
      FUNCTION OXYGEN1(E)
 
      REAL E, Elog, X, OXYGEN1
 
        Elog = Alog(E)
        If(E.LT.531.7) Then
            X = 2.57264 + (10.9321*Elog) + (-1.79383*Elog*Elog) +
     +          (0.102619*Elog*Elog*Elog)
        Else
            X = 16.53869 + (0.6428144 + 3.)*Elog - 0.3177744*Elog*Elog
     &          + 7.9471897e-3*(Elog**3)
 
        Endif
      X=X-3.*Elog
      OXYGEN1 = 0.
      IF(abs(X).lt.70.) OXYGEN1 = Exp(X) 
 
      RETURN
      END
C--------------------------------------------------------------------------
C--------------------------------------------------------------------------
C       FUNCTION SILICN
C
C       Source: Atomic and Nuclear Data Tables, January 1982
C
C       Description: SILICN calculates the mass absorption cross section
C                for silicon in cm**2/g.
C       History: updated Aug 30, 1991 (MBC)
C                updated March 4, 1992 (MBC)
C
C       COMMON BLOCKS:
C               NONE
C
C       IMPLICIT
C               none
C
C       SUBROUTINES CALLED:
C               NONE
C
C---------------------------------------------------------------------------
        FUNCTION SILICN1(E)
 
          REAL E,Elog,X,SILICN1
C
        Elog = Alog(E)
 
        If(E.LT.100.6) Then
            X = -3.066295 + (7.006248 + 3.)*Elog - 0.9627411*Elog*Elog
        Else if(E.LT.1840.0) Then
                X = -182.7217 + (125.061 + 3.)*Elog - 29.47269*Elog*Elog
     &              + 3.03284*(Elog**3) - 0.1173096*(Elog**4)
            Else
                X = -33.39074 + (18.42992 + 3.)*Elog
     &              - 2.385117*Elog*Elog + 8.887583e-2*(Elog**3)
 
        Endif

      X=X-3.*Elog
      SILICN1 = 0.
      IF(abs(X).lt.70.) SILICN1 = Exp(X) 
        Return
        End
C-------------------------------------------------------------------------
C------------------------------------------------------------------------
C     Real Function: SODIUM
C     Source: Atomic & Nuclear Data Tables, January 1982
C
C     Description:
C     Calculates mass absorption coefficient (mu/rho) for sodium
C
C     History: original Aug 6, 1991 (MBC)
C
C     Usage:   FUNCTION SODIUM(E)
C              E = Energy in eV
C              SODIUM returns mu/rho in cm**2/g
C
C     COMMON Blocks:
C            none
C
C     IMPLICIT:
C            none
C-------------------------------------------------------------------------
C
      FUNCTION SODIUM1(E)
C
C     IMPLICIT NONE
C
      REAL E,ELOG,X,SODIUM1
C
C
      Elog = Alog(E)
 
      IF(E.LT.1071.7) THEN
      X= -2737.598 + (2798.704 + 3.)*Elog - 1009.892*Elog*Elog +
     &  87.16455*(Elog**3) + 43.20644*(Elog**4) - 15.27259*(Elog**5)
     & + 2.180531*(Elog**6) - 0.1526546*(Elog**7) +
     &   4.3137977E-3*(Elog**8)
 
            ELSE
            X=1.534019 + (6.261744 + 3.)*Elog
     &        - 0.9914126*Elog*Elog +
     &        3.5278253E-2*(Elog**3)
 
 
      ENDIF
      X=X-3.*Elog 
      SODIUM1 = 0.
      IF(abs(X).lt.70.) SODIUM1 = Exp(X)
 
 
      RETURN
      END
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C     Real Function: SULFUR
C     Source: Atomic & Nuclear Data Tables, January 1982
C
C     Description:
C     Calculates mass absorption coefficient (mu/rho) for sulfur
C
C     History: original Aug 6, 1991 (MBC)
C
C     Usage:   FUNCTION SULFUR(E)
C              E = Energy in eV
C              SULFUR returns mu/rho in cm**2/g
C
C     COMMON Blocks:
C            none
C
C     IMPLICIT
C            none
C-------------------------------------------------------------------------
C
      FUNCTION SULFUR1(E)
C
C
      REAL E,ELOG,X,SULFUR1
C
C
      Elog = Alog(E)
 
 
      IF(E.LT.165.0) THEN
      X=598.2911 + (3. - 678.2265)*Elog + 308.1133*Elog*Elog
     &  - 68.99324*(Elog**3) + 7.62458*(Elog**4)
     &  - 0.3335031*(Elog**5)
 
            ELSE IF(E.LT.2470.5) THEN
            X=3994.831 + (3. - 3693.886)*Elog +
     &        1417.287*Elog*Elog - 287.9909*(Elog**3) +
     &        32.70061*(Elog**4) - 1.968987*(Elog**5) +
     &        4.9149349E-2*(Elog**6)
 
                ELSE
                X=-22.49628 + (14.24599+ 3.)*Elog -
     &             1.848444*Elog*Elog + 6.6506132E-2*(Elog**3)
 
      ENDIF
      X=X-3.*Elog
      SULFUR1 = 0.
      IF(abs(X).lt.70.) SULFUR1 = Exp(X)
 
 
      RETURN
      END
C
C--------------------------------------------------------------------------
C--------------------------------------------------------------------------
