# Maple code to define Hermite interpolating functions/coords/coeffs/mols
# $Header$

#
# Note:
# interpolation order 2 <==> fn order 3, 3-point (2nd order) derivative mols
# interpolation order 3 <==> fn order 3, 5-point (4th order) derivative mols
# interpolation order 4 <==> fn order 5, 5-point (4th order) derivative mols
#

################################################################################
################################################################################
################################################################################

#
# derivative primitives
# (argument is a procedure which should map m into the DATA() reference)
#

dx_3point :=
proc(f::procedure(integer))
(1/2) * (-f(-1) + f(+1))
end proc;

dx_5point :=
proc(f::procedure(integer))
(1/12) * (f(-2) - 8*f(-1) + 8*f(+1) - f(+2))
end proc;

################################################################################
################################################################################
################################################################################

#
# 1-D interpolating functions
#

fn_1d_order3 :=
proc(x)
+ c0 + c1*x + c2*x^2 + c3*x^3
end proc;

fn_1d_order5 :=
proc(x)
+ c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4 + c5*x^5
end proc;

################################################################################

# coordinates for 1-D interpolating functions
coord_list_1d := [x];

################################################################################

#
# coefficients in 1-D interpolating functions
#

coeff_set_1d_order3 := {c0, c1, c2, c3};
coeff_set_1d_order5 := {c0, c1, c2, c3, c4, c5};

################################################################################

#
# 1-D derivative molecules (argument = molecule center)
#

deriv_1d_dx_3point := proc(i::integer)
		      dx_3point(proc(mi::integer) DATA(i+mi) end proc)
		      end proc;
deriv_1d_dx_5point := proc(i::integer)
		      dx_5point(proc(mi::integer) DATA(i+mi) end proc)
		      end proc;

################################################################################
################################################################################
################################################################################

#
# 2-D interpolating functions
#

fn_2d_order3 :=
proc(x,y)
+ c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3
+ c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2
+ c01*y   + c11*x*y   + c21*x^2*y   + c31*x^3*y
+ c00     + c10*x     + c20*x^2     + c30*x^3
end proc;

fn_2d_order5 :=
proc(x,y)
+ c05*y^5 + c15*x*y^5 + c25*x^2*y^5 + c35*x^3*y^5 + c45*x^4*y^5 + c55*x^5*y^5
+ c04*y^4 + c14*x*y^4 + c24*x^2*y^4 + c34*x^3*y^4 + c44*x^4*y^4 + c54*x^5*y^4
+ c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3 + c43*x^4*y^3 + c53*x^5*y^3
+ c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2 + c42*x^4*y^2 + c52*x^5*y^2
+ c01*y   + c11*x*y   + c21*x^2*y   + c31*x^3*y   + c41*x^4*y   + c51*x^5*y
+ c00     + c10*x     + c20*x^2     + c30*x^3     + c40*x^4     + c50*x^5
end proc;

################################################################################

# coordinates for 2-D interpolating functions
coord_list_2d := [x,y];

################################################################################

#
# coefficients in 2-D interpolating functions
#

coeff_set_2d_order3 := {
		       c03, c13, c23, c33,
		       c02, c12, c22, c32,
		       c01, c11, c21, c31,
		       c00, c10, c20, c30
		       };

coeff_set_2d_order5 := {
		       c05, c15, c25, c35, c45, c55,
		       c04, c14, c24, c34, c44, c54,
		       c03, c13, c23, c33, c43, c53,
		       c02, c12, c22, c32, c42, c52,
		       c01, c11, c21, c31, c41, c51,
		       c00, c10, c20, c30, c40, c50
		       };

################################################################################

#
# 2-D derivative molecules (arguments = molecule center)
#

deriv_2d_dx_3point := proc(i::integer, j::integer)
		      dx_3point(
			 proc(mi::integer) DATA(i+mi,j) end proc
			       )
		      end proc;
deriv_2d_dy_3point := proc(i::integer, j::integer)
		      dx_3point(
			 proc(mj::integer) DATA(i,j+mj) end proc
			       )
		      end proc;
deriv_2d_dxy_3point := proc(i::integer, j::integer)
		       dx_3point(
			  proc(mi::integer)
			  dx_3point(proc(mj::integer) DATA(i+mi,j+mj) end proc)
			  end proc
				)
		       end proc;

