# 	Generate the functions used in sinc3n3.cc, and put them in
# 	s3mo.cc  
# 
# 
interface(errorbreak = 2);

read `sinc3.mas`;

cf := proc(expr);
   convert(expr, float);
end:

# 	Make subroutines for the splines used near the endpoints:
w_1 := unapply(poly(0,1,0,0, 0,1), x);
w_2 := unapply(poly(0,0,1,0, 0,1), x);

if (macroc_loaded <> true) then;
   read `macroc`;
fi;

optimized := true;
precision := double;
autodeclare := double;

#	first spline and external variables:
prog := [
      [includeC, `<math.h>`],
      [textC, `extern "C" {`],
      [includeC, `"matrix.h"`],
      [includeC, `"matrix2.h"`],
      [textC, `}`],
      [includeC, `"s3glob.inc"`],
      [includeC, `"s3mo.inc"`],
      [functionPP, double, ` w_1`, [[double, x]]],
      [[textC, `{`],
	    [returnC, cf(w_1(x))]],
      [textC, `}`]
      ];
writeto(`s3mo.cc`):
genC(prog):
lprint(` `):
writeto(terminal);

#	second spline:
prog := [
      [commentC, ` `],
      [commentC, ` second spline:`],
      [functionPP, double, ` w_2`, [[double, x]]],
      [[textC, `{`],
	    [declareC, double, dtmp],
	    [equalC, dtmp, cf(w_2(x))],
	    [returnC, dtmp]],
      [textC, `}`]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);


#	The shift for the BCs:
Q := poly(1,0,0,1,0,1);
P := poly(0,0,0,1,0,1);
PHI := a*Q + (b-a)*P;

