program derivative C C nx = # of points in X, nq = # of points in Q^2 C have to be adjusted (in 2 places) in the program depending on the C values you have! C implicit none integer ordn integer nx,nq parameter(nx=58,nq=10) real*8 xx(nx),qq(nq) real*8 dum1,dum2,err integer i,j,k1,iknl real*8 x,q,deriv,reufunc,imufunc,redfunc,imdfunc,resfunc,imsfunc common /kins/ xx,qq real*8 reu(nx,nq),imu(nx,nq),red(nx,nq),imd(nx,nq) real*8 res(nx,nq),ims(nx,nq) common /amps/ reu,imu,red,imd,res,ims external deriv,reufunc,imufunc,redfunc,imdfunc,resfunc,imsfunc do k1 = 1,2 C C computes derivative for unpolarized (iknl=1) and polarized (iknl=2) C amplitudes. C do iknl = 1,2 C C computes derivative for unpolarized amplitudes only! C c do iknl = 1,1 if (iknl.eq.1.and.k1.eq.1) then open(unit =11,file='luamptw3.dat',status='unknown') do j = 1, nx read(11,*) x xx(j) = x do i = 1,nq read(11,*) q,dum1,dum2 qq(i) = q reu(j,i) = dum1 imu(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='ldamptw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 red(j,i) = dum1 imd(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='lsamptw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 res(j,i)= dum1 ims(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='luamptw3d.dat',status='unknown') open(unit =12,file='ldamptw3d.dat',status='unknown') open(unit =13,file='lsamptw3d.dat',status='unknown') elseif(iknl.eq.1.and.k1.eq.2) then open(unit =11,file='luampetw3.dat',status='unknown') do j = 1, nx read(11,*) x xx(j) = x do i = 1,nq read(11,*) q,dum1,dum2 qq(i) = q reu(j,i) = dum1 imu(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='ldampetw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 red(j,i) = dum1 imd(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='lsampetw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 res(j,i)= dum1 ims(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='luampetw3d.dat',status='unknown') open(unit =12,file='ldampetw3d.dat',status='unknown') open(unit =13,file='lsampetw3d.dat',status='unknown') elseif(iknl.eq.2.and.k1.eq.1) then open(unit =11,file='luamppoltw3.dat',status='unknown') do j = 1, nx read(11,*) x xx(j) = x do i = 1,nq read(11,*) q,dum1,dum2 qq(i) = q reu(j,i) = dum1 imu(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='ldamppoltw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 red(j,i) = dum1 imd(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='lsamppoltw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 res(j,i)= dum1 ims(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='luamppoltw3d.dat',status='unknown') open(unit =12,file='ldamppoltw3d.dat',status='unknown') open(unit =13,file='lsamppoltw3d.dat',status='unknown') elseif(iknl.eq.2.and.k1.eq.2) then open(unit =11,file='luamppoletw3.dat',status='unknown') do j = 1, nx read(11,*) x xx(j) = x do i = 1,nq read(11,*) q,dum1,dum2 qq(i) = q reu(j,i) = dum1 imu(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='ldamppoletw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 red(j,i) = dum1 imd(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='lsamppoletw3.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 res(j,i)= dum1 ims(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='luamppoletw3d.dat',status='unknown') open(unit =12,file='ldamppoletw3d.dat',status='unknown') open(unit =13,file='lsamppoletw3d.dat',status='unknown') endif do i = 1,nx write(11,*) xx(i) write(12,*) xx(i) write(13,*) xx(i) do j = 1,nq write(11,102) qq(j),xx(i)*deriv(reufunc,j,xx(i),xx(i)/10D0,err), > xx(i)*deriv(imufunc,j,xx(i),xx(i)/10D0,err) write(12,102) qq(j),xx(i)*deriv(redfunc,j,xx(i),xx(i)/10D0,err), > xx(i)*deriv(imdfunc,j,xx(i),xx(i)/10D0,err) write(13,102) qq(j),xx(i)*deriv(resfunc,j,xx(i),xx(i)/10D0,err), > xx(i)*deriv(imsfunc,j,xx(i),xx(i)/10D0,err) enddo enddo close(unit=13) close(unit=12) close(unit=11) enddo enddo 102 FORMAT(3(E15.8,1X)) end C C interpolation routine to determine value of a function at point x C real*8 function reufunc(x,nqtemp) implicit none integer nx,nq,i,j,count,nqtemp real*8 temp,err,x real*8 imufunc,redfunc,imdfunc,resfunc,imsfunc parameter(nx=58,nq=10) real*8 xx(nx),qq(nq) real*8 tempv(5),tempx(5) real*8 reu(nx,nq),imu(nx,nq),red(nx,nq),imd(nx,nq) real*8 res(nx,nq),ims(nx,nq) common /amps/ reu,imu,red,imd,res,ims common /kins/ xx,qq do j=1,nx-1 if(x.ge.xx(j).and.x.lt.xx(j+1)) count = j enddo do i=1,5 if (count.lt.3) then tempv(i) = reu(i,nqtemp) tempx(i) = xx(i) elseif(count.gt.nx-2) then tempv(i) = reu(nx-i+1,nqtemp) tempx(i) = xx(nx-i+1) else tempv(i) = reu(count+3-i,nqtemp) tempx(i) = xx(count+3-i) endif enddo call polint(tempx,tempv,5,x,temp,err) reufunc = temp c print *, x,tempv,temp c read(*,*) return entry imufunc(x,nqtemp) do j=1,nx-1 if(x.ge.xx(j).and.x.lt.xx(j+1)) count = j enddo do i=1,5 if (count.lt.3) then tempv(i) = imu(i,nqtemp) tempx(i) = xx(i) elseif(count.gt.nx-2) then tempv(i) = imu(nx-i+1,nqtemp) tempx(i) = xx(nx-i+1) else tempv(i) = imu(count+3-i,nqtemp) tempx(i) = xx(count+3-i) endif enddo call polint(tempx,tempv,5,x,temp,err) imufunc = temp c print *,x,tempv,temp c read(*,*) return entry redfunc(x,nqtemp) do j=1,nx-1 if(x.ge.xx(j).and.x.lt.xx(j+1)) count = j enddo do i=1,5 if (count.lt.3) then tempv(i) = red(i,nqtemp) tempx(i) = xx(i) elseif(count.gt.nx-2) then tempv(i) = red(nx-i+1,nqtemp) tempx(i) = xx(nx-i+1) else tempv(i) = red(count+3-i,nqtemp) tempx(i) = xx(count+3-i) endif enddo call polint(tempx,tempv,5,x,temp,err) redfunc = temp return entry imdfunc(x,nqtemp) do j=1,nx-1 if(x.ge.xx(j).and.x.lt.xx(j+1)) count = j enddo do i=1,5 if (count.lt.3) then tempv(i) = imd(i,nqtemp) tempx(i) = xx(i) elseif(count.gt.nx-2) then tempv(i) = imd(nx-i+1,nqtemp) tempx(i) = xx(nx-i+1) else tempv(i) = imd(count+3-i,nqtemp) tempx(i) = xx(count+3-i) endif enddo call polint(tempx,tempv,5,x,temp,err) imdfunc = temp return entry resfunc(x,nqtemp) do j=1,nx-1 if(x.ge.xx(j).and.x.lt.xx(j+1)) count = j enddo do i=1,5 if (count.lt.3) then tempv(i) = res(i,nqtemp) tempx(i) = xx(i) elseif(count.gt.nx-2) then tempv(i) = res(nx-i+1,nqtemp) tempx(i) = xx(nx-i+1) else tempv(i) = res(count+3-i,nqtemp) tempx(i) = xx(count+3-i) endif enddo call polint(tempx,tempv,5,x,temp,err) resfunc = temp return entry imsfunc(x,nqtemp) do j=1,nx-1 if(x.ge.xx(j).and.x.lt.xx(j+1)) count = j enddo do i=1,5 if (count.lt.3) then tempv(i) = ims(i,nqtemp) tempx(i) = xx(i) elseif(count.gt.nx-2) then tempv(i) = ims(nx-i+1,nqtemp) tempx(i) = xx(nx-i+1) else tempv(i) = ims(count+3-i,nqtemp) tempx(i) = xx(count+3-i) endif enddo call polint(tempx,tempv,5,x,temp,err) imsfunc = temp return end C C function to compute the derivative of a function C real*8 function deriv(func,nqtemp,x,h,err) implicit none integer ntab,nqtemp real*8 err,h,x,func,con,con2,big,safe parameter(con=1.4D0,con2=con*con,big=1.E30,ntab=10,safe=2D0) external func integer i,j real*8 errt,fac,hh,a(ntab,ntab) hh=h a(1,1) = (func(x+hh,nqtemp)-func(x-hh,nqtemp))/(2D0*hh) err = big do i=2,ntab hh=hh/con a(1,i) = (func(x+hh,nqtemp)-func(x-hh,nqtemp))/(2D0*hh) fac = con2 do j =2,i a(j,i) = (a(j-1,i)*fac-a(j-1,i-1))/(fac-1D0) fac = fac*con2 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1))) if (errt.le.err) then err = errt deriv = a(j,i) endif enddo if (abs(a(i,i)-a(i-1,i-1)).ge.safe*err) return enddo return end SUBROUTINE POLINT (XA,YA,N,X,Y,DY) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (NMAX=10) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) NS=1 DIF=ABS(X-XA(1)) DO 11 I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.0.) PAUSE c THEN c print *,X c PAUSE c endif DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END