deriv_2d_dx_5point := proc(i::integer, j::integer)
		      dx_5point(
			 proc(mi::integer) DATA(i+mi,j) end proc
			       )
		      end proc;
deriv_2d_dy_5point := proc(i::integer, j::integer)
		      dx_5point(
			 proc(mj::integer) DATA(i,j+mj) end proc
			       )
		      end proc;
deriv_2d_dxy_5point := proc(i::integer, j::integer)
		       dx_5point(
			  proc(mi::integer)
			  dx_5point(proc(mj::integer) DATA(i+mi,j+mj) end proc)
			  end proc
				)
		       end proc;

################################################################################
################################################################################
################################################################################

#
# 3-D interpolating functions
#

fn_3d_order3 :=
proc(x,y,z)
# z^3 ---------------------------------------------------------------
+ c033*y^3*z^3 + c133*x*y^3*z^3 + c233*x^2*y^3*z^3 + c333*x^3*y^3*z^3
+ c023*y^2*z^3 + c123*x*y^2*z^3 + c223*x^2*y^2*z^3 + c323*x^3*y^2*z^3
+ c013*y  *z^3 + c113*x*y  *z^3 + c213*x^2*y  *z^3 + c313*x^3*y  *z^3
+ c003    *z^3 + c103*x    *z^3 + c203*x^2    *z^3 + c303*x^3    *z^3
# z^2 ---------------------------------------------------------------
+ c032*y^3*z^2 + c132*x*y^3*z^2 + c232*x^2*y^3*z^2 + c332*x^3*y^3*z^2
+ c022*y^2*z^2 + c122*x*y^2*z^2 + c222*x^2*y^2*z^2 + c322*x^3*y^2*z^2
+ c012*y  *z^2 + c112*x*y  *z^2 + c212*x^2*y  *z^2 + c312*x^3*y  *z^2
+ c002    *z^2 + c102*x    *z^2 + c202*x^2    *z^2 + c302*x^3    *z^2
# z^1 ---------------------------------------------------------------
+ c031*y^3*z   + c131*x*y^3*z   + c231*x^2*y^3*z   + c331*x^3*y^3*z
+ c021*y^2*z   + c121*x*y^2*z   + c221*x^2*y^2*z   + c321*x^3*y^2*z
+ c011*y  *z   + c111*x*y  *z   + c211*x^2*y  *z   + c311*x^3*y  *z
+ c001    *z   + c101*x    *z   + c201*x^2    *z   + c301*x^3    *z
# z^0 ---------------------------------------------------------------
+ c030*y^3     + c130*x*y^3     + c230*x^2*y^3     + c330*x^3*y^3
+ c020*y^2     + c120*x*y^2     + c220*x^2*y^2     + c320*x^3*y^2
+ c010*y       + c110*x*y       + c210*x^2*y       + c310*x^3*y
+ c000         + c100*x         + c200*x^2         + c300*x^3
end proc;

