program fdstitch

c   To checkout dfrmint operation

      integer xpoints

      real*8        eta,beta,f12,f32,f52,df12_eta,df32_eta,
     $              df52_eta,df12_beta,df32_beta,df52_beta,
     $              f12th,f32th,f52th,df12_eth,df32_eth,
     $              df52_eth,df12_bth,df32_bth,df52_bth
      common /dummy/         f12,f32,f52,df12_eta,df32_eta,
     $              df52_eta,df12_beta,df32_beta,df52_beta,
     $              f12th,f32th,f52th,df12_eth,df32_eth,
     $              df52_eth,df12_bth,df32_bth,df52_bth

      real*8 xmx, xmn, ymx, ymn, fff(18), fa, fb, etaa, etab, betaa,
     $       betab, maxerr(9), fffa(9), fffb(9), err
      real*8 ff(3), dff(6), ffth(3), dffth(6)
      real*8 fs(3), eps
      equivalence (fff(1), f12)
      equivalence (ff(1), f12)
      equivalence (dff(1), df12_eta)
      equivalence (ffth(1), f12th)
      equivalence (dffth(1), df12_eth)

      integer i, j, k, l, n, ivar, if, nc, igo, legoff, inc, nsym
      logical logf

      real lbstart, lbend, etastart, etaend, dlbeta, deta, pen,
     $     symht, legx, legy, z, delta
      character*6 flag

      integer eval, evala, evalb, interp_order
      common /evalt/ eval, interp_order

      interp_order = 6
c     write (6,'(a,i2)') ' interpolation order = ',interp_order
c     write (6,'(a)') ' eps? '
c     read (5,*) eps
      eps = 1.e-4
      igo = 1
    1 continue
c     if (igo.ne.2.and.igo.ne.3) then
c         write (6,'(a,$)') ' interpolation order ? '
c         read (5,*) interp_order
c     endif
c     if (interp_order.gt.6.or.interp_order.lt.2) then
c           write (6,'(a)') ' bad order'
c           stop
c     endif
      delta = .000001
      xpoints = 160
c     xpoints = 320
c     xpoints = 480
      if (igo.ne.2) then
          write (6,'(a,$)') ' 1 beta boundary, 2 eta boundary ? '
c         write (6,'(a,$)') ' 1 beta boundary, 2 eta boundary; delta ? '
c         read (5,*) ivar, delta
          read (5,*) ivar
          if (ivar.eq.1) then
              write (6,'(a,$)') ' LOG 10 beta value ? '
              read (5,*) lbstart
              etaend = 80
              etastart = -40
              if (lbstart.ge.-2.8.and.lbstart.le.1.0)
     $            xpoints = 3.*xpoints
              deta = (etaend - etastart)/(xpoints - 1)
          elseif (ivar.eq.2) then
              write (6,'(a,$)') ' Eta value ? '
              read (5,*) etastart
              lbend = 5
              lbstart = -7
              dlbeta = (lbend - lbstart)/(xpoints - 1)
          endif
      endif
      if (igo.ne.3) then
          write (6,'(a)') ' Value to compare: '
          write (6,'(a)')'  1    2    3      4          5        6'
          write (6,'(a)')' f12  f32  f52  df12_eta  df32_eta  df52_eta'
          write (6,'(a)')'                7             8        9'
          write (6,'(a)')'             df12_beta  df32_beta  df52_beta'
          read (5,*) if
      endif
      do i = 1,9
          maxerr(i) = 0.
      enddo
      do j = 1,xpoints
          if (ivar.eq.1) then 
              beta = lbstart
              betaa = 10.**(beta - delta)
              betab = 10.**(beta + delta)
              beta = 10.**beta
              eta = etastart + deta*(j-1)
              etaa = eta
              etab = eta
              if (lbstart.ge.-2.8.and.lbstart.le.1.0.and.
     $                                          eta.gt.-28) go to 2
          elseif (ivar.eq.2) then
              eta = etastart
              etaa = eta - delta
              etab = eta + delta
              beta = 10.**(lbstart + dlbeta*(j-1))
              betaa = beta
              betab = beta
c             if (eta.eq.-30.and.(beta.lt.10.**(-1.8).or.
c    $                            beta.gt.10.**(2))) go to 8
c             if (eta.eq.-30.and.beta.gt.10.**(3.7)) go to 8
          endif

          call calcdfi (etaa,betaa,
     $                    ff,
     $                    dff,
     $                    ffth,
     $                    dffth)
          fa = fff(if)
          do i = 1,9
              fffa(i) = fff(i)
          enddo
          evala = eval
