      real*8 temparr(45),rhoarr(45),kappa5arr(45),rho
      real*8 tauarr(45),pgarr(45),pearr(45),tau,temp,pg,pe
      real*8 zlamda,kappa,kappastd(45),taustd
      real*8 dtaustd,flux,zlamdastd,bnu,getkappa,e2
      real*8 vturb,h,k
      real*8 comp(45)
      real*8 pi,dtau
      integer i,init

c line opacity initialization
      real*8 param(11,4),elem(50,5),avoigt
      integer ip,jp,ie,je,nlines,j
      real*8 opac,znu,c
      c=2.9979245e18
      ip=11
      jp=4
      ie=50
      je=5
      open(unit=11,file='specline',status='old')

c      print *,'enter number of spectral lines = 2'
c      read *, nlines
      nlines=11
      do 100 i=1,nlines
      read(11,101) (param(i,j),j=1,4)
 100  param(i,2) = c/param(i,2)
 101  format(f2.0,f9.2,2f7.2)

c      print *,'enter voigt a, turbulent velocity'
c      read *, avoigt,vturb
      avoigt=0.05
      vturb=2e5

      call atparam(elem,ie,je)

c      print *,'do you want to change abundance if not reenter value=
c     *',elem(11,1)
c      read *,elem(11,1)
c      print *,'do you want to change gf values if not reenter values=
c     *',param(1,4),'  ',param(2,4)
c      read *,param(1,4),param(2,4)

      open(unit=10,file='nelsontable',status='old')
      init=1
      zlamdastd=5000.
      k=1.38062e-16
      pi=acos(-1.)

      do i=1,45
         read(10,*) temparr(i),tauarr(i),pgarr(i),pearr(i),rhoarr(i),
     l h,kappa5arr(i)
         kappastd(i)=getkappa(zlamdastd,temparr(i),pearr(i),pgarr(i),
     l rhoarr(i))
         kappastd(i)=kappastd(i)
      enddo

c     diagnostic for tau=1
c      temp=6469
c      pe=10**1.821
c      pg=10**5.144
      
c      do zlamda=3000,15000,100
      do zlamda=5885,5900,0.05
c      do zlamda=5758,5765,0.1
         flux=0.
         tau=0.
         znu=c/zlamda
         do i=2,45
            taustd=tauarr(i)
            temp=temparr(i)
            pe=pearr(i)
            pg=pgarr(i)
            rho=rhoarr(i)
            dtaustd=tauarr(i)-tauarr(i-1)
            kappa=getkappa(zlamda,temp,pe,pg,rho)
            call opacline(temp,znu,avoigt,vturb,param,ip,jp,opac,
     *           nlines,elem,ie,je,pg,rho)
            dtau=(kappa+opac)/kappastd(i)*dtaustd
            tau=tau+dtau
            comp(i)=bnu(zlamda,temp,init)*e2(tau)*dtau*2.*pi
            flux=flux+comp(i)
         enddo
         write(13,98) zlamda,flux