phi := x -> ln(x/(1-x));		# from [0,1] to R
z = phi(u):
solve(", u):
psi := unapply(", z);			# from R to [0,1]

#	(reciprocal) Weight function:
# p1 := unapply(diff(phi(x),x),x);  	# "standard" weight
# p1 := unapply(1,x);  			# no weight
p1 := unapply(1/(x**(1/3)*(1-x)),x);	# better suited for the given
	# exact answer (x**(1/2) at left, (1-x)**(3/2) at right)

Sphi := (i,h,x) -> sinc((phi(x)-i*h)/h);

#	the main differential equation:
define(L, forall([g], L(g) = -Diff(g,x,x) + p(x)*Diff(g,x) + q(x)*g));
vv := proc(g);		# value(value(arg))
   value(value(g));
end;

p   := x -> -1/(6*x);		# the functions in the D.E.
q   := x -> 1/x**2;
f   := x -> -1/12*(-14*x^2+7*x^3+16*x^4+21*x^(5/2)*a-10*x^(5/2)*(1-x)^(1/2)*b-12*x^(3/2)*a)/x^(7/2)/(1-x)^(1/2);
exa := x -> (a+x^(1/2))*(1-x)^(3/2)+b*x;	# exact answer, for checking
check_should_be_zero := factor(vv(L(exa(x)) - f(x)));

f(x) - vv(L(PHI)):
f := unapply(", x);

# 	Define the "coefficient functions" of the actual coefficients:
# 
k := 'k';
j := 'j';	# clear for the following definitions:
# 	=======================================================
coe_spline1 := factor(vv(L(w_1(x)))*Sphi(j,h,x)/p1(x));
# 	in C:
prog := [
      [commentC, `            `],
      [commentC, ` the first spline:         `],
      [functionmPP, double, ` coe_spline1`, [[double, x]],
	    [
		  [declareC, double, dtmp],
		  [equalC, dtmp, cf(coe_spline1)],
		  [returnC, dtmp]]
	    ]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);

#	=======================================================
coe_sincs := factor(vv(L(sinc((phi(x)-k*h)/h)/p1(x)))*Sphi(j,h,x)/p1(x));
prog := [
      [commentC, `            `],
      [commentC, ` the main sinc sum:         `],
      [functionmPP, double, ` coe_sincs`, [[double, x]],
	    [
		  [declareC, double, dtmp],
		  [equalC, dtmp, cf(coe_sincs)],
		  [returnC, dtmp]]
	    ]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);


#	=======================================================
coe_spline2 := factor(vv(L(w_2(x)))*Sphi(j,h,x)/p1(x));
prog := [
      [commentC, `            `],
      [commentC, ` the second spline:         `],
      [functionmPP, double, ` coe_spline2`, [[double, x]],
	    [
		  [declareC, double, dtmp],
		  [equalC, dtmp, cf(coe_spline2)],
		  [returnC, dtmp]]
	    ]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);


#	=======================================================
# 	need psi(x):
prog := [
      [commentC, `            `],
      [commentC, `            `],
      [functionmPP, double, ` psi`, [[double, x]],
	    [
		  [declareC, double, dtmp],
		  [equalC, dtmp, cf(psi(x))],
		  [returnC, dtmp]]
	    ]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);


#	=======================================================
#	 need Sphi(j,h,x)
prog := [
      [commentC, `            `],
      [commentC, `            `],
      [functionmPP, double, ` Sphi`, [[int, j], [double, h], [double, x]],
	    [
		  [declareC, double, dtmp],
		  [equalC, dtmp, cf(Sphi(j,h,x))],
		  [returnC, dtmp]]
	    ]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);


# 	=======================================================
# 	for u, need:
piece_2:= (factor(sinc((phi(x)-k*h)/h)/p1(x))); # u_vec[k+M+2-1]
# PHI 
piece_1:= (factor(w_1(x)));			# u_vec[1-1]
piece_3:= (factor(w_2(x)));			# u_vec[m+2-1]

# 	need double u(double x)
prog := [
      [commentC, `            `],
      [commentC, `            `],
      [functionPP, double, ` u`, [[double, x]]], 
      [[textC, `{`],
	    [declarem, ` `, double, [dtmp]],
	    
	    [textC, `{ `],
	    [equalC, dtmp, cf(piece_1)*`u_vec->ve[0]`],
	    [textC, `} `],
	    
	    [textC, `{ `],
	    [equalC, dtmp, dtmp + cf(PHI)],
	    [textC, `} `],
	    
	    [textC, `{ `],
	    [equalC, dtmp, dtmp + cf(piece_3)*`u_vec->ve[m+1]`],
	    [textC, `} `],
	    
	    [form, `k=-M`, `k<=N`, `k++`,
		  [equalC, dtmp, dtmp + `u_vec->ve[k+M+1]`*cf(piece_2)]],
	    [returnC, dtmp]],
      [textC, `}`]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);


#	=======================================================
# need double u_p(double x)
fd := proc(yy);
   factor(diff(yy, x));
end;
prog := [
      [functionPP, double, ` u_p`, [[double, x]]], 
      [[textC, `{`],
	    [declarem, ` `, double, [dtmp]],
	    [textC, `{ `],
	    [equalC, dtmp, cf(fd(piece_1))*`(u_vec->ve[0])`],
	    [textC, `} `],
	    
	    [textC, `{ `],
	    [equalC, dtmp, dtmp + cf(fd(PHI))],
	    [textC, `} `],
	    
	    [textC, `{ `],
	    [equalC, dtmp, dtmp + cf(fd(piece_3))*`(u_vec->ve[m+1])`],
	    [textC, `} `],
	    
	    [form, `k=-M`, `k<=N`, `k++`,
		  [equalC, dtmp, dtmp + `(u_vec->ve[k+M+1])`*cf(fd(piece_2))]],
	    [returnC, dtmp]],
      [textC, `}`]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);

#	=======================================================
#	 need f(x)
prog := [
      [functionmPP, double, ` coe_rhs`, [[double, x]],
	    [
		  [declareC, double, dtmp],
		  [equalC, dtmp, cf(f(x)*Sphi(j,h,x)/p1(x))],
		  [returnC, dtmp]]
	    ]
      ];
appendto(`s3mo.cc`):
genC(prog):
lprint(` `):
appendto(terminal);

#	function prototypes,  included in sinc3n3.cc via s3mo.inc
prog := [
      [textC, `double sinc(double x);`],
      [textC, `double sinc_p(double x);`],
      [textC, `double sinc_pp(double x);`],
      [protoPP, double, ` coe_rhs`, [[double, x]]],
      [protoPP, double, ` w_1`, [[double, x]]],
      [protoPP, double, ` w_2`, [[double, x]]],
      [protoPP, double, ` coe_spline1`, [[double, x]]],
      [protoPP, double, ` coe_sincs`, [[double, x]]],
      [protoPP, double, ` coe_spline2`, [[double, x]]],
      [protoPP, double, ` psi`, [[double, x]]],
      [protoPP, double, ` Sphi`, [[int, j], [double, h], [double, x]]],
      [protoPP, double, ` u`, [[double, x]]], 
      [protoPP, double, ` u_p`, [[double, x]]]
      ]:
writeto(`s3mo.inc`):
genC(prog):
lprint(` `):
writeto(terminal);

interface(errorbreak = 1);