fn_3d_order5 :=
proc(x,y,z)
# z^5
+ c055*y^5*z^5 + c155*x*y^5*z^5 + c255*x^2*y^5*z^5 + c355*x^3*y^5*z^5 + c455*x^4*y^5*z^5 + c555*x^5*y^5*z^5
+ c045*y^4*z^5 + c145*x*y^4*z^5 + c245*x^2*y^4*z^5 + c345*x^3*y^4*z^5 + c445*x^4*y^4*z^5 + c545*x^5*y^4*z^5
+ c035*y^3*z^5 + c135*x*y^3*z^5 + c235*x^2*y^3*z^5 + c335*x^3*y^3*z^5 + c435*x^4*y^3*z^5 + c535*x^5*y^3*z^5
+ c025*y^2*z^5 + c125*x*y^2*z^5 + c225*x^2*y^2*z^5 + c325*x^3*y^2*z^5 + c425*x^4*y^2*z^5 + c525*x^5*y^2*z^5
+ c015*y  *z^5 + c115*x*y  *z^5 + c215*x^2*y  *z^5 + c315*x^3*y  *z^5 + c415*x^4*y  *z^5 + c515*x^5*y  *z^5
+ c005    *z^5 + c105*x    *z^5 + c205*x^2    *z^5 + c305*x^3    *z^5 + c405*x^4    *z^5 + c505*x^5    *z^5
# z^4
+ c054*y^5*z^4 + c154*x*y^5*z^4 + c254*x^2*y^5*z^4 + c354*x^3*y^5*z^4 + c454*x^4*y^5*z^4 + c554*x^5*y^5*z^4
+ c044*y^4*z^4 + c144*x*y^4*z^4 + c244*x^2*y^4*z^4 + c344*x^3*y^4*z^4 + c444*x^4*y^4*z^4 + c544*x^5*y^4*z^4
+ c034*y^3*z^4 + c134*x*y^3*z^4 + c234*x^2*y^3*z^4 + c334*x^3*y^3*z^4 + c434*x^4*y^3*z^4 + c534*x^5*y^3*z^4
+ c024*y^2*z^4 + c124*x*y^2*z^4 + c224*x^2*y^2*z^4 + c324*x^3*y^2*z^4 + c424*x^4*y^2*z^4 + c524*x^5*y^2*z^4
+ c014*y  *z^4 + c114*x*y  *z^4 + c214*x^2*y  *z^4 + c314*x^3*y  *z^4 + c414*x^4*y  *z^4 + c514*x^5*y  *z^4
+ c004    *z^4 + c104*x    *z^4 + c204*x^2    *z^4 + c304*x^3    *z^4 + c404*x^4    *z^4 + c504*x^5    *z^4
# z^3
+ c053*y^5*z^3 + c153*x*y^5*z^3 + c253*x^2*y^5*z^3 + c353*x^3*y^5*z^3 + c453*x^4*y^5*z^3 + c553*x^5*y^5*z^3
+ c043*y^4*z^3 + c143*x*y^4*z^3 + c243*x^2*y^4*z^3 + c343*x^3*y^4*z^3 + c443*x^4*y^4*z^3 + c543*x^5*y^4*z^3
+ c033*y^3*z^3 + c133*x*y^3*z^3 + c233*x^2*y^3*z^3 + c333*x^3*y^3*z^3 + c433*x^4*y^3*z^3 + c533*x^5*y^3*z^3
+ c023*y^2*z^3 + c123*x*y^2*z^3 + c223*x^2*y^2*z^3 + c323*x^3*y^2*z^3 + c423*x^4*y^2*z^3 + c523*x^5*y^2*z^3
+ c013*y  *z^3 + c113*x*y  *z^3 + c213*x^2*y  *z^3 + c313*x^3*y  *z^3 + c413*x^4*y  *z^3 + c513*x^5*y  *z^3
+ c003    *z^3 + c103*x    *z^3 + c203*x^2    *z^3 + c303*x^3    *z^3 + c403*x^4    *z^3 + c503*x^5    *z^3
# z^2
+ c052*y^5*z^2 + c152*x*y^5*z^2 + c252*x^2*y^5*z^2 + c352*x^3*y^5*z^2 + c452*x^4*y^5*z^2 + c552*x^5*y^5*z^2
+ c042*y^4*z^2 + c142*x*y^4*z^2 + c242*x^2*y^4*z^2 + c342*x^3*y^4*z^2 + c442*x^4*y^4*z^2 + c542*x^5*y^4*z^2
+ c032*y^3*z^2 + c132*x*y^3*z^2 + c232*x^2*y^3*z^2 + c332*x^3*y^3*z^2 + c432*x^4*y^3*z^2 + c532*x^5*y^3*z^2
+ c022*y^2*z^2 + c122*x*y^2*z^2 + c222*x^2*y^2*z^2 + c322*x^3*y^2*z^2 + c422*x^4*y^2*z^2 + c522*x^5*y^2*z^2
+ c012*y  *z^2 + c112*x*y  *z^2 + c212*x^2*y  *z^2 + c312*x^3*y  *z^2 + c412*x^4*y  *z^2 + c512*x^5*y  *z^2
+ c002    *z^2 + c102*x    *z^2 + c202*x^2    *z^2 + c302*x^3    *z^2 + c402*x^4    *z^2 + c502*x^5    *z^2
# z^1
+ c051*y^5*z   + c151*x*y^5*z   + c251*x^2*y^5*z   + c351*x^3*y^5*z   + c451*x^4*y^5*z   + c551*x^5*y^5*z
+ c041*y^4*z   + c141*x*y^4*z   + c241*x^2*y^4*z   + c341*x^3*y^4*z   + c441*x^4*y^4*z   + c541*x^5*y^4*z
+ c031*y^3*z   + c131*x*y^3*z   + c231*x^2*y^3*z   + c331*x^3*y^3*z   + c431*x^4*y^3*z   + c531*x^5*y^3*z
+ c021*y^2*z   + c121*x*y^2*z   + c221*x^2*y^2*z   + c321*x^3*y^2*z   + c421*x^4*y^2*z   + c521*x^5*y^2*z
+ c011*y  *z   + c111*x*y  *z   + c211*x^2*y  *z   + c311*x^3*y  *z   + c411*x^4*y  *z   + c511*x^5*y  *z
+ c001    *z   + c101*x    *z   + c201*x^2    *z   + c301*x^3    *z   + c401*x^4    *z   + c501*x^5    *z
# z^0
+ c050*y^5     + c150*x*y^5     + c250*x^2*y^5     + c350*x^3*y^5     + c450*x^4*y^5     + c550*x^5*y^5
+ c040*y^4     + c140*x*y^4     + c240*x^2*y^4     + c340*x^3*y^4     + c440*x^4*y^4     + c540*x^5*y^4
+ c030*y^3     + c130*x*y^3     + c230*x^2*y^3     + c330*x^3*y^3     + c430*x^4*y^3     + c530*x^5*y^3
+ c020*y^2     + c120*x*y^2     + c220*x^2*y^2     + c320*x^3*y^2     + c420*x^4*y^2     + c520*x^5*y^2
+ c010*y       + c110*x*y       + c210*x^2*y       + c310*x^3*y       + c410*x^4*y       + c510*x^5*y
+ c000         + c100*x         + c200*x^2         + c300*x^3         + c400*x^4         + c500*x^5
end proc;