c         write(14,*) zlamda,(comp(j),j=2,45) This prints the matrix
         write(22,97) opac,kappa,zlamda
      enddo

 98   format(f8.2,' ',e10.4)
 97   format(e10.4,' ',e10.4,' ',e12.6)
 96   format(e10.4,' ',e10.4,' ',e10.4,' ',e10.4,' ',e10.4,' ',e10.4)

      end

      real*8 function bnu(zlam,temp,i)

      real*8 znu,c1,c2,h,k,c,zlam,temp
      integer i
      
      h=6.62620e-27
      c=2.9979250e10
      k=1.38062e-16
      znu=c/(zlam*1e-8)

      if (i .eq. 1) then
        c1 = 2.*h/c/c
        c2 = h/k
        i=2
      endif

       bnu=c1*znu**3/(exp (c2*znu/temp) - 1.)

      return
      end

      real*8 function e2(x)
      real*8 x,num,den,a(5),b(5)
      integer*4 i

      data a/1.d0,8.5733287401,18.0590169730,8.6347608925,
     1        .2677737343/
      data b/1.d0,9.5733223454,25.6329561486,21.0996530827,
     1       3.9584969228/

      num=0.d0
      den=0.d0
      do 100 i=1,5
      num=a(i)*x**(5-i)+num
      den=b(i)*x**(5-i)+den
 100  continue

      e1=(num/den)/(x*dexp(x))
      e2=(dexp(-x)-x*e1)

      end

      subroutine voigt(a,z,h)

      real*8 a, z, h, y, delta, pi

      pi = 3.14159265
      h = 0.0
      delta = 0.01
      
      if (a.eq.0.0) then
         h = dexp(-(z*z))
         goto 10
      endif

      do y=-20.,20.,0.01
      h = h + dexp(-(y*y))*delta/((z - y)**2 + a*a)
      enddo
            
      h = h*(a/pi)
 10   return

      end

      real*8 function getkappa(lamda,t,pe,pg,rho)
      real*8 khbf,lamda,t,pe,pg,alpha0,r,theta,ii,chi3,sum,loge
      real*8 khff,chilamda,gf,temporary,rho,nh
      real*8 khbfminus,a0,a1,a2,a3,a4,a5,a6,alphabf
      real*8 khffminus,loglamda,f0,f1,f2
      real*8 ke,alphae
      real*8 u1,u0,phi
      real*8 alpha(5),edge(5)
      real*8 kappa
      real gnprime
      integer n,n0

c Neutral h bound-free
      edge(1)=912
      edge(2)=3647
      edge(3)=8206
      edge(4)=14588
      edge(5)=22790
      alpha0=1.044e-26
      r=109677.576*1e-8
      theta=5039.77/t
      ii=13.60
      no=2
      chi3=ii*(1.-1./(n0+3)**2)
      loge=0.434294481903
      sum=0.
      do n=no,no+2
         call compgn(n,lamda,r,gnprime)
         alpha(n)=0.
         if(lamda.le.edge(n)) then
            sum=sum+gnprime/n**3*10**(-theta*chi(n))
            alpha(n)=gnprime/n**5*lamda**3*alpha0
         endif
      enddo
      temporary=sum+loge/2./theta/ii*(10**(-chi3*theta)-10**(-ii*theta))
      khbf=alpha0*lamda**3*temporary

c Neutral h free-free
      chilamda=1.2398e4/lamda
      gf=1+0.3456*(lamda*r)**(-1./3.)*(loge/theta/chilamda+0.5)
      khff=alpha0*lamda**3*gf*loge/2./theta/ii*10**(-theta*ii)
      
c H- bound-free
      a0=1.99654
      a1=-1.18267e-5
      a2=2.64243e-6
      a3=-4.40524e-10
      a4=3.23992e-14
      a5=-1.39568e-18
      a6=2.78701e-23
      alphabf=(a0+a1*lamda+a2*lamda**2+a3*lamda**3+a4*lamda**4+
     l a5*lamda**5+a6*lamda**6)*1e-18
      khbfminus=4.158e-10*alphabf*pe*theta**2.5*10**(0.754*theta)

c H- free-free
      loglamda=log10(lamda)
      f0=-2.2763-1.6850*loglamda+0.76661*loglamda**2-
     l 0.0533464*loglamda**3
      f1=15.2827-9.2846*loglamda+1.99381*loglamda**2-
     l 0.142631*loglamda**3
      f2=-197.789+190.266*loglamda-67.9775*loglamda**2+
     l 10.6913*loglamda**3-0.625151*loglamda**4
      khffminus=1e-26*pe*10**(f0+f1*log10(theta)+f2*log10(theta)**2)
      alphaff=1e-26*10**(f0+f1*log10(theta)+f2*log10(theta)**2)

c Scattering
      alphae=0.6648e-24
      ke=alphae*pe*(1.+6.3e-2)/(pg-pe)

