Digits := 12;

read `sinc8.mas`;	# some needed subroutines

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

#	and the weight function:
## g := unapply(simplify(1/diff(phi(x),x)),x);	# the "standard"
 g := unapply(x**(7/10)*(1-x)**(7/10), x);
## g := unapply(1, x);

# 	pick a function for testing:
F := x-> (1-x)+x**(3/2)+x**2 + sin(5*x);
prime_F := unapply(diff(F(x),x), x);

#	remap boundaries via splines:
l_spline := poly(1,0,0,0,0,1):
r_spline := poly(0,0,1,0,0,1):
prime_r_spline := unapply(diff(r_spline,x), x);
prime_l_spline := unapply(diff(l_spline,x), x);
F_new := unapply(F(x) - l_spline*F(0) - r_spline*F(1), x);
prime_F_new := unapply(diff(F_new(x),x), x);


#	Pick some parameters:
N := 7;
d := evalf(Pi);
alpha := 1;
h := evalf(sqrt(Pi*d/alpha/N));


#	define the interpolant:
f_sinc := unapply(sum(F_new(psi(k*h))/g(psi(k*h))*sinc((phi(x) - k*h)/h)*g(x), k=-N..N), x):
prime_f_sinc := unapply(diff(f_sinc(x),x), x):


#	some plots for comparison:
# approximate remapping:
p1 := plot(f_sinc(x), x=0..1, style=point):
p2 := plot(F_new(x), x=0..1, style=line):
plots[display]({p1,p2}, title = cat(`Remapped F,  N = `, N, `     Digits = `, Digits));

p1 := plot(prime_f_sinc(x), x=0..1, style=point):
p2 := plot(prime_F_new(x), x=0..1, style=line):
plots[display]({p1,p2}, title = cat(`Remapped F',  N = `, N, `    Digits = `, Digits));

# approximate the original function:
p1 := plot(f_sinc(x) + l_spline(x)*F(0) + r_spline(x)*F(1), x=0..1, style=point):
p2 := plot(F(x), x=0..1, style=line):
plots[display]({p1,p2}, title = cat(`F,   N = `, N, `    Digits = `, Digits));

p1 := plot(prime_f_sinc(x) + prime_l_spline(x)*F(0) + prime_r_spline(x)*F(1), x=0..1, style=point):
p2 := plot(prime_F(x), x=0..1, style=line):
plots[display]({p1,p2}, title = cat(`F',  N = `, N, `    Digits = `, Digits));