c                      Check out subroutine for functions only
          call calcfi (etaa,betaa, fs)
          do i = 1,3
              if (fs(i).gt.(1.+eps)*ff(i).or.fs(i).lt.(1.-eps)*ff(i))
     $           write (6,'(i2, f12.7, f8.3, 1pg16.9)')
     $           i, etaa, log10(betaa), fs(i)/ff(i)
          enddo
c
          call calcfi (etab,betab, fs)
          call calcdfi (etab,betab,
     $                    ff,
     $                    dff,
     $                    ffth,
     $                    dffth)
          fb = fff(if)
          do i = 1,9
              fffb(i) = fff(i)
          enddo
          evalb = eval
c                      Check out subroutine for functions only
          call calcfi (etab,betab, fs)
          do i = 1,3
              if (fs(i).gt.(1.+eps)*ff(i).or.fs(i).lt.(1.-eps)*ff(i))
     $           write (6,'(i2, f12.7, f8.3, 1pg16.9)')
     $           i, etab, log10(betab), fs(i)/ff(i)
          enddo
c
          do i = 1,9
              if (fffb(i).ne.0.) maxerr(i) = 
     $                         max(maxerr(i), fffa(i)/fffb(i))
              if (fffa(i).ne.0.) maxerr(i) = 
     $                         max(maxerr(i), fffb(i)/fffa(i))
          enddo
          if (evalb.eq.evala) then
              if (fb.ne.0.) then 
                  if (abs(1.-fa/fb).gt.1.e-2) then
                      flag = ' #####'
                  elseif (abs(1.-fa/fb).gt.1.e-3) then
                      flag = '  <<<<'
                  elseif (abs(1.-fa/fb).gt.1.e-4) then
                      flag = '  ....'
                  else
                      flag = '      '
                  endif
                  write (6,'(2f8.3,a,i1,a,i1,2(1pg13.5),0p,f12.8,a)') 
     $                 eta, log10(beta), '**',evala, '**',evalb, 
     $                 fa, fb, fa/fb, flag
              else
                  flag = '  !!!!'
                  write (6,'(2f8.3,a,i1,a,i1,2(1pg13.5),a)') 
     $                 eta, log10(beta), 
     $                 '**',evala, '**',evalb, fa, fb, flag
              endif 
          else
              if (fb.ne.0.) then 
                  if (abs(1.-fa/fb).gt.1.e-2) then
                      flag = ' #####'
                  elseif (abs(1.-fa/fb).gt.1.e-3) then
                      flag = '  <<<<'
                  elseif (abs(1.-fa/fb).gt.1.e-4) then
                      flag = '  ....'
                  else
                      flag = '      '
                  endif
                  write (6,'(2f8.3,2i3,2(1pg13.5),0p,f12.8,a)') 
     $                         eta, log10(beta), evala, evalb, 
     $                 fa, fb, fa/fb, flag
              else
                  flag = '  !!!!'
                  write (6,'(2f8.3,2i3,2(1pg13.5),a)') 
     $                         eta, log10(beta), evala, evalb, 
     $                 fa, fb, flag
              endif 
          endif
    8 continue
      enddo
    2 continue
      write (6,'(a)') ' '
      write (6,'(a)') '      Maximum Relative Errors '
      write (6,'(a)') '   f 1/2        f 3/2        f 5/2'    
      write (6,'(3(1pg13.5))') (maxerr(l),l=1,3)
      write (6,'(2a)') 
     $         'dF1/2/deta   dF3/2/deta  dF5/2/deta  dF1/2/dbeta',
     $         '    dF3/2/dbeta   dF5/2/dbeta'
      write (6,'(6(1pg13.5))') (maxerr(l),l=4,9)
      write (6,'(a)') ' all F      all dF/ deta all dF/d beta'
      write (6,'(9(1pg13.5))') (maxerr(1)+maxerr(2)+maxerr(3))/3.,
     $                         (maxerr(4)+maxerr(5)+maxerr(6))/3.,
     $                         (maxerr(7)+maxerr(8)+maxerr(9))/3. 
      write (6,'(a)') '   all F        all dF'
      write (6,'(9(1pg13.5))') (maxerr(1)+maxerr(2)+maxerr(3))/3.,
     $                         (maxerr(4)+maxerr(5)+maxerr(6)+
     $                          maxerr(7)+maxerr(8)+maxerr(9))/6. 
      write (6,'(a)') ' '
      write (6,'(a,i1,a)') 
     $ '0 stop, 1 again, 2 same boundary, new FD value (current is ',
     $                 if,')'
      write (6,'(a)') 
     $ '3 same FD, new boundary'
      read (5,*) igo
      if (igo.ne.0) go to 1
      stop
      end