From 2995148248d0375113932157ba86d2b45e513c26 Mon Sep 17 00:00:00 2001 From: jthorn Date: Fri, 20 May 2005 13:19:18 +0000 Subject: This file is unused; move it from src/ to new directory archive/ . [I chose to keep it in archive/, rather than only in the CVS attic, because the CVS attic isn't usable with normal Unix tools ("cat, awk, and grep" et al).] git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@67 0a4070d5-58f5-498f-b6c0-2693e757fa0f --- archive/interp2.F | 213 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/interp2.F | 213 ------------------------------------------------------ 2 files changed, 213 insertions(+), 213 deletions(-) create mode 100644 archive/interp2.F delete mode 100644 src/interp2.F diff --git a/archive/interp2.F b/archive/interp2.F new file mode 100644 index 0000000..006287f --- /dev/null +++ b/archive/interp2.F @@ -0,0 +1,213 @@ +c /*@@ +c@routine interp2d +c@date Fri Feb 14 08:46:53 1997 +c@author Paul Walker +c@desc +c Interpolates from 2D data var with coordinates x and y and +c sizes nx and ny onto 1D data out with position outx and outy +c and nout points. +c

+c This has linear, quadratic and cubic interpolators in it. +c Or will one day very soon. +c@enddesc +c@calls +c@calledby numerical_axig +c@@*/ + + subroutine interp2d(var,x,y,nx,ny,out,outx,outy,nout,order) + implicit none + integer nx,ny,nout + real*8 var(nx,ny), x(nx), y(ny) + real*8 out(nout),outx(nout),outy(nout) + integer order +c Interpolation goes from ibelow to ibelow+[1,2,3] depending on order + integer i,j,ibelow,jbelow,pt + real*8 xsym,ysym,findx,findy,frac + real*8 ydir(order+1) + real*8 ft(10), xt(10) + real*8 poly2inter, quad_2d, cubic_2d + real*8 dx, dy, PI + integer twobhjsad + + PI = 3.14159265 + +c Set up the grid spacings + dx = x(2) - x(1) + dy = y(2) - y(1) + +c Loop over all out points + do pt=1,nout + ysym = 1.0D0 + xsym = 1.0D0 +c Check bounds + findx = outx(pt) + if (findx .lt. x(1)) then + write (*,*) "Below inner bound at ",pt,outx(pt) + STOP + endif + if (findx .gt. x(nx)) then + write (*,*) "Above x bounds at ",pt,outx(pt),x(nx) + STOP + endif + findy = outy(pt) + if (findy .lt. y(1)) then + write (*,*) "Below y inner bound at ",pt,outy(pt) + STOP + endif + if (findy .gt. y(ny)) then + write (*,*) "Below y inner bound at ",pt,outy(pt) + STOP + endif + +c Locate ourselves in i,j space +c do i=1,nx +c if (x(i) .lt. findx) then +c ibelow = i +c endif +c enddo +c Assume a regular grid + ibelow = (findx-x(1))/dx+1 + + if (order .eq. 3 .and. ibelow .gt. 1) then + ibelow = ibelow - 1 + endif + if (ibelow + order .gt. nx) then + ibelow = nx - order + endif + +c do i=1,ny +c if (y(i) .lt. findy) then +c jbelow = i +c endif +c enddo +c Assume a regular grid + jbelow = (findy-y(1))/dy+1 + if (order .eq. 3 .and. jbelow .gt. 1) then + jbelow = jbelow - 1 + endif + if (jbelow + order .gt. ny) then + jbelow = ny - order + endif + +c write (*,*) "PT :",findx,findy +c write (*,*) "SYM:",sym +c write (*,*) "BOUND X ",ibelow,x(ibelow),x(ibelow+1) +c write (*,*) "BOUND Y ",jbelow,y(jbelow),y(jbelow+1) + + +c So do the interpolation + if (order .eq. 1) then +c Interp in the x direction + frac = (findx-x(ibelow))/(x(ibelow+1)-x(ibelow)) + ydir(1) = frac * var(ibelow+1,jbelow) + + $ (1.0 - frac)*var(ibelow,jbelow) + ydir(2) = frac * var(ibelow+1,jbelow+1) + + $ (1.0 - frac)*var(ibelow,jbelow+1) +c And now in the y direction + frac = (findy-y(jbelow))/(y(jbelow+1)-y(jbelow)) + out(pt) = xsym * ysym * + $ (frac * ydir(2) + (1.0 - frac) * ydir(1)) + else if (order .eq. 2) then +c Load up for calls to poly2inter + do j=1,3 + do i=1,3 + ft(i) = var(ibelow+i-1,jbelow+j-1) + xt(i) = x(ibelow+i-1) + enddo + ydir(j) = quad_2d(ft,xt(1),dx,findx) + enddo + do j=1,3 + xt(j) = y(jbelow+j-1) + enddo + out(pt) = xsym * ysym*quad_2d(ydir,xt(1),dy,findy) + else if (order .eq. 3) then +c Load up for calls to cubic + do j=1,4 + do i=1,4 + ft(i) = var(ibelow+i-1,jbelow+j-1) + xt(i) = x(ibelow+i-1) + enddo + ydir(j) = cubic_2d(ft,xt(1),dx,findx) + enddo + do j=1,4 + xt(j) = y(jbelow+j-1) + enddo + out(pt) = xsym * ysym*cubic_2d(ydir,xt(1),dy,findy) + else + write (*,*) "ORDER set wrong in interp2d",order + stop + endif + enddo + + return + end + + real*8 function linear_2d(f, x0, dx, findx) + implicit none + real*8 f(2),x0,dx,findx + real*8 frac + + frac = (findx-x0)/dx + linear_2d = (frac)*f(2) + (1.0-frac)*f(1) + + return + end + + real*8 function quad_2d(f, x0, dx, findx) + implicit none + real*8 f(3),x0, dx, findx + real*8 f0,f1,f2 + real*8 a,b,c, dx2, x02, o2dx2 +c Mathematica tells us +c - List(List(Rule(c,(2*dx**2*f0 + 3*dx*f0*x0 - 4*dx*f1*x0 +c - + dx*f2*x0 + f0*x0**2 - 2*f1*x0**2 + f2*x0**2) +c - /(2*dx**2)), Rule(b,(-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0 +c - + 4*f1*x0 - 2*f2*x0)/(2*dx**2)),Rule(a,(f0 - 2*f1 + +c - f2)/(2*dx**2)))) + + f0 = f(1) + f1 = f(2) + f2 = f(3) + dx2 = dx**2 + x02 = x0**2 + o2dx2 = 1.0D0/(2.0D0*dx2) + + c = (2.0D0*dx2*f0 + dx*x0*(3.0D0*f0 - 4.0D0*f1 + f2) + + $ x02*(f0 - 2.0D0*f1 + f2))*o2dx2 + b = (dx * (-3.0D0*f0 + 4.0D0*f1 - f2) + x0 * (- 2.0D0*f0 + + $ 4.0D0*f1 - 2.0D0*f2))*o2dx2 + + a = (f0 - 2.0D0*f1 + f2)*o2dx2 + + quad_2d = a*findx**2 + b*findx + c + + return + end + + real*8 function cubic_2d(f, x0, dx, findx) + implicit none + real*8 a,b,c,d + real*8 f(4),x0,dx,findx + + a = -(f(1)-3.0*f(2)+3.0*f(3)-f(4)) / (6.0*(dx**3)) + + b = (f(1)-2.0*f(2)+f(3))/(2.0*(dx**2)) + + $ (f(1)-3.0*f(2)+3.0*f(3)-f(4))*(dx+x0)/(2.0*(dx**3)) + + c = ((dx**2)*(-11.0*f(1) + 18.0*f(2) - 9.0*f(3) + 2.0*f(4)) + + $ dx*x0* (-12.0*f(1) + 30.0*f(2) - 24.0*f(3) + 6.0*f(4)) + + $ (x0**2)*( -3.0*f(1) + 9.0*f(2) - 9.0*f(3) + 3.0*f(4))) / + $ (6.0*(dx**3)) + + d = ((dx**3)* ( 6.0*f(1) ) + + $ (dx**2)*x0*( 11.0*f(1) - 18.0*f(2) + 9.0*f(3) - 2.0*f(4)) + + $ (x0**2)*dx*( 6.0*f(1) - 15.0*f(2) + 12.0*f(3) - 3.0*f(4)) + + $ (x0**3)* ( 1.0*f(1) - 3.0*f(2) + 3.0*f(3) - 1.0*f(4)))/ + $ (6.0*(dx**3)) + + cubic_2d = ((a*findx + b)*findx + c)*findx + d + + return + end + + diff --git a/src/interp2.F b/src/interp2.F deleted file mode 100644 index 006287f..0000000 --- a/src/interp2.F +++ /dev/null @@ -1,213 +0,0 @@ -c /*@@ -c@routine interp2d -c@date Fri Feb 14 08:46:53 1997 -c@author Paul Walker -c@desc -c Interpolates from 2D data var with coordinates x and y and -c sizes nx and ny onto 1D data out with position outx and outy -c and nout points. -c

-c This has linear, quadratic and cubic interpolators in it. -c Or will one day very soon. -c@enddesc -c@calls -c@calledby numerical_axig -c@@*/ - - subroutine interp2d(var,x,y,nx,ny,out,outx,outy,nout,order) - implicit none - integer nx,ny,nout - real*8 var(nx,ny), x(nx), y(ny) - real*8 out(nout),outx(nout),outy(nout) - integer order -c Interpolation goes from ibelow to ibelow+[1,2,3] depending on order - integer i,j,ibelow,jbelow,pt - real*8 xsym,ysym,findx,findy,frac - real*8 ydir(order+1) - real*8 ft(10), xt(10) - real*8 poly2inter, quad_2d, cubic_2d - real*8 dx, dy, PI - integer twobhjsad - - PI = 3.14159265 - -c Set up the grid spacings - dx = x(2) - x(1) - dy = y(2) - y(1) - -c Loop over all out points - do pt=1,nout - ysym = 1.0D0 - xsym = 1.0D0 -c Check bounds - findx = outx(pt) - if (findx .lt. x(1)) then - write (*,*) "Below inner bound at ",pt,outx(pt) - STOP - endif - if (findx .gt. x(nx)) then - write (*,*) "Above x bounds at ",pt,outx(pt),x(nx) - STOP - endif - findy = outy(pt) - if (findy .lt. y(1)) then - write (*,*) "Below y inner bound at ",pt,outy(pt) - STOP - endif - if (findy .gt. y(ny)) then - write (*,*) "Below y inner bound at ",pt,outy(pt) - STOP - endif - -c Locate ourselves in i,j space -c do i=1,nx -c if (x(i) .lt. findx) then -c ibelow = i -c endif -c enddo -c Assume a regular grid - ibelow = (findx-x(1))/dx+1 - - if (order .eq. 3 .and. ibelow .gt. 1) then - ibelow = ibelow - 1 - endif - if (ibelow + order .gt. nx) then - ibelow = nx - order - endif - -c do i=1,ny -c if (y(i) .lt. findy) then -c jbelow = i -c endif -c enddo -c Assume a regular grid - jbelow = (findy-y(1))/dy+1 - if (order .eq. 3 .and. jbelow .gt. 1) then - jbelow = jbelow - 1 - endif - if (jbelow + order .gt. ny) then - jbelow = ny - order - endif - -c write (*,*) "PT :",findx,findy -c write (*,*) "SYM:",sym -c write (*,*) "BOUND X ",ibelow,x(ibelow),x(ibelow+1) -c write (*,*) "BOUND Y ",jbelow,y(jbelow),y(jbelow+1) - - -c So do the interpolation - if (order .eq. 1) then -c Interp in the x direction - frac = (findx-x(ibelow))/(x(ibelow+1)-x(ibelow)) - ydir(1) = frac * var(ibelow+1,jbelow) + - $ (1.0 - frac)*var(ibelow,jbelow) - ydir(2) = frac * var(ibelow+1,jbelow+1) + - $ (1.0 - frac)*var(ibelow,jbelow+1) -c And now in the y direction - frac = (findy-y(jbelow))/(y(jbelow+1)-y(jbelow)) - out(pt) = xsym * ysym * - $ (frac * ydir(2) + (1.0 - frac) * ydir(1)) - else if (order .eq. 2) then -c Load up for calls to poly2inter - do j=1,3 - do i=1,3 - ft(i) = var(ibelow+i-1,jbelow+j-1) - xt(i) = x(ibelow+i-1) - enddo - ydir(j) = quad_2d(ft,xt(1),dx,findx) - enddo - do j=1,3 - xt(j) = y(jbelow+j-1) - enddo - out(pt) = xsym * ysym*quad_2d(ydir,xt(1),dy,findy) - else if (order .eq. 3) then -c Load up for calls to cubic - do j=1,4 - do i=1,4 - ft(i) = var(ibelow+i-1,jbelow+j-1) - xt(i) = x(ibelow+i-1) - enddo - ydir(j) = cubic_2d(ft,xt(1),dx,findx) - enddo - do j=1,4 - xt(j) = y(jbelow+j-1) - enddo - out(pt) = xsym * ysym*cubic_2d(ydir,xt(1),dy,findy) - else - write (*,*) "ORDER set wrong in interp2d",order - stop - endif - enddo - - return - end - - real*8 function linear_2d(f, x0, dx, findx) - implicit none - real*8 f(2),x0,dx,findx - real*8 frac - - frac = (findx-x0)/dx - linear_2d = (frac)*f(2) + (1.0-frac)*f(1) - - return - end - - real*8 function quad_2d(f, x0, dx, findx) - implicit none - real*8 f(3),x0, dx, findx - real*8 f0,f1,f2 - real*8 a,b,c, dx2, x02, o2dx2 -c Mathematica tells us -c - List(List(Rule(c,(2*dx**2*f0 + 3*dx*f0*x0 - 4*dx*f1*x0 -c - + dx*f2*x0 + f0*x0**2 - 2*f1*x0**2 + f2*x0**2) -c - /(2*dx**2)), Rule(b,(-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0 -c - + 4*f1*x0 - 2*f2*x0)/(2*dx**2)),Rule(a,(f0 - 2*f1 + -c - f2)/(2*dx**2)))) - - f0 = f(1) - f1 = f(2) - f2 = f(3) - dx2 = dx**2 - x02 = x0**2 - o2dx2 = 1.0D0/(2.0D0*dx2) - - c = (2.0D0*dx2*f0 + dx*x0*(3.0D0*f0 - 4.0D0*f1 + f2) + - $ x02*(f0 - 2.0D0*f1 + f2))*o2dx2 - b = (dx * (-3.0D0*f0 + 4.0D0*f1 - f2) + x0 * (- 2.0D0*f0 + - $ 4.0D0*f1 - 2.0D0*f2))*o2dx2 - - a = (f0 - 2.0D0*f1 + f2)*o2dx2 - - quad_2d = a*findx**2 + b*findx + c - - return - end - - real*8 function cubic_2d(f, x0, dx, findx) - implicit none - real*8 a,b,c,d - real*8 f(4),x0,dx,findx - - a = -(f(1)-3.0*f(2)+3.0*f(3)-f(4)) / (6.0*(dx**3)) - - b = (f(1)-2.0*f(2)+f(3))/(2.0*(dx**2)) + - $ (f(1)-3.0*f(2)+3.0*f(3)-f(4))*(dx+x0)/(2.0*(dx**3)) - - c = ((dx**2)*(-11.0*f(1) + 18.0*f(2) - 9.0*f(3) + 2.0*f(4)) + - $ dx*x0* (-12.0*f(1) + 30.0*f(2) - 24.0*f(3) + 6.0*f(4)) + - $ (x0**2)*( -3.0*f(1) + 9.0*f(2) - 9.0*f(3) + 3.0*f(4))) / - $ (6.0*(dx**3)) - - d = ((dx**3)* ( 6.0*f(1) ) + - $ (dx**2)*x0*( 11.0*f(1) - 18.0*f(2) + 9.0*f(3) - 2.0*f(4)) + - $ (x0**2)*dx*( 6.0*f(1) - 15.0*f(2) + 12.0*f(3) - 3.0*f(4)) + - $ (x0**3)* ( 1.0*f(1) - 3.0*f(2) + 3.0*f(3) - 1.0*f(4)))/ - $ (6.0*(dx**3)) - - cubic_2d = ((a*findx + b)*findx + c)*findx + d - - return - end - - -- cgit v1.2.3