ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c formulae.f c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine formulae ( t, r, qn, q, n) implicit real*8 (a-z) tn = 1.d1**t rn = 1.d1**r rdouble = dlog10(2*rn) x = ( 1.75d1 + rdouble - 3.d0*t )/6.d0 y = ( -2.45d1 + rdouble + 3.d0*t )/6.d0 if ( dabs(x).gt.7.d-1.or.y.lt.0.d0 ) & then fxy = 1.d0 elseif(y-1.6d0+1.25d0*x.lt.0.d0) & then fxy = 1.05d0 +(0.39d0 -1.25d0*x - 0.35d0*dsin(4.5d0*x) & -0.3d0*dexp(-1.d0*(4.5d0*x + 0.9d0)**2.d0)) & *dexp(-1.d0*((y-1.6d0+1.25d0*x) & /(0.57d0-0.25d0*x))**2.d0) else fxy = 1.05d0 +(0.39d0 - 1.25d0*x - 0.35d0*dsin(4.5d0*x) & -0.3d0*dexp(-1.d0*(4.5d0*x + 0.9d0)**2.d0)) endif lamb = (1.686d-10)*tn gamd = 1.1095d11*rn & /((tn**2.d0)*((1.d0+(1.019d-6*rn)**(2.d0/3.d0))**0.5d0)) ft = 2.4d0 + 0.6d0*(gamd**0.25d0) + 0.51d0*(gamd**0.5d0) & + 1.25d0*(gamd**0.75d0) fl = (8.6d0*gamd + 1.35d0*(gamd**(7.d0/4.d0))) & /(2.25d2 - 1.7d1*(gamd**0.5d0) + gamd ) if(dabs(gamd**0.5d0).gt.700) & then qap=0.d0 else qap = fxy*(fl+ft)*dexp(-1.d0*gamd**0.5d0)*(gamd**3.d0) & *(lamb**9.d0)*3.d21 endif cvd = ( 0.5d0 + 2.d0*0.2319d0)**2.d0 & + n*(0.5d0 - 2.d0*0.2319d0)**2.d0 if ( qap*cvd.le.0.d0 ) & then q = 9.9999999d-99 else q = dlog10(qap*cvd) qn = qap*cvd endif return end