program intrat3D

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
c      
c    The  program  is  a  fortran  77  code provided  as  an  accessory  to
c    interrogate  the  provided datasets,  which  consist  of emission  and
c    recombination  coefficients of  hydrogen  calculated for  Case B  with
c    electron  energies described  by  kappa distributions.   The data  are
c    tabulated  on  a  grid   of  electron  number  density,  Ne,  electron
c    temperature, Te, and kappa.   If necessary the program interpolates to
c    the  user  defined  values  of  the three  parameters  using  Lagrange
c    interpolation.  Interpolated  values  are  given  for  polynomials  of
c    orders: 1, 3 and 5. The  purpose of multiple interpolations is to test
c    the convergence and hence check the reliability of the interpolation.
c    
c    Two main input data files are  required to run the program: e1bk.d and
c    t1bk.d.  The first file contains emission coefficients for transitions
c    between states with  upper quantum number nu and  lower quantum number
c    nl,  where the quantum  numbers range  between 2-99  with nl<nu,  as a
c    function of Ne, Te and  kappa. The second file contains total hydrogen
c    recombination  coefficients and recombination  coefficients to  the 2s
c    state as a  function of Ne, Te and kappa. The  structure of these data
c    files  is fully  explained in  the  documentation of  the dataset  and
c    associated paper (see ReadMe.dat file  and [1] in the References; also
c    refer to [2]).
c    
c    During the  execution of the program,  options are given  to guide the
c    user and obtain the required values.  All input and output data are in
c    the cgs system.  There are three main options:
c    
c    1)  Obtaining the emission coefficient in  units  of   erg.cm^3.s^{-1} 
c    for a given transition. The required inputs are the  upper  and  lower  
c    quantum numbers  of  the transition, as well as Ne, Te and kappa.
c    
c    2) Obtaining the total recombination  coefficient of H for a given set
c    of physical parameters. The required inputs are Ne, Te and kappa.
c    
c    3) Obtaining the total recombination coefficient to the 2s state of  H
c    for a given set of physical parameters. The required inputs are Ne, Te 
c    and kappa.
c    
c    As indicated already, the raw  and processed data are calculated based
c    on the  assumption of kappa  electron energy distribution  rather than
c    Maxwell-Boltzmann   (MB)   distribution.    However,  data   with   MB
c    distribution can also be obtained  by using large values of kappa. For
c    this purpose, the  use of the maximum available  kappa value, which is
c    10^6, is recommended.
c    
c    ====================================================
c    
c    Program data sheet:
c    -Name: intrat3D
c    -Version: 1.0
c    -Date: 2014
c    -Language: fortran 77
c    -Tested compilers: gfortran, f77, Intel fortran
c    -Compilation from command line (e.g. for f77): $f77 intrat3D.f
c    -Tested platforms: Ubuntu 12.04, Scientific Linux
c    
c    ====================================================
c    
c    Acknowledgement:
c    
c    -P.J. Storey: University College London, pjs@star.ucl.ac.uk
c    -Taha Sochi:  University College London, t.sochi@ucl.ac.uk
c    
c    ====================================================
c    
c    Notes:
c    
c    (1) The  intrat code and documentation  of Hummer and  Storey [2] were
c    consulted  and  used  as  a   model  for  writing  intrat3D  code  and
c    documentation.
c    
c    (2)  Error  reports,  suggestions  and  enquiries should  be  sent  to
c    P.J. Storey or Taha Sochi.
c    
c    ====================================================
c    
c    References:
c    
c    [1]  P.J.  Storey,  Taha   Sochi  (2014)  Emission  and  recombination
c    coefficients for  hydrogen with kappa  electron  energy  distribution.
c    MNRAS X(X): X-X.
c    
c    [2] P.J.   Storey, D.G.  Hummer (1995)  Recombination line intensities
c    for  hydrogenic  ions  -  IV.  Total  recombination  coefficients  and
c    machine-readable tables for Z=1 to 8. MNRAS 272(1): 41-48.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      write(*,5)

      call hdatk

      call getinput

      stop

