CCC# c/*@@ CCC# c @file lineInterpers.F CCC# c @author Paul Walker CCC# c @date Halloween 1996! CCC# c @desc CCC# c Line interpolaters which are called by the various CCC# c prolong routines, @seeroutine prolong1, @seeroutine CCC# c prolong2, and @seeroutine prolong3. CCC# c @enddesc CCC# c @rcsid $Id$ CCC# c@@*/ CCC# CCC# c/*@@ CCC# c @routine lineInterp CCC# c @author Paul Walker CCC# c @date Halloween 1996! CCC# c @desc various 1D interpolators CCC# c These are various 1D inteprolators of different CCC# c orders. With the new interpolators with the sweep methods CCC# c to add a higher order interpolator, all you need to do CCC# c is add a high order 1D interpolator, and it will be called CCC# c with the appropriate temporaries in 2 and 3D. CCC# c

CCC# c Currently orders are ... CCC# c @enddesc CCC# c @calledby prolong1, prolong2, prolong3 CCC# c @history CCC# c @endhistory CCC# c @comment CCC# c The input arguments are the dagh defaults; @seeroutine CCC# c prolong1 is simply a passthrough to this. CCC# c @endcomment CCC# c CCC# c @par prolong_order CCC# c @pdesc Prolongation order. Use this sensibly if you CCC# c add a new prolong! CCC# c @ptype int CCC# c @pvalues 1 CCC# c @endpar CCC# c CCC# c@@*/ subroutine lineInterp(data_from, lb_from, ub_from, shape_from, $ data_to, lb_to, ub_to, shape_to, $ lb_r, ub_r, shape_r) implicit none integer int_order c DAGH Variables integer lb_from(1), ub_from(1), shape_from(1) real*8 data_from(shape_from(1)) integer lb_to(1), ub_to(1), shape_to(1) real*8 data_to(shape_to(1)) integer lb_r(1), ub_r(1), shape_r(1) integer i, targ, idxf, ct,cf real*8 frac c Integer differentials integer dif, dit c Integer start and end point of to in from integer to_start, to_end, loc_from, loc_to, j real*8 f(3),x(3),x0,quad,xf,dx,cubic c integer geti save int_order data int_order/-1/ if (int_order .eq. -1) then c int_order = geti("prolong_order") int_order = 1 endif dit = (ub_to(1) -lb_to(1)) / (shape_to(1)-1) dif = (ub_from(1)-lb_from(1))/ (shape_from(1)-1) to_start = (lb_r(1)-lb_to(1))/dit + 1 to_end = (ub_r(1) - lb_to(1))/dit + 1 c Bounds check if (to_start .lt. 1) then write (*,*) "Wrong lower bound!" to_start = 1 endif if (to_end .gt. shape_to(1)) then write (*,*) "Wrong upper bound!" to_end = shape_to(1) endif if (shape_from(1) .lt. 2) then write (*,*) "Nothing to interpolate!" return endif c Linear if (int_order.eq.1 .or. shape_from(1) .eq. 2) then do i=to_start, to_end ct = (i-1)*dit + lb_to(1) idxf = (ct-lb_from(1))/dif + 1 if (idxf+1 .gt. shape_from(1)) then idxf = idxf-1 endif cf = (idxf-1)*dif + lb_from(1) frac = 1.0D0*(ct-cf)/(1.0D0*dif) c Bounds check if (idxf .ge. shape_from(1)) then write (*,*) "IDX Bound Fail!" idxf = shape_from(1) endif if (idxf .lt. 1) then write (*,*) "IDX Bound lower Fail!" idxf = 1 endif data_to(i) = (1.0-frac)*data_from(idxf) + $ frac*data_from(idxf+1) enddo else c Quadratic if (int_order .eq. 2) then do i=to_start, to_end targ = (((i-1)*dit+lb_to(1))-lb_from(1))/dif + 1 if (targ .gt. shape_from(1)) then targ = targ - 1 endif if (targ+1 .gt. shape_from(1)) then targ = targ -1 endif if (targ-1 .lt. 1) then targ = targ + 1 endif do j=1,3 f(j) = data_from(targ+j-2) if (targ+j-2 .lt. 1) then write (*,*) "LOWER BOUND INTERP ERROR" stop endif if (targ+j-2 .gt. shape_from(1)) then write (*,*) "UPPER BOUND INTERP ERROR" write (*,*) j, targ write (*,*) lb_from(1), ub_from(1), shape_from(1) stop endif enddo loc_from = (targ-2)*dif + lb_from(1) loc_to = (i-1)*dit + lb_to(1) x0 = 1.0 * loc_from dx = 1.0 * dif xf = 1.0 * loc_to data_to(i) = quad(f,x0,dx,xf) enddo else if (int_order .eq. 3) then do i=to_start,to_end targ = (((i-1)*dit+lb_to(1))-lb_from(1))/dif + 1 if (targ .gt. shape_from(1)) then targ = targ - 1 endif if (targ+1 .gt. shape_from(1)) then targ = targ -1 endif if (targ-1 .lt. 1) then targ = targ + 1 endif do j=1,4 f(j) = data_from(targ+j-2) enddo loc_from = (targ-2)*dif + lb_from(1) loc_to = (i-1)*dit + lb_to(1) x0 = 1.0 * loc_from dx = 1.0 * dif xf = 1.0 * loc_to data_to(i) = cubic(f, x0, dx, xf) enddo endif endif endif return end real*8 function quad(f, x0, dx, xf) implicit none real*8 f(3),x0, dx, xf 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 + dx*f2*x0 + c - f0*x0**2 - 2*f1*x0**2 + f2*x0**2)/(2*dx**2)), c - Rule(b,(-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0 + 4*f1*x0 - 2*f2*x0)/ c - (2*dx**2)),Rule(a,(f0 - 2*f1 + 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*dx2*f0 + dx*x0*(3*f0 - 4*f1 + f2) + $ f0*x0**2 - 2*f1*x0**2 + f2*x0**2)*o2dx2 b = (-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0 + $ 4*f1*x0 - 2*f2*x0)*o2dx2 a = (f0 - 2*f1 + f2)*o2dx2 quad = a*xf**2 + b*xf + c end function cubic(f, x0, dx, xf) implicit none real*8 a,b,c,d,cubic real*8 f(4),x0,dx,xf 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 = ((a*xf + b)*xf + c)*xf + d return end