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