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