***************************************************************************
*                                                                         *
*     analytic fitting formula for logK(gamma,epsilon_r)                  *
*                                                                         * 
*          g=(log(gamma)-1.1505)/1.1505                                   *
*          u=(epsilon_r-5.0)/5.0                                          * 
*                                                                         *
*          logK=(a(1)g^0+a(2)g^1+...+a(11)g^10)u^0                        *
*               +(a(12)g^0+a(13)g^1+...+a(22)g^10)u^1                     *
*               +...                                                      *
*               +(a(110)g^0+a(111)g^1+...+a(121)g^10)u^10                 * 
*                                                                         *
*              for   1 <= gamma <=200 , 0 <= epsilon_r <= 10              *
*                                                                         * 
***************************************************************************
      program fitting_formula

      implicit none
      real*8 ga,er,lk

      write(6,*)'gamma=? (1<=gamma<=200)'
      read(5,*)ga

      write(6,*)'epsilon_r=? (0<=epsilon_r<=10)'
      read(5,*)er

      write(6,100)'logK(',ga,',',er,')=',lk(ga,er)
 100  format(a5,f6.2,a1,f5.2,a2,f13.6)
      
      end

c--------fitting formula of logK-------------------------------------------
      double precision function lk(ga,er)

      implicit none
      integer i,j
      real*8 g,u,ga,er,a(121)

      data a/  -0.365927E+01, 0.155574E-02, 0.620281E-02,-0.753523E-03,
     $     0.118877E-01, -0.410209E-01,  0.189159E-01,  0.510796E-01, 
     $     -0.426759E-01, -0.197340E-01,  0.200491E-01, -0.126734E+01,
     $     -0.455740E-02, -0.179841E-01,  0.365732E-01,  0.111979E-01,
     $     -0.973313E-01,  0.588188E-01,  0.110523E+00, -0.108244E+00,
     $     -0.438567E-01,  0.510687E-01,  0.607584E+00,  0.258322E-01, 
     $     0.104425E+00, -0.262102E+00, -0.105106E+00,  0.896704E+00, 
     $     -0.547015E+00, -0.109263E+01,  0.102972E+01,  0.451951E+00, 
     $     -0.493422E+00,  0.630602E-01,  0.700779E-01,  0.279543E+00, 
     $     -0.731969E+00, -0.462001E+00,  0.255956E+01, -0.122948E+01,
     $     -0.299190E+01,  0.253992E+01,  0.119137E+01, -0.122559E+01, 
     $     -0.778163E+00, -0.264033E+00, -0.107405E+01,  0.274892E+01,  
     $     0.132482E+01, -0.959895E+01,  0.538904E+01,  0.116478E+02,
     $     -0.104973E+02, -0.479595E+01,  0.506476E+01, -0.179849E+01, 
     $     -0.301686E+00, -0.121880E+01,  0.319750E+01,  0.226167E+01,
     $     -0.113011E+02,  0.480222E+01,  0.132623E+02, -0.105685E+02, 
     $     -0.530281E+01,  0.518690E+01,  0.353738E+01,  0.884245E+00, 
     $     0.362099E+01, -0.925572E+01, -0.495541E+01,  0.324488E+02, 
     $     -0.170172E+02, -0.392927E+02,  0.341813E+02,  0.161526E+02, 
     $     -0.166155E+02,  0.290085E+01,  0.472177E+00,  0.192629E+01, 
     $     -0.512232E+01, -0.403418E+01,  0.184055E+02, -0.683208E+01, 
     $     -0.216796E+02,  0.161786E+02,  0.869953E+01, -0.808293E+01, 
     $     -0.498720E+01, -0.118370E+01, -0.487413E+01,  0.125061E+02,  
     $     0.740638E+01, -0.441392E+02,  0.214258E+02,  0.533707E+02,
     $     -0.446244E+02, -0.219164E+02,  0.218847E+02, -0.180622E+01, 
     $     -0.224453E+00, -0.925366E+00,  0.253278E+01,  0.226404E+01, 
     $     -0.935145E+01,  0.291262E+01,  0.110456E+02, -0.764260E+01,
     $     -0.444246E+01,  0.390128E+01,  0.262678E+01,  0.527141E+00, 
     $     0.218368E+01, -0.566738E+01, -0.372743E+01,  0.202577E+02, 
     $     -0.899871E+01, -0.244789E+02,  0.195906E+02,  0.100465E+02, 
     $     -0.971024E+01/

      g=(log10(ga)-1.1505d0)/1.1505d0
      u=(er-5d0)/5d0
      
      lk=0d0
      
      do i=1,10
         do j=1,10
            lk=lk+a(j*11d0+i+1d0)*g**i*u**j
         end do
      end do
      
      i=0d0
      do j=1,10
         lk=lk+a(j*11d0+i+1d0)*u**j
      end do
         
      j=0d0
      do i=1,10
         lk=lk+a(j*11d0+i+1d0)*g**i
      end do

      lk=lk+a(1)
      
      end