################################################################################

# coordinates for 3-D interpolating functions
coord_list_3d := [x,y,z];

################################################################################

#
# coefficients in 3-D interpolating functions
#

coeff_set_3d_order3 := {
		       # z^3
		       c033, c133, c233, c333,
		       c023, c123, c223, c323,
		       c013, c113, c213, c313,
		       c003, c103, c203, c303,
		       # z^2
		       c032, c132, c232, c332,
		       c022, c122, c222, c322,
		       c012, c112, c212, c312,
		       c002, c102, c202, c302,
		       # z^1
		       c031, c131, c231, c331,
		       c021, c121, c221, c321,
		       c011, c111, c211, c311,
		       c001, c101, c201, c301,
		       # z^0
		       c030, c130, c230, c330,
		       c020, c120, c220, c320,
		       c010, c110, c210, c310,
		       c000, c100, c200, c300
		       };

coeff_set_3d_order5 := {
		       # z^5
		       c055, c155, c255, c355, c455, c555,
		       c045, c145, c245, c345, c445, c545,
		       c035, c135, c235, c335, c435, c535,
		       c025, c125, c225, c325, c425, c525,
		       c015, c115, c215, c315, c415, c515,
		       c005, c105, c205, c305, c405, c505,
		       # z^4
		       c054, c154, c254, c354, c454, c554,
		       c044, c144, c244, c344, c444, c544,
		       c034, c134, c234, c334, c434, c534,
		       c024, c124, c224, c324, c424, c524,
		       c014, c114, c214, c314, c414, c514,
		       c004, c104, c204, c304, c404, c504,
		       # z^3
		       c053, c153, c253, c353, c453, c553,
		       c043, c143, c243, c343, c443, c543,
		       c033, c133, c233, c333, c433, c533,
		       c023, c123, c223, c323, c423, c523,
		       c013, c113, c213, c313, c413, c513,
		       c003, c103, c203, c303, c403, c503,
		       # z^2
		       c052, c152, c252, c352, c452, c552,
		       c042, c142, c242, c342, c442, c542,
		       c032, c132, c232, c332, c432, c532,
		       c022, c122, c222, c322, c422, c522,
		       c012, c112, c212, c312, c412, c512,
		       c002, c102, c202, c302, c402, c502,
		       # z^1
		       c051, c151, c251, c351, c451, c551,
		       c041, c141, c241, c341, c441, c541,
		       c031, c131, c231, c331, c431, c531,
		       c021, c121, c221, c321, c421, c521,
		       c011, c111, c211, c311, c411, c511,
		       c001, c101, c201, c301, c401, c501,
		       # z^0
		       c050, c150, c250, c350, c450, c550,
		       c040, c140, c240, c340, c440, c540,
		       c030, c130, c230, c330, c430, c530,
		       c020, c120, c220, c320, c420, c520,
		       c010, c110, c210, c310, c410, c510,
		       c000, c100, c200, c300, c400, c500
		       };

