c===================================================================c
c    analytic fitting formula for the numerical result              c
c                                                                   c
c  J(r_s,R)= a(0,0)+a(1,0)*s+a(2,0)*s^2+ �c +a(10,0)*s^10           c
c           +a(0,1)*u+a(1,1)*s*u+a(2,1)*s^2*u+ ... +a(10,1)*s^10*u  c
c           + ... + a(10,10)*s^10*u^10                              c
c                                                                   c 
c   for  0.00001 <= r_s <= 0.1 ,  0 < R <= 50                       c
c                                                                   c
c===================================================================c
      program fitting_formulae_I

      implicit none
      real*8 rs,kfr,fitj,ffj

      write(6,*)'r_s=?   (0.00001< r_s <0.1)'
      read(5,*)rs

      write(6,*)'R=?     (0< R < 50)'
      read(5,*)kfr
      
      fitj=ffj(rs,kfr)
      
      write(6,'(a6,f7.5,a3,f6.3,a2,f10.6)')
     $     'J(r_s=',rs,';R=',kfr,')=',fitj

      end

c ==== FITTING FORMULA OF J(rs,R) =========
      DOUBLE PRECISION FUNCTION ffj(rs,kfr)

      implicit none
      integer i,j
      real*8 rs,kfr,fj,c(0:10,0:10),s,u,rs0
      data c/ 2.80066E-02, 2.91425E-04, 3.71730E-04,-3.40043E-04,
     $        8.38363E-03, 2.98675E-02, 1.44775E-02,-3.96957E-02,
     $       -3.32815E-02, 1.51205E-02, 1.52622E-02,
     $       -1.34650E-02,-4.77037E-04,-1.11205E-03, 1.76471E-03,
     $       -3.40534E-03,-3.93466E-02,-4.06027E-02, 4.02288E-02,
     $        5.69571E-02,-1.18443E-02,-2.14839E-02,
     $        4.70157E-03, 3.24480E-04, 1.11414E-03,-7.27685E-04,
     $       -1.66683E-03, 1.81463E-02, 3.23311E-02,-5.91508E-03,
     $       -3.05845E-02,-2.51457E-03, 8.10578E-03,
     $       -1.62773E-03,-1.84976E-04,-7.49702E-04,-2.02979E-03,
     $       -4.20858E-03,-6.64180E-03,-8.65196E-03,-8.13296E-03,
     $       -1.66852E-03, 4.25544E-03, 2.56649E-03,
     $        3.57498E-04, 1.26065E-04, 4.08871E-04, 3.21763E-03,
     $        1.00310E-02, 9.00760E-03,-2.12560E-03,-1.31089E-04,
     $        7.99359E-03, 2.57406E-03,-1.54054E-03,
     $        2.77894E-03, 3.14205E-04, 7.01182E-04, 4.64229E-03,
     $        2.03169E-02, 2.90509E-02, 1.34492E-03,-2.02386E-02,
     $       -1.02361E-03, 1.07350E-02, 3.10744E-03,
     $       -1.84849E-03,-4.25228E-04,-7.64487E-04,-7.73692E-03,
     $       -3.47978E-02,-4.92450E-02,-2.51385E-03, 3.68866E-02,
     $        5.80492E-03,-2.02777E-02,-8.44178E-03,
     $       -5.69873E-03,-5.70927E-04,-1.83728E-03,-1.08840E-02,
     $       -3.75841E-02,-4.46316E-02, 3.37337E-03, 2.53148E-02,
     $       -1.04420E-02,-1.85917E-02,-2.77226E-03,
     $        4.05787E-03, 6.77612E-04, 1.70005E-03, 1.41242E-02,
     $        5.56369E-02, 7.19123E-02,-1.61833E-03,-5.16877E-02,
     $       -3.04257E-05, 3.19260E-02, 1.11309E-02,
     $        4.57584E-03, 3.74132E-04, 1.40476E-03, 6.32205E-03,
     $        1.91593E-02, 2.38135E-02, 3.64315E-03,-8.55152E-03,
     $        3.95689E-03, 7.62300E-03, 1.36109E-03,
     $       -3.62782E-03,-4.07058E-04,-1.30179E-03,-7.84898E-03,
     $       -2.75063E-02,-3.55769E-02,-2.90671E-03, 2.08971E-02,
     $       -2.15731E-04,-1.40661E-02,-5.00155E-03/

c~~~~~~ ANALYTIC FITTING FORMULAE ~~~~~~~~~~
      rs0=0
      if(rs.lt.0.00001)then
         rs0=1-rs
         rs=0.00001
      endif
      s=(log10(rs)+3d0)/2d0
      u=(kfr-25d0)/25d0
      fj=0
      if(s.ne.0.and.u.ne.0)then
         do j=0,10
            do i=0,10
               fj=c(i,j)*s**i*u**j+fj
            enddo
         enddo
      else if(s.ne.0)then
         j=0
         do i=0,10
            fj=c(i,j)*s**i+fj
         enddo
      else if(u.ne.0)then
         i=0
         do j=0,10
            fj=c(i,j)*u**j+fj
         enddo
      else
         fj=c(0,0)
      endif
      if(rs0.ne.0)then
         rs=1-rs0
      endif
      ffj=fj
      end