      implicit none

      integer iopt,iopt1
      integer istop
      real*8 wmin,wmax

* read theoretical emissivities from individual files or binary file "OIIlines_AB" and pick
* out the relevant lines


      istop=0
      if(istop.eq.1) stop

 1    write(6,10)
 10   format(/1x,'Choose one of three options:' //
     : '1: Generate a line list' // 
     : '2: Extract emission coefficients for a user supplied list of'
     : ' lines in user'/'   specified conditions' //
     : '3: Exit' //
     : 1x,'Enter an integer:')

      read(5,*) iopt

      if(iopt.ne.1.and.iopt.ne.2.and.iopt.ne.3) go to 1

      if(iopt.eq.3) stop

      call reademiss

      if(iopt.eq.1) then
        
 2      write(6,20)
 20     format(//1x,'Further options:'//'1: List all lines'
     :  ' (8889 between 367.97 and 3.90713e+07 Angstrom)'//
     :   '2: Specify a wavelength range for the'
     :  ' list'//1x,'Enter an integer')

        read(5,*) iopt1

        if(iopt1.ne.1.and.iopt1.ne.2) go to 2

        wmin=367.97d0
        wmax= 3.907135E+07

        if(iopt1.eq.2) then
 3         write(6,30) 
 30       format(/1x,'Enter lower and upper wavelength in Angstrom')
          read(5,*) wmin, wmax
          if(wmin.gt.wmax) then
            write(6,31) 
 31         format(1x,'Lower wavelength greater than upper')
            go to 3
          endif 
        endif

      endif

      if(iopt.eq.1) then

        call opt1(wmin,wmax)

      elseif(iopt.eq.2) then

        call opt2

      endif

      go to 1

      stop

      end

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

      subroutine opt1(wmin,wmax)

      implicit none
    
* mxl = maximum number of line emission coeffs
* mxt = maximum number of temperatures
* mxd = maximum number of densities
* mxout = maximum number of user requested emission coeffs

      integer mxl,mxt,mxd,mxout

      parameter (mxl=10000,mxt=25,mxd=16,mxout=30)
      integer iopt,istop

      character*1 eflag(mxl),eflagc(mxl),eflagcs,ul,Case
      character*84 txt1(mxl),txt1c(mxl),txt1cs
      character*32 txt2(mxl),txt2c(mxl),txt2cs
      character*19 woutfile,eoutfile

      real*8 em(mxt,mxd,mxl),wl(mxl),
     :       eu(mxl),el(mxl),temlg(mxt),
     :       denslg(mxd),wlc(mxl),emc(mxl),
     :       emcs,wmin,wmax,ws,es,emmax,wv1,ev1

      integer nt,nd,il,inl,jl,nl,nw

      common/alltheory/wl,em,eflag,eu,el,txt1,txt2,temlg,denslg,
     :                 nt,nd,nw,Case

      data ul/'_'/
      data nl/8889/

      woutfile='OIIdata'//ul//'list'//ul//'worder'
      eoutfile='OIIdata'//ul//'list'//ul//'eorder'

      open(unit=16,status='unknown',file=woutfile)
      open(unit=17,status='unknown',file=eoutfile)

* Load line list for canonical temperature and density into local storage
* and select those in specified range

      inl=0

      emmax=0.d0
      wv1=4649.13d0

      do il=1,nl
         if(abs(wl(il)-wv1).lt.1.d-2) then
           ev1=em(21,11,il)
         endif
         if(wl(il).ge.wmin.and.wl(il).le.wmax) then
           inl=inl+1
           wlc(inl)=wl(il)
           emc(inl)=em(21,11,il)
           txt1c(inl)=txt1(il)
           txt2c(inl)=txt2(il)
           eflagc(inl)=eflag(il)
           if(emc(inl).gt.emmax) emmax=emc(inl)
c           write(6,*) inl,il,wlc(inl)
         endif
      enddo

      write(16,1) inl,wmin,wmax
      write(17,1) inl,wmin,wmax
      write(6,1) inl,wmin,wmax
 1    format(/1x,'Found ',i5,' lines in interval ',f12.2,
     :    ' to ',f12.2,' Angstrom'/)
      if(inl.eq.0) return
      write(6,11) eoutfile,woutfile
c      write(16,11) eoutfile,woutfile
 11   format(1x,'File ',a19,
     :  ' contains the line list in order of decreasing emission'
     :  ' coefficient'//
     :  1x,'File ',a19,' contains the full list preceded by a list'
     :  ' of only the strongest lines, in order of increasing'
     :  ' wavelength')
 

* Sort into increasing wavelength order

      do il=1,inl-1
         ws=wlc(il)
         do jl=il+1,inl
           if(wlc(jl).lt.ws) then
             wlc(il)=wlc(jl)          
             wlc(jl)=ws
             ws=wlc(il)
             emcs=emc(il)
             txt1cs=txt1c(il)
             txt2cs=txt2c(il)
             emc(il)=emc(jl)
             emc(jl)=emcs
             txt1c(il)=txt1c(jl)
             txt1c(jl)=txt1cs
             txt2c(il)=txt2c(jl)
             txt2c(jl)=txt2cs
             eflagcs=eflagc(il)
             eflagc(il)=eflagc(jl)
             eflagc(jl)=eflagcs
           endif
         enddo
      enddo

      write(16,4) Case
 4    format(' This file is in two sections:'/
     :       ' Section 1 - Lines with emission coefficient at'
     :       ' least 0.1% of strongest line in this wavelength range,'/
     :   12x,' in order of increasing wavelength'/
     :       ' Section 2 - All lines in this wavelength range in order'
     :       ' of increasing wavelength'//
     :       ' All emission coefficients at a temperature of 10000K and'
     :          ' an electron number density of 10000/cm^3 in Case 'a1/)

      write(16,41)
 41   format(1x,'I          = Line index in full list for this'
     :     ' wavelength range (see second section below)')
      write(16,42)
 42   format(
     :  1x,'X/T        = Wavelength from only experimental energies (X)'
     :     ' or at least one theoretical energy (T)'/
     :  1x,'             T implies wavelength not' 
     :     ' spectroscopically useful'/
     :  1x,'W          = Wavelength[A],'/
     :  1x,'EM         = Emission coefficient[cgs],'/
     :  1x,'E/EMAX     = EM relative to the maximum value in this'
     :     ' wavelength range,'/
     :  1x,'E/EV1      = EM rel to V1 4649.13'/
     :  1x,'IPI, IPF   = parent level index'/
     :  1x,'NLI, NLF   = valence electron n, l of initial (I)'
     :     ' and final (F) states'/
     :  1x,'2JI, 2JF   = 2*J (2*total angular momentum)'/
     :  1x,'PI PF      = parity (0=EVEN, 1=ODD)'/)

      write(16,22) 
 22   format(' Section 1'/)
      write(16,21)

      do il=1,inl
 
        if((emc(il)/emmax).gt.1.d-3) then

          if(wlc(il).le.1.e5) then   

              write(16,5) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :                emc(il)/ev1,txt1c(il)
 5            format(1x,i5,f12.2,2x,a1,1x,1pe12.3,0p2f12.5,a84)
          else

            write(16,51) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :                emc(il)/ev1,txt1c(il)
 51         format(1x,i5,1pe12.5,2x,a1,1x,1pe12.3,0p2f12.5,a84)
          endif

        endif

      enddo

      write(16,2) 
 2    format(//' Section 2'/
     :         ' Full list in order of increasing wavelength'//)
      write(16,21) 
 21   format(5x,'I',11x,'W',14x,'EM',6x,'E/EMAX',7x,'E/EV1',3x,
     :       'Initial Term',12x,'Final Term',13x,
     :       ' IPI NLI 2JI PI',3x,'  IPF NLF 2JF PF')

      do il=1,inl

         if(wlc(il).le.1.e5) then
            if(emc(il)/ev1.gt.1.d-5) then 
              write(16,3) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :               emc(il)/ev1,txt1c(il)
 3            format(1x,i5,f12.2,2x,a1,1x,1p3e12.3,a84)
            else
              write(16,31) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :                emc(il)/ev1,txt1c(il)
 31           format(1x,i5,f12.2,2x,a1,1x,1p3e12.3,a84)
            endif
         else
            if(emc(il)/ev1.gt.1.d-5) then 
              write(16,32) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :               emc(il)/ev1,txt1c(il)
 32           format(1x,i5,1pe12.5,2x,a1,1x,1p3e12.3,a84)
            else
              write(16,33) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :                emc(il)/ev1,txt1c(il)
 33           format(1x,i5,1pe12.5,2x,a1,1x,1p3e12.3,a84)
            endif
          endif

      enddo


* Sort into decreasing emission coefficient order

      do il=1,inl-1
         es=emc(il)
         do jl=il+1,inl
           if(emc(jl).gt.es) then
             emc(il)=emc(jl)          
             emc(jl)=es
             es=emc(il)
             ws=wlc(il)
             wlc(il)=wlc(jl)
             wlc(jl)=ws
             txt1cs=txt1c(il)
             txt2cs=txt2c(il)
             txt1c(il)=txt1c(jl)
             txt1c(jl)=txt1cs
             txt2c(il)=txt2c(jl)
             txt2c(jl)=txt2cs
             eflagcs=eflagc(il)
             eflagc(il)=eflagc(jl)
             eflagc(jl)=eflagcs
           endif
         enddo
      enddo

      write(17,7) Case
 7    format(1x,'Lines in order of decreasing emission coefficient'/
     :       1x,'Emission coefficients at a temperature of 10000K and'
     :          ' an electron number density of 10000/cm^3 in Case 'a1/)
      write(17,43)
 43   format(1x,'I          = Line index')
      write(17,42)
      write(17,21)

c      do il=1,inl
 
c        if((emc(il)/emmax).gt.1.d-3) then  

c          if(wlc(il).le.1.e5) then   
c            write(17,5) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
c     :                emc(il)/ev1,txt1c(il)
c          else
c            write(17,51) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
c     :                emc(il)/ev1,txt1c(il)
c          endif

c        endif

c      enddo

c      write(17,6) 
c 6    format(//' Full list'//' In order of decreasing emissivity'//)
c      write(17,21)

      do il=1,inl

         if(wlc(il).le.1.e5) then
            if(emc(il)/ev1.gt.1.d-5.and.emc(il)/emmax.gt.1.d-5) then 
              write(17,5) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :               emc(il)/ev1,txt1c(il)
            else
              write(17,31) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :                emc(il)/ev1,txt1c(il)
            endif
         else
            if(emc(il)/ev1.gt.1.d-5.and.emc(il)/emmax.gt.1.d-5) then 
              write(17,51) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :               emc(il)/ev1,txt1c(il)
            else
              write(17,33) il,wlc(il),eflagc(il),emc(il),emc(il)/emmax,
     :                emc(il)/ev1,txt1c(il)
            endif
          endif

      enddo

      call flush(16)
      call flush(17)
      return

      end

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

      subroutine opt2

      implicit none

* mxl = maximum number of line emission coefficients
* mxt = maximum number of temperatures
* mxd = maximum number of densities
* mxout = maximum number of user requested emission coefficients

      integer mxl,mxt,mxd,mxout

      parameter (mxl=10000,mxt=25,mxd=16,mxout=30)

      character*84 txt1(mxl)
      character*32 txt2(mxl)
      character*30 infile
      character*34 outfile
      character*22 datname
      character*1 eflag(mxl),ef(mxl)
      character*1 ul,blk,Case

      real*8 dens(mxd),emms(mxt,mxd,mxout),wave(mxout),
     :  em(mxt,mxd,mxl),wl(mxl),wv(mxout),
     :  temlg(mxt),denslg(mxd),temlog(mxt),logdens(mxd),
     :  tin(mxt),din(mxd),eps2d(mxt,mxd,mxout),eu(mxl),
     :  el(mxl)
      real*8 xn,yt,TwoDnt,d,t,dt,dd,eps(3)

      integer nden,idens,nline,iline,k,ic,np,nt,nd,ntdin,itdin,ip,
     : ntem,item,iprint,i,n,nf(mxl),kc,isline,nw,
     : nsline,kline,ntin,ndin,itin,idin,status,access
      integer iopt,istop


      common/obs/wave,nline

      common/alltheory/wl,em,eflag,eu,el,txt1,txt2,temlg,denslg,
     :                 nt,nd,nw,Case

      common/theory/temlog,logdens,emms,ntem,nden

      data ul/'_'/, blk/' '/

      write(6,*) 'List of OII wavelengths is required in file',
     :             ' OIIdata_xxxx'
      write(6,*) 'Enter xxxx (1-22 characters)'
      read(5,*) datname
      infile='OIIdata'//ul//datname
      status=access(infile,'')

      if(status.ne.0) then
        write(6,*) 'File  ',infile,' does not exist'
        return
      endif

      k=index(infile,blk)
      infile=infile(1:k-1)
      outfile=infile(1:k-1)//ul//'out'

      open(unit=15,status='unknown',file=infile)
      open(unit=20,status='unknown',file=outfile)

      rewind(15)

* read 
* number of lines for which emission coeffs required - nline

      read(15,*) nline
      write(6,*) 'number of lines = ',nline

* read line wavelengths
* note that lines are identified by their wavelength as 
* printed in the lines.. files and given to 2 d.p.


      do iline=1,nline

        read(15,*) wave(iline)
        write(6,10) iline,wave(iline)
 10     format(1x'line number',i5,' wavelength [A] ',f12.2)
      enddo

      close(15)

      istop=0
      if(istop.eq.1) stop

      ntem=nt
      nden=nd
    
      do item=1,ntem
        temlog(item)=temlg(item)
      enddo
      do idens=1,nden
        logdens(idens)=denslg(idens)
      enddo

      do item=1,ntem
        do idens=1,nden

* mark all lines as not found

          do iline=1,nline
              nf(iline)=1
          enddo

          kline=0
          do iline=1,nline
            do isline=1,nw
               if(abs(wl(isline)-wave(iline))
     :                     .lt.0.001) then
                kline=kline+1 
                if(nf(iline).eq.0) then
                 write(6,*) 'found two lines that match ',
     :                       wave(iline)
c                 stop
                endif
                emms(item,idens,kline)=log10(em(item,idens,isline))
                ef(kline)=eflag(isline)
                wv(kline)=wave(iline)
                nf(iline)=0
               endif
            enddo
          enddo

* check that all lines were found

          do iline=1,nline
              if(nf(iline).gt.0) then
                write(6,*) 'line ',wave(iline),' not found'
              endif
          enddo

        enddo
      enddo


* and print for checking

*      do iline=1,nline
*          write(6,*) wave(iline)
*          do item=1,ntem
*            write(6,*) temlog(item),(emms(item,idens,iline),
*     :                idens=1,ndens)
*          enddo
*      enddo
      
      istop=0
      if(istop.eq.1) stop

 1    write(6,20) 
 20   format(/'Choose one of two options:' //
     :  '1: Provide a list of temperature/density pairs in a file'
     :  ' TDpairs_xxxx' //
     :  '2: Provide a range of temperatures and densities defining a'
     :  ' 1 or 2 dimensional grid in a file TDranges_xxxx'//
     :  'Enter an integer') 

      read(5,*) iopt

      if(iopt.ne.1.and.iopt.ne.2) go to 1

      write(6,*) 'Enter xxxx (1-21 characters)'

      read(5,*) datname

      write(6,21) outfile
 21   format(/1x,'Output will be in file  ',a34/)

      if(iopt.eq.1) then   
        infile='TDpairs'//ul//datname
      elseif(iopt.eq.2) then
        infile='TDranges'//ul//datname
      endif

      k=index(infile,blk)
      infile=infile(1:k-1)

      status=access(infile,'')
      if(status.ne.0) then
        write(6,*) 'File  ',infile,' does not exist'
        return
      endif

 30   format(1x,'Interpolation to log10(N[/cm^3]) = ',f7.3,
     :   '  log10(T[K]) = ',f7.3///5x,'line','wavelength[A]',
     :   ' Interp pts ',' coefficient [cm^3/s]')


      open(unit=18,status='unknown',file=infile)
      rewind(18)

      if(iopt.eq.1) then

        write(20,24)
 24     format(/1x,'I  =  Line index'/
     :          1x,'W  =  Wavelength[A],'/
     :          1x,'EM =  Emission coefficient[cgs],'/
     :          1x,'X  =  (air)wavelength from experimental energies'/
     :          1x,'T  =  one or both energies from theory,'
     :                   ' wavelength inaccurate')
        write(20,23)
 23     format(/1x,'Where necessary, the emission coefficients'
     :             ' are interpolated from'/
     :          1x,'the theory mesh using 6-point Lagrange in'
     :             ' temperature and density'/)  

        read(18,*) ntdin
        do itdin=1,ntdin
           read(18,*) t,d
           xn=log10(d)
           yt=log10(t)
c           write(6,25) t,d
           write(20,25) t,d,Case
 25        format(/1x,'Temperature = ',1pe12.3,'   Density = ',
     :          1pe12.3,'   Case = ',a1)
c           write(6,26) yt,xn
           write(20,26) yt,xn
           write(20,27)
 26        format(1x,'log10(T) =  ',f10.3,'       log10(N) = ',f7.3/)
 27        format(8x,'I',12x,'W',4x,'X/T',13x,'EM')
           if(xn.ge.denslg(1).and.xn.le.denslg(nden).and.
     :       yt.ge.temlg(1).and.yt.le.temlg(ntem)) then
             do iline=1,kline
               do ip=1,3
                 eps(ip)=10.d0**TwoDnt(iline,xn,yt,ip)
               enddo
c               write(6,40) iline,wv(iline),ef(iline),(eps(ip),ip=1,3)
               write(20,40) iline,wv(iline),ef(iline),eps(3)
             enddo
           else
             write(6,35)
             write(20,35)
 35          format(/1x,'Temperature or density outside available',
     :         ' theory range'/)
           endif
        enddo
 40     format(5x,i4,f13.2,5x,a2,1p3e15.3)

        elseif(iopt.eq.2) then

          write(20,23)

          read(18,*) ntin,ndin
ctest
          write(6,*) 'ntin,ndin',ntin,ndin

          if(ntin.eq.1) then
            read(18,*) tin(1)
            dt=0.d0
          elseif(ntin.gt.1) then
            read(18,*) tin(1),dt
ctest
            write(6,*) tin(1),dt

          endif
          if(ndin.eq.1) then
            read(18,*) din(1)
            dd=0.d0
          elseif(ndin.gt.1) then
            read(18,*) din(1),dd
          endif

          do itin=1,ntin
             tin(itin)=tin(1)+dt*(itin-1)
             yt=tin(itin)
             do idin=1,ndin
               din(idin)=din(1)+dd*(idin-1) 
               xn=din(idin)
               if(xn.ge.logdens(1).and.xn.le.logdens(nden).and.
     :             yt.ge.temlog(1).and.yt.le.temlog(ntem)) then
                 do iline=1,kline
                   ip=3
                   eps2d(itin,idin,iline)=
     :             10.d0**TwoDnt(iline,xn,yt,ip)
                 enddo
               else
                 write(6,35)
                 write(20,35)
                 return
               endif
             enddo
          enddo

          do iline=1,kline
      
             write(20,50) iline,wv(iline),Case
 50          format(//1x,'Emission coefficients [cgs] for line ',i2,
     :                   '   Wavelength ',f10.2,'   Case ',a1//)
             write(20,60) (din(idin),idin=1,ndin)
 60          format(1x,'log10(T[K])',3x,'log10(N[cgs])=',
     :              f7.3,29(4x,f6.3))

             do itin=1,ntin

               write(20,70) tin(itin),(eps2d(itin,idin,iline),
     :             idin=1,ndin)
 70            format(5x,f7.3,14x,1p30e10.3)

             enddo
          enddo

        endif

      call flush(20)
      return

      end

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

      SUBROUTINE LAGR(F,X,FV,XV,NP,NI,ND,IER)
      IMPLICIT REAL*8(A-H,O-Z)
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  X(I)  = ARRAY OF CORRESPONDING X 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(ND),X(ND),XIQ(20)

      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 5999 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 5998 ITAU=1,LG
            XIQ(ITAU)=X(ILOW+ITAU-1)
 5998     CONTINUE
C
C  CONSTRUCT INTERPOLATION COEFFICIENTS
          FV=0.D0
          DO 5997 ITAU=1,LG
            FL=PHI(XIQ,LG,ITAU,XV)/PHI(XIQ,LG,ITAU,XIQ(ITAU))
            FV=FV+FL*F(ILOW+ITAU-1)
 5997     CONTINUE
          RETURN
        ENDIF
 5999 CONTINUE
      END

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

      DOUBLE PRECISION FUNCTION PHI(XIQ,LG,ITAU,X)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XIQ(20)
      FK=1.D0
      DO 5999 L=1,LG
        IF(L.NE.ITAU) THEN
          FK=FK*(X-XIQ(L))
        ENDIF
 5999 CONTINUE
      PHI=FK
      RETURN
      END
      
******************************************************************************

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

      implicit none

      integer mxt,mxd,mxout
      parameter (mxout=30,mxt=25,mxd=16)

      real*8 temlog(mxt),logdens(mxd),emms(mxt,mxd,mxout)
      real*8 ft(mxt),xt(mxt),fvt,yt
      integer il,in,it,ntem,nden,ier,np

      common/theory/temlog,logdens,emms,ntem,nden

*      write(6,*) 'enter 1Dt, yt ',yt

      do it=1,ntem
         ft(it)=emms(it,in,il)
         xt(it)=temlog(it)
      enddo

*      write(6,*) 'in 1Dt ntem,ft,xt',ntem,ft,xt

      call lagr(ft,xt,fvt,yt,np,ntem,mxt,ier)
      OneDt=fvt

*      write(6,*) 'in 1Dt yt,fvt,1Dt',yt,fvt,OneDt

      return

      end

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

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

      implicit none

      integer mxt,mxd,mxout
      parameter (mxout=30,mxt=25,mxd=16)

      real*8 temlog(mxt),logdens(mxd),fn(mxd),an(mxd),
     :       emms(mxt,mxd,mxout)
      real*8 fvn,xn,yt,OneDt
      integer il,np,ntem,nden,in,ier

      common/theory/temlog,logdens,emms,ntem,nden

*      write(6,*) 'enter 1Dn, il, xn, yt',il,xn,yt

      do in=1,nden
         fn(in)=OneDt(il,in,yt,np)
         an(in)=logdens(in)
      enddo

*        write(6,*) 'in 1Dn before call lagr, xn=',xn
*        write(6,*) 'fn,an',fn,an

      call lagr(fn,an,fvn,xn,np,nden,mxd,ier)
      OneDn=fvn


*      write(6,*) 'return'

      return

      end

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

      double precision function TwoDnt(il,xn,yt,np)

      implicit none

      real*8 xn,yt,OneDn
      integer il,np

*      write(6,*) 'enter TwoDnt, il,xn,yt ',il,xn,yt
      
      TwoDnt=OneDn(il,xn,yt,np)

      return
      end

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

      subroutine reademiss

* iopt = 1 read data from binary file 'lines.bin'
* iopt = 2 read data from individual lines_xxx_yyy_dr files

      implicit none

      integer mxl,mxt,mxd
      parameter (mxl=10000,mxt=25,mxd=16)

      character*1 eflag(mxl),caseA,caseB,caseC,Case,keepCase
      character*4 ktem(mxt),kdens(mxd)
      character*5 lines
      character*12 namedata
      character*25 name
      character*80 text
      character*84 txt1(mxl)
      character*32 txt2(mxl)

      real*8 em(mxt,mxd,mxl),wl(mxl),
     :  temlg(mxt),denslg(mxd),
     :  eu(mxl),el(mxl),temp(mxd)

      integer ic,it,nt,id,nd,i,icase,status,nw,access

      common/alltheory/wl,em,eflag,eu,el,txt1,txt2,temlg,denslg,
     :                 nt,nd,nw,Case

      data temlg/2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,
     :           3.0,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,
     :           4.0,4.1,4.2,4.3,4.4/
      data denslg/2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
     :             4.2,4.4,4.6,4.8,5.0/
      data  caseA/'A'/,caseB/'B'/,caseC/'C'/,keepCase/' '/
      data nt/25/, nd/16/, nw/8889/

* read line list and theoretical emission coeffs from binary 
* file 'OIIlines_AB' and store wavelengths and emission coeffs


        namedata='OIIlines_ABC'

        status=access(namedata,'')
        if(status.ne.0) then
          write(6,*) 'Data file  ',namedata,' does not exist'
          stop
        endif


        open(unit=12,status='unknown',file=namedata)

 1      write(6,10) 
 10   format(1x,'Specify Case (A, B or C)')
      read(5,*) Case

      if(Case.ne.caseA.and.Case.ne.caseB.and.Case.ne.caseC) go to 1

      if(Case.eq.keepCase) return

        write(6,*) 'open and read file  ',namedata,', please wait'

        do ic=1,18
          read(12,*)
        enddo

        do ic=1,nw

          read(12,1001) txt1(ic),wl(ic),eflag(ic),eu(ic),el(ic),
     :                  txt2(ic)

c          if(ic.eq.8889) write(6,*) txt1(ic),wl(ic)

        enddo
 1001   format(6x,a84,f11.2,4x,a1,1x,2f11.2,a32)

        do ic=1,nw

          do icase=1,3

            do i=1,3
              read(12,*)
            enddo
 
            do it=1,nt

              if(icase.eq.1.and.Case.eq.caseA.or.
     :           icase.eq.2.and.Case.eq.caseB.or.
     :           icase.eq.3.and.Case.eq.caseC) then             
                 read(12,1002) (em(it,id,ic),id=1,nd)
              else
                 read(12,1002) (temp(id),id=1,nd)
              endif

            enddo
 1002       format(26x,16(e10.3))
          enddo

      enddo
 
      write(6,40) temlg(1),temlg(nt),denslg(1),denslg(nd)
c      write(20,10) temlg(1),temlg(nt),denslg(1),denslg(nd)
 40   format(/1x,'Theoretical emission coefficients read for '/1x,
     :  'Case  ',a1/1x,
     :  'log10(Temperature[K]) = ',16x,f6.2,' to ',f6.2,
     :  ' with step 0.1'/1x,
     :  'log10(Electron number density[/cm^3]) = ',f6.2,
     :  ' to ',f6.2,' with step 0.2'/)

      keepCase=Case
      close(12)

      return
      end

 