################################################################################

#
# 3-D derivative molecules (arguments = molecule center)
#

deriv_3d_dx_3point := proc(i::integer, j::integer, k::integer)
		      dx_3point(
			 proc(mi::integer) DATA(i+mi,j,k) end proc
			       )
		      end proc;
deriv_3d_dy_3point := proc(i::integer, j::integer, k::integer)
		      dx_3point(
			 proc(mj::integer) DATA(i,j+mj,k) end proc
			       )
		      end proc;
deriv_3d_dz_3point := proc(i::integer, j::integer, k::integer)
		      dx_3point(
			 proc(mk::integer) DATA(i,j,k+mk) end proc
			       )
		      end proc;
deriv_3d_dxy_3point := proc(i::integer, j::integer, k::integer)
		       dx_3point(
			  proc(mi::integer)
			  dx_3point(
			     proc(mj::integer) DATA(i+mi,j+mj,k) end proc
				   )
			  end proc
				)
		       end proc;
deriv_3d_dxz_3point := proc(i::integer, j::integer, k::integer)
		       dx_3point(
			  proc(mi::integer)
			  dx_3point(
			     proc(mk::integer) DATA(i+mi,j,k+mk) end proc
				   )
			  end proc
				)
		       end proc;
deriv_3d_dyz_3point := proc(i::integer, j::integer, k::integer)
		       dx_3point(
			  proc(mj::integer)
			  dx_3point(
			     proc(mk::integer) DATA(i,j+mj,k+mk) end proc
				   )
			  end proc
				)
		       end proc;
deriv_3d_dxyz_3point := proc(i::integer, j::integer, k::integer)
			dx_3point(
			   proc(mi::integer)
			   dx_3point(
			      proc(mj::integer)
				 dx_3point(
				    proc(mk::integer)
				    DATA(i+mi,j+mj,k+mk)
				    end proc
					  )
			      end proc
				    )
			   end proc
				 )
			end proc;

deriv_3d_dx_5point := proc(i::integer, j::integer, k::integer)
		      dx_5point(
			 proc(mi::integer) DATA(i+mi,j,k) end proc
			       )
		      end proc;
deriv_3d_dy_5point := proc(i::integer, j::integer, k::integer)
		      dx_5point(
			 proc(mj::integer) DATA(i,j+mj,k) end proc
			       )
		      end proc;
deriv_3d_dz_5point := proc(i::integer, j::integer, k::integer)
		      dx_5point(
			 proc(mk::integer) DATA(i,j,k+mk) end proc
			       )
		      end proc;
deriv_3d_dxy_5point := proc(i::integer, j::integer, k::integer)
		       dx_5point(
			  proc(mi::integer)
			  dx_5point(
			     proc(mj::integer) DATA(i+mi,j+mj,k) end proc
				   )
			  end proc
				)
		       end proc;
deriv_3d_dxz_5point := proc(i::integer, j::integer, k::integer)
		       dx_5point(
			  proc(mi::integer)
			  dx_5point(
			     proc(mk::integer) DATA(i+mi,j,k+mk) end proc
				   )
			  end proc
				)
		       end proc;
deriv_3d_dyz_5point := proc(i::integer, j::integer, k::integer)
		       dx_5point(
			  proc(mj::integer)
			  dx_5point(
			     proc(mk::integer) DATA(i,j+mj,k+mk) end proc
				   )
			  end proc
				)
		       end proc;
deriv_3d_dxyz_5point := proc(i::integer, j::integer, k::integer)
			dx_5point(
			   proc(mi::integer)
			   dx_5point(
			      proc(mj::integer)
				 dx_5point(
				    proc(mk::integer)
				    DATA(i+mi,j+mj,k+mk)
				    end proc
					  )
			      end proc
				    )
			   end proc
				 )
			end proc;

################################################################################
################################################################################
################################################################################