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