c******************************************************************************
c  This file contains the functions coulx, coulbf1s, and coulff; these 
c  routines are used by the continuous opacity routines.
c******************************************************************************



      real*8 function coulx (n,freq,z)
c******************************************************************************
c  This is function COULX from ATLAS.
c******************************************************************************`

      implicit real*8 (a-h,o-z)
      real*8 a(6), b(6), c(6)
      data a/ 0.9916, 1.105, 1.101, 1.101, 1.102, 1.0986/
      data b/  2.719d13, -2.375d14, -9.863d13, -5.765d13, -3.909d13,
     .        -2.704d13/
      data c/ -2.268d30,  4.077d28,  1.035d28,  4.593d27,  2.371d27,
     .         1.229d27/

      if (freq < z*z*3.28805d15/dfloat(n*n)) then
         coulx = 0.
      else
         coulx = 2.815d29/freq**3/dfloat(n**5)*z**4
         if (n > 6) return
         if (n == 1) then
            coulx = coulx*coulbf1s(freq,z)
         else
            coulx = coulx*(a(n)+(b(n)+c(n)*(z*z/freq))*(z*z/freq))
         endif
      endif

      return
      end





      real*8 function coulbf1s (freq,z)
c******************************************************************************
c  This is function COULBF1S from ATLAS.
c******************************************************************************`

      implicit real*8 (a-h,o-z)
      real*8 gaunt1s(151)
      data gaunt1s/
     . 0.7973,0.8094,0.8212,0.8328,0.8439,0.8548,0.8653,0.8754,0.8852,
     . 0.8946,0.9035,0.9120,0.9201,0.9278,0.9351,0.9420,0.9484,0.9544,
     . 0.9601,0.9653,0.9702,0.9745,0.9785,0.9820,0.9852,0.9879,0.9903,
     . 0.9922,0.9938,0.9949,0.9957,0.9960,0.9960,0.9957,0.9949,0.9938,
     . 0.9923,0.9905,0.9884,0.9859,0.9832,0.9801,0.9767,0.9730,0.9688,
     . 0.9645,0.9598,0.9550,0.9499,0.9445,0.9389,0.9330,0.9269,0.9206,
     . 0.9140,0.9071,0.9001,0.8930,0.8856,0.8781,0.8705,0.8627,0.8546,
     . 0.8464,0.8381,0.8298,0.8213,0.8128,0.8042,0.7954,0.7866,0.7777,
     . 0.7685,0.7593,0.7502,0.7410,0.7318,0.7226,0.7134,0.7042,0.6951,
     . 0.6859,0.6767,0.6675,0.6584,0.6492,0.6401,0.6310,0.6219,0.6129,
     . 0.6039,0.5948,0.5859,0.5769,0.5680,0.5590,0.5502,0.5413,0.5324,
     . 0.5236,0.5148,0.5063,0.4979,0.4896,0.4814,0.4733,0.4652,0.4572,
     . 0.4493,0.4415,0.4337,0.4261,0.4185,0.4110,0.4035,0.3962,0.3889,
     . 0.3818,0.3749,0.3680,0.3611,0.3544,0.3478,0.3413,0.3348,0.3285,
     . 0.3222,0.3160,0.3099,0.3039,0.2980,0.2923,0.2866,0.2810,0.2755,
     . 0.2701,0.2648,0.2595,0.2544,0.2493,0.2443,0.2394,0.2345,0.2298,
     . 0.2251,0.2205,0.2160,0.2115,0.2072,0.2029,0.1987/

      if (freq/z**2 < 3.28805d15) then
         coulbf1s = 0.
      else
         elog = log10(freq/z**2/3.28805d15)
         i = elog/0.02
         i = max(min(i+1,150),1)
         coulbf1s = gaunt1s(i) +
     .             (gaunt1s(i+1)-gaunt1s(i))/0.02*(elog-(i-1)*0.02)
      endif
      
      return
      end






      real*8 function coulff (nz,tlog,freq)
c******************************************************************************
c  This is function COULFF from ATLAS.
c******************************************************************************`

      implicit real*8 (a-h,o-z)
      real*8 z4log(6), a(11,12)
      data z4log/ 0., 1.20412, 1.90849, 2.40824, 2.79588, 3.11261/
      data a/
     . 5.53,5.49,5.46,5.43,5.40,5.25,5.00,4.69,4.48,4.16,3.85,
     . 4.91,4.87,4.84,4.80,4.77,4.63,4.40,4.13,3.87,3.52,3.27,
     . 4.29,4.25,4.22,4.18,4.15,4.02,3.80,3.57,3.27,2.98,2.70,
     . 3.64,3.61,3.59,3.56,3.54,3.41,3.22,2.97,2.70,2.45,2.20,
     . 3.00,2.98,2.97,2.95,2.94,2.81,2.65,2.44,2.21,2.01,1.81,
     . 2.41,2.41,2.41,2.41,2.41,2.32,2.19,2.02,1.84,1.67,1.50,
     . 1.87,1.89,1.91,1.93,1.95,1.90,1.80,1.68,1.52,1.41,1.30,
     . 1.33,1.39,1.44,1.49,1.55,1.56,1.51,1.42,1.33,1.25,1.17,
     . 0.90,0.95,1.00,1.08,1.17,1.30,1.32,1.30,1.20,1.15,1.11,
     . 0.55,0.58,0.62,0.70,0.85,1.01,1.15,1.18,1.15,1.11,1.08,
     . 0.33,0.36,0.39,0.46,0.59,0.76,0.97,1.09,1.13,1.10,1.08,
     . 0.19,0.21,0.24,0.28,0.38,0.53,0.76,0.96,1.08,1.09,1.09/

      gamlog = 10.39638 - tlog/1.15129 + z4log(nz)
      igam = max0(min0(int(gamlog+7.),10),1)
      hvktlg = (dlog(freq) - tlog)/1.15129 - 20.63764
      ihvkt = max0(min0(int(hvktlg+9.),11),1)
      p = gamlog - dfloat(igam-7)
      q = hvktlg - dfloat(ihvkt-9)
      coulff = (1.-p)*((1.-q)*a(igam,ihvkt) + q*a(igam,ihvkt+1)) +
     .         p*((1.-q)*a(igam+1,ihvkt) + q*a(igam+1,ihvkt+1))

      return
      end