5     format(//'           ****** Welcome to Intrat3D ******'//,
     * ' This code is provided as an accessory to the data associated'/
     * ' with the following paper:'/
     * ' P.J. Storey, Taha Sochi. Emission and recombination'/
     * ' coefficients for hydrogen with kappa electron energy'/
     * ' distribution. MNRAS XXX.'// 
     * ' Error reports, suggestions and enquiries should be sent to'/
     * ' P.J. Storey (UCL, pjs@star.ucl.ac.uk) or'/
     * ' Taha Sochi  (UCL, t.sochi@ucl.ac.uk).'//
     * ' Interpolation of emission and recombination coefficinets'/
     * ' for hydrogen calculated in Case B with kappa-distributed'/
     * ' electron energies. Maxwell-Boltzmann results can be'/
     * ' regained by specifying a large value of kappa (e.g. 10^6).'//     
     * ' Results are given for first, third and fifth order'/
     * ' Lagrange interpolation in the grid of computed'/
     * ' coefficients so that convergence can be monitored.'//
     * ' Results are a function of electron number density,'/
     * ' electron temperature and kappa'/)

      end

*******************************************************************************

      subroutine getinput

      implicit none

      integer mxl,mxt,mxn,mxk
      parameter (mxl=5000,mxn=20,mxt=20,mxk=50)

      character char*1,opts*1
      integer nu,nl,il,i,maxx,np,iopt,ni(3),istop
      real*8 xn,yt,zk,ri(3),ThreeDntk

      real*8 densx,tempx,kappax,ex,a,a2s,skap
      integer ntempx,ndensx,nkappax,nmin,ncut

      common/theory/densx,tempx,kappax,skap,ex,a,a2s,ntempx,ndensx,
     *              nkappax,nmin,ncut

      dimension densx(mxn),tempx(mxt),kappax(mxk),ex(mxl,mxn,mxt,mxk),
     *     a(mxn,mxt,mxk),a2s(mxn,mxt,mxk),skap(mxk)
      dimension opts(5)

      data ni/1,2,3/, maxx/3/   ! interpolation parameters
      data opts/'e','i','n','c','r'/

      open(unit=16,file='intrat3D.d',status='unknown')

      write(*,1) 
 2    write(*,*) 'Enter iopt'
      read(*,*) iopt 

      nu=0
      nl=0
 3    if(iopt.eq.1) then
         write(*,50) ncut 
         read(*,*) nu,nl
         if(nu.le.nl) then
           write(*,51)
           go to 3
         endif    
         if(nu.gt.ncut.or.nl.gt.ncut.or.nu.lt.3.or.nl.lt.2.and.
     *         nl.lt.nu) then
           write(*,52) ncut 
           go to 3
         endif 
       elseif(iopt.ne.2.and.iopt.ne.3) then
         write(*,*) 'Invalid value of iopt'
         go to 2
      endif

 4    write(*,*) 'Enter electron number density, temperature and kappa'
      read(*,*) xn,yt,zk
         
      istop=0
      if(xn.lt.10.d0**densx(1).or.xn.gt.10.d0**densx(ndensx)) then
         write(*,5) 10.d0**densx(1),10.d0**densx(ndensx)
         istop=1
      endif
      if(yt.lt.10.d0**tempx(1).or.yt.gt.10.d0**tempx(ntempx)) then 
         write(*,6) 10.d0**tempx(1),10.d0**tempx(ntempx)
         istop=1
      endif
      if(zk.lt.10.d0**kappax(1).or.zk.gt.10.d0**kappax(nkappax)) 
     *     then 
         write(*,7) 10.d0**kappax(1),10.d0**kappax(nkappax)
         istop=1
      endif
      if(istop.eq.1) go to 8 

      il=(((ncut-nu)*(ncut+nu-1))/2)+nl

      do i=1,maxx
         np = ni(i)
         ri(i) = ThreeDntk(il,xn,yt,zk,np,iopt)
      enddo
      
      if(iopt.eq.1) then
      
        write(*,200)
        write(16,200)
        write(*,201) (2*ni(i)-1,i=1,maxx)
        write(16,201) (2*ni(i)-1,i=1,maxx) 
        write(*,202) nu,nl,xn,yt,zk,(ri(i),i=1,maxx)
        write(16,202) nu,nl,xn,yt,zk,(ri(i),i=1,maxx)
      elseif(iopt.eq.2) then
        write(*,203)
        write(16,203)
        write(*,201) (2*ni(i)-1,i=1,maxx)
        write(16,201) (2*ni(i)-1,i=1,maxx) 
        write(*,204) xn,yt,zk,(ri(i),i=1,maxx)
        write(16,204) xn,yt,zk,(ri(i),i=1,maxx)
      elseif(iopt.eq.3) then
        write(*,205)
        write(16,205)
        write(*,201) (2*ni(i)-1,i=1,maxx)
        write(16,201) (2*ni(i)-1,i=1,maxx) 
        write(*,204) xn,yt,zk,(ri(i),i=1,maxx)
        write(16,204) xn,yt,zk,(ri(i),i=1,maxx)
      endif
 8    if(iopt.eq.1) then
        write(*,*) 'exit(e), new iopt(i), new transition(n),', 
     *  ' new physical conditions(c)'
      elseif(iopt.eq.2.or.iopt.eq.3) then
        write(*,*) 'exit(e), new iopt(i), new physical conditions(c)'
      endif
      read(*,*) char
      if(char.eq.opts(1)) stop
      if(char.eq.opts(2)) go to 2
      if(char.eq.opts(3)) go to 3
      if(char.eq.opts(4)) go to 4
      write(*,*) 'not recognised'
      go to 8

      return

 1    format(1x,/' Input and units (in cgs):'/
     *   ' Electron number density: cm^{-3}'/
     *   ' Electron temperature: K'/
     *   ' Kappa: no units'//
     *   ' Input options (iopt)'/
     *   ' iopt = 1 for emission coefficient for a specific transition'/
     *   ' iopt = 2 for total recombination coefficient for H in'
     *   ' Case B'/
     *   ' iopt = 3 for recombination coefficient to the H 2s'
     *   ' state'//
     *   ' Output and units (in cgs):'/
     *   ' iopt = 1  emission coefficient in erg.cm^3.s^{-1}'/
     *   ' iopt = 2 and 3 recombination coefficients in cm^3.s^{-1}'//)
 5    format(' Density outside allowable range ',1pd10.2,
     *   ' to ',1pd10.2)
 6    format(/' Temperature outside allowable range ',1pd10.2,
     *   ' to ',1pd10.2)
 7    format(/' kappa outside allowable range ',1pd10.2,
     *   ' to ',1pd10.2)
 50   format('Enter upper and lower n, between ',i2,' and 2')
 51   format(' upper n may not be less or equal lower n')
 52   format(' upper n or lower n outside allowed range',
     *       ' (',i2,',2)')
 200  format(//53x,'Emission coefficient')
 201  format(10x,'density',5x,'temperature'1x,'kappa',4x,
     *    'interp. order = ',i3,2i11)
 202  format(1x,2i3,1p3d12.3,8x,1p3d11.3//)
 203  format(//53x,'Recombination coefficient (total)')
 204  format(1x,6x,1p3d12.3,8x,1p3d11.3//)
 205  format(//53x,'Recombination coefficient (2s)')

      end
      
*******************************************************************************
      DOUBLE PRECISION FUNCTION PHI(XIQ,LG,ITAU,X)
      IMPLICIT NONE
      REAL*8 XIQ,FK,X
      INTEGER LG,ITAU,L
      DIMENSION XIQ(10)
      FK=1.
      DO L=1,LG
        IF(L.NE.ITAU) THEN
          FK=FK*(X-XIQ(L))
        ENDIF
      ENDDO
      PHI=FK
      RETURN
      END

***********************************************************************

      SUBROUTINE LAGR(F,X,FV,XV,NP,NI,NDIM,IER)
      IMPLICIT NONE
      REAL*8 F,X,FL,FV,XV,XIQ,PHI
      INTEGER NP,NI,NDIM,IER,LG,IN,ILOW,ITAU,I
C
C  SUBROUTINE TO PERFORM LAGRANGIAN INTERPOLATION OF A SET OF
C  TABULATED FUNCTION VALUES F(X(I))
C
C  INPUTS
C  F(I)  = ARRAY OF TABULATED FUNCTION VALUES
C  NI    = NUMBER OF TABULATED VALUES
C  2*NP  = NUMBER OF POINTS TO BE USED IN INTERPOLATION
C  XV    = VALUE OF X AT WHICH FUNCTION VALUE IS REQUIRED
C
C  OUTPUTS
C  FV    = INTERPOLATED VALUE
C  IER   = ERROR CODE (.NE.0 IF XV OUTSIDE RANGE OF X)
C
      DIMENSION F(NDIM),X(NDIM),XIQ(10)

*      write(*,*) 'in LAGR: np, xv',np,xv

      LG=2*NP
      IER=0
      IF (XV.LT.X(1).OR.XV.GT.X(NI)) THEN
        IER=1
        RETURN
      ENDIF
C
C  LOCATE XV IN X
      DO I=2,NI
        IF(X(I).GE.XV) THEN
C
C  SELECT SUBSET OF X FOR INTERPOLATION, TAKING CARE AT ENDS
          IN=NI-LG+1
          ILOW=I-NP
          IF(ILOW.LT.1) ILOW=1
          IF(ILOW.GT.IN) ILOW=IN
          DO ITAU=1,LG
            XIQ(ITAU)=X(ILOW+ITAU-1)
          ENDDO

C  CONSTRUCT INTERPOLATION COEFFICIENTS

          FV=0.
          DO ITAU=1,LG
            FL=PHI(XIQ,LG,ITAU,XV)/PHI(XIQ,LG,ITAU,XIQ(ITAU))
            FV=FV+FL*F(ILOW+ITAU-1)
          ENDDO
          RETURN
        ENDIF
      ENDDO

      END

c***********************************************************************

      subroutine hdatk

c***********************************************************************
c
c     Read data from file e1bk.d and t1bk.d and store them in the
c     common block hdatak.
c
c***********************************************************************

      implicit none

      integer mxl,mxt,mxn,mxk
      parameter (mxl=5000,mxn=20,mxt=20,mxk=50)

      real*8 densx,tempx,kappax,ex,a,a2s,skap
      integer ntempx,ndensx,nkappax,nmin,ncut

      real*8 zk,gammln
      integer ia,ib,j,i,ne,ik,k,maxx,ni

      character*3 name

      common/theory/densx,tempx,kappax,skap,ex,a,a2s,ntempx,ndensx,
     *              nkappax,nmin,ncut

      dimension densx(mxn),tempx(mxt),kappax(mxk),ex(mxl,mxn,mxt,mxk),
     *     a(mxn,mxt,mxk),a2s(mxn,mxt,mxk),skap(mxk)

      dimension ni(3)

      data maxx/3/
      data ni/1,2,3/   ! interpolation parameters

      name = '1bk'

      write(*,*) 'reading dataset e' // name // '.d'
c
      open(unit=15,file='e'//name//'.d',status='old')
      open(unit=16,file='t'//name//'.d',status='old')
c
c read e1bk.d dataset
c
      read(15,*) ndensx,ntempx,nkappax
      do 102 ik=1,nkappax
         do 101 ia=1,ntempx
           do 100 ib=1,ndensx
                read(15,25) densx(ib),tempx(ia),kappax(ik),
     *                      nmin,ncut
*                kappax(ik)=log10(abs(kappax(ik)))
*                tempx(ia)=log10(tempx(ia))
*                densx(ib)=log10(densx(ib))
                 ne=ncut*(ncut-1)/2
                read(15,30) (ex(j,ib,ia,ik),j=1,ne)
 100       continue
 101    continue
 102  continue

      write(*,*) 'read completed'

c generate kappa scale factors

      do ik=1,nkappax
         zk=10.d0**kappax(ik)
         skap(ik)=1.5d0*dlog(zk-1.5d0)+gammln(zk-0.5d0)-
     *                       gammln(zk+1.d0)
         skap(ik)=exp(skap(ik))/(1.d0-1.5d0/zk)
      enddo

c read t1bk.d dataset

      write(*,*) 'reading dataset t' // name // '.d'
      read(16,*) (((a(i,k,j),i=1,ndensx),k=1,ntempx),j=1,nkappax)
      read(16,*) (((a2s(i,k,j),i=1,ndensx),k=1,ntempx),j=1,nkappax)
      write(*,*) 'read completed'

      write(*,'(/A,i4 )') ' Minimum principal quantum number =   2'
      write(*,'(A,i4 )') ' Maximum principal quantum number =',ncut

      write(*,32) 10.d0**densx(1),10.d0**densx(ndensx)

      write(*,33) 10.d0**tempx(1),10.d0**tempx(ntempx)

      write(*,34) 10.d0**kappax(1),10.d0**kappax(nkappax)

      return

 25   format(6x,f10.3,f10.3,4x,f10.3,5x,2i5)
 30   format((8e10.3))
 32   format(/' Range of densities:',3x,1pe10.3,' - ',e10.3)
 33   format(' Range of temperatures:',1pe10.3,' - ',e10.3)
 34   format(' Range of kappa:',7x,1pe10.3,' - ',e10.3)


      end


********************************************************************************************

      double precision function ThreeDntk(il,xn,yt,zk,np,iopt)

      implicit none

      real*8 xn,yt,zk,OneDn,gg,gammln
      integer il,np,iopt

      ThreeDntk=10.d0**OneDn(il,xn,yt,zk,np,iopt)

* remove kappa scaling

      gg=1.5d0*dlog(zk-1.5d0)+gammln(zk-0.5d0)
     :                          -gammln(zk+1.d0)
      gg=exp(gg)/(1.d0-3.d0/(2.d0*zk))
      ThreeDntk=ThreeDntk/gg

      return
      end

******************************************************************************

      double precision function OneDn(il,xn,yt,zk,np,iopt)

      implicit none

      integer mxl,mxn,mxt,mxk
      parameter (mxl=5000,mxn=20,mxt=20,mxk=50)

      real*8 densx,tempx,kappax,ex,a,a2s,skap
      integer ntempx,ndensx,nkappax,nmin,ncut

      real*8 fn,an,yt,zk,fvn,xn,xnl,OneDt
      integer il,nden,in,np,ier,iopt

      common/theory/densx,tempx,kappax,skap,ex,a,a2s,ntempx,ndensx,
     *              nkappax,nmin,ncut

      dimension densx(mxn),tempx(mxt),kappax(mxk),ex(mxl,mxn,mxt,mxk),
     *     a(mxn,mxt,mxk),a2s(mxn,mxt,mxk),skap(mxk)
      dimension fn(mxn),an(mxn)

      do in=1,ndensx
         fn(in)=OneDt(il,in,yt,zk,np,iopt)
         an(in)=densx(in)
      enddo

      xnl=log10(xn)
      call lagr(fn,an,fvn,xnl,np,ndensx,mxn,ier)
      OneDn=fvn

      return

      end

**************************************************************************

      double precision function OneDt(il,in,yt,zk,np,iopt)

      implicit none 

      integer mxl,mxn,mxt,mxk
      parameter (mxl=5000,mxn=20,mxt=20,mxk=50)

      real*8 densx,tempx,kappax,ex,a,a2s,skap
      integer ntempx,ndensx,nkappax,nmin,ncut

      real*8 xt,ft,yt,zk,ytl,fvt,OneDk
      integer il,in,np,ier,it,iopt

      common/theory/densx,tempx,kappax,skap,ex,a,a2s,ntempx,ndensx,
     *              nkappax,nmin,ncut

      dimension densx(mxn),tempx(mxt),kappax(mxk),ex(mxl,mxn,mxt,mxk),
     *     a(mxn,mxt,mxk),a2s(mxn,mxt,mxk),skap(mxk)

      dimension ft(mxt),xt(mxt)

      do it=1,ntempx
         xt(it)=tempx(it)
         ft(it)=OneDk(il,in,it,zk,np,iopt)
      enddo

      ytl=log10(yt)
      call lagr(ft,xt,fvt,ytl,np,ntempx,mxt,ier)
      OneDt=fvt

      return

      end

***********************************************************************

      double precision function OneDk(il,in,jt,zk,np,iopt)

      implicit none

      integer mxl,mxn,mxt,mxk
      parameter (mxl=5000,mxn=20,mxt=20,mxk=50)

      real*8 densx,tempx,kappax,ex,a,a2s,skap
      integer ntempx,ndensx,nkappax,nmin,ncut

      real*8 xk,fk,zk,zkl,fvk
      integer il,in,jt,ik,ier,np,iopt

      common/theory/densx,tempx,kappax,skap,ex,a,a2s,ntempx,ndensx,
     *              nkappax,nmin,ncut

      dimension densx(mxn),tempx(mxt),kappax(mxk),ex(mxl,mxn,mxt,mxk),
     *     a(mxn,mxt,mxk),a2s(mxn,mxt,mxk),skap(mxk)

      dimension fk(mxk),xk(mxk)

      do ik=1,nkappax

         xk(ik)=kappax(ik)
         if(iopt.eq.1) then
           fk(ik)=ex(il,in,jt,ik)*skap(ik)
         elseif(iopt.eq.2) then
           fk(ik)=a(in,jt,ik)*skap(ik)
         elseif(iopt.eq.3) then
           fk(ik)=a2s(in,jt,ik)*skap(ik)
         endif
         fk(ik)=log10(fk(ik))
      enddo

      zkl=log10(zk)
      call lagr(fk,xk,fvk,zkl,np,nkappax,mxk,ier)
      OneDk=fvk

      return

      end

***********************************************************************************

      FUNCTION gammln(xx)
c
      implicit none
c
      INTEGER j
      DOUBLE PRECISION gammln,xx,ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     *-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END