c      cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      cc           mkfit routine                                  cc
c      cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
       subroutine mkfit( qpai, temp, rhom )
       implicit real*8 (a-h,o-y)
       real*8 lam
c
       dimension rhom( 1 ) ,temp( 1 ) ,qpai( 1 ) 
       common /param/ ndat, n
c
            dn = float( n )
        snthtw = 0.2319d0 
            cv = 0.5d0 + 2.d0*snthtw
            ca = 0.5d0
           cvp = 1.d0  - cv
           cap = 1.d0  - ca
c
          coep = (cv**2.d0 + ca**2.d0) + dn*(cvp**2.d0 + cap**2.d0)
          coem = (cv**2.d0 - ca**2.d0) + dn*(cvp**2.d0 - cap**2.d0)
c
        do 100 i = 1, ndat
           lam = temp(i)/5.9302d9
            xi = ( rhom(i)*1.d-9 )**( 1.d0/3.d0 )/lam
           qp1 = 10.7480d0*lam**2.d0 
     &          + 0.3967d0*lam**0.5d0 + 1.005d0
           qp2 = rhom(i)/(7.692d7*lam**3.d0 + 9.715d6*lam**0.5d0)
            qp = qp1**(-1.d0)*(1.d0 + qp2)**(-0.3d0)
c
         if ( temp(i) .lt. 1.d10 )
     &    then
            a0 =  6.002d19 
            a1 =  2.084d20
            a2 =  1.872d21
            b1 =  9.383d-1
            b2 = -4.141d-1
            b3 =  5.829d-2
             c =  5.5924d0
          else
            a0 =  6.002d19 
            a1 =  2.084d20
            a2 =  1.872d21
            b1 =  1.2383d0
            b2 = -0.8141d0
            b3 =  0.d0
             c =  4.9924d0
         endif
c
c---------------------------------------------------------------
c          fp1 = ( a0 + a1*xi + a2*xi**2.d0 )*exp(-c*xi)
c
         if ( xi .gt. 1.d2 )
     &     then
             fp1 = 0.d0
           else
             fp1 = ( a0 + a1*xi + a2*xi**2.d0 )*exp(-c*xi)
         endif
c---------------------------------------------------------------
           fp2 = xi**3.d0 + b1/lam + b2/( lam**2.d0 )
     &                             + b3/( lam**3.d0 )
            fp = fp1/fp2
c
             g = 1.d0 - 13.04d0*lam**2.d0 + 133.5d0*lam**4.d0
     &                +  1534d0*lam**6.d0 + 918.6d0*lam**8.d0
c
c------------------------------------------------------------------
c       qpai(i) = 0.5d0*coep*(1.d0 + coem/coep*qp)
c    &                                   *g*exp(-2.d0/lam)*fp
c
       if( lam .lt. 7.d-3 )
     &   then
           qpai(i) = 0.d0
         else
           qpai(i) = 0.5d0*coep*(1.d0 + coem/coep*qp)
     &                                   *g*exp(-2.d0/lam)*fp
       endif 
c-----------------------------------------------------------------
 100   continue
        return
        end