c Total
      u1=1.
      u0=10**0.30
      phi=1.2020e9*u1/u0*theta**(-2.5)*10**(-theta*ii)
      kappa=(khbf+khff+khbfminus)*(1.-10**(-chilamda*theta))+
     l khffminus*(1.+phi/pe)**(-1)+ke

      call state(theta,dlog10(pg),nh)
      getkappa=kappa*nh/rho
c      write (14,98)lamda,alphabf,alphaff,(alpha(j),j=2,4),kappa
 98   format(f6.0,' ',e10.4,' ',e10.4,' ',e10.4,' ',e10.4,' 'e10.4
     l ,' ',e10.4)
c      print 97, khbf,khff,khbfminus,khffminus,ke
 97   format(e10.4,' ',e10.4,' ',e10.4,' ',e10.4,' 'e10.4)
      
      return
      end

      subroutine compgn(n,lamda,r,gnprime)
      real*8 lamda,r
      integer n
      gnprime=1.-0.3456*(lamda*r)**(-1./3.)*(lamda*r/n**2-0.5)
      return
      end

      function chi(n)
      chi=13.60*(1.-1./n**2)
      return
      end



	SUBROUTINE STATE(THETA,LOGPgas,Nh)
c
c	This will calculate the Nh, Number of Neutral Hydrogen
c	 k = 1.3806e-16 erg/K   	NHYDROG = 10**12.00 
C	AHe = NHelium/NHydorg
c       Stolen from Sue Lederer with her consent
c
	REAL*8 LOGPGAS,Nh,K,THETA,T,AHe
	T = 5040./THETA
	AHe = 10**10.93/(10**12.00)		
	K = 1.3806e-16
	Nh = 10**LOGPgas/(K*T*(1.+Ahe))
	RETURN
	END


c Reta's opacity procedures begin here

      Subroutine atparam(elem,ie,je)
      real*8 elem(ie,je)
      elem(11,1) = 1.65959E-10
      elem(11,2) = 2.04
      elem(11,3) = 0.00017
      elem(11,4) = 5040.0
      elem(11,5) = 22.991
      elem(28,1) = 1.99526e-6
      elem(28,2) = 29.512
      elem(28,3) = 0.00204
      elem(28,4) = 5040.0
      elem(28,5) = 58.71
      elem(22,1) = 1.349e-7
      elem(22,2) = 30.200
      elem(22,3) = 0.00906
      elem(22,4) = 5040.0
      elem(22,5) = 47.90
      elem(26,1) = 3.981e-5
      elem(26,2) = 42.658
      elem(26,3) = 0.00406
      elem(26,4) = 5040.0
      elem(26,5) = 55.847
      return
      end

      Subroutine opacline(temp, znu, avoigt, vturb, param,ip,jp,
     *opac,nlines,elem,ie,je,pg,rho)
      real*8 param(ip,jp),elem(ie,je)
      real*8 temp,znu,avoigt,vturb,opac
      real*8 vturb2,c1,c2,c3,c,znuz,iparam,u,ut
      real*8 H,pg,theta,rho,nh
      integer ip,jp,nlines,ie,je

      vturb2 = vturb * vturb
      c1 = 0.01497365
      c2 = 0.86167E-04
      c3 = 1.663395E+08
      c = 2.997929E+10
      opac=0.0
      theta=5039.77/temp
      do 1 i=1,nlines
         iparam = param(i,1)
         znuz = param(i,2)*(sqrt((c3*temp)/elem(iparam,5) + vturb2))/c
         u = (znu-param(i,2))/znuz
         if(u.lt.-60.0) go to 1
         if(u.gt.60.0) go to 1
         ielem = param(i,1)
         ut = elem(ielem,2) + elem(ielem,3)*(temp-elem(ielem,4))
         call voigt(avoigt,u,H)
         opac=opac+(elem(ielem,1)*c1*param(i,4)*
     l        dexp(-param(i,3)/(c2*temp))*H)/(ut*znuz)
         call state(theta,dlog10(pg),nh)
 1    continue
      opac=opac*nh/rho
 2    return
      end
