#	Useful routines for sinc function implementations
#	Made by Mike Hohn, hohn@math.utah.edu

`help/text/poly` := TEXT(
`poly	-- return 3rd degree polynomial satisfying the specified `,
`	end conditions`):

# ==================================================================
# Rational approximation on [a,infinity); give arguments in the form 
# (ffrom, which), where which = 1,2,3 or 4, and 
# 1 = 1,0,0,0
# 2 = 0,1,0,0  
# 3 = 0,0,1,0
# 4 = 0,0,x,1

a2inf := proc(ffrom, which)
   local c1;
   if which =1 then;
      # 	case 1:
      c1 := subs(x=x-ffrom, (x+1)/(x^2+x+1));
   elif which =2 then;
      # 	case 2:
      c1 := subs(x=x-ffrom, x/(x^2+x+1));
   elif which =3 then;
      # 	case 3:
      c1 := subs(x=x-ffrom, x^2/(x^2+x+1))
   elif which =4 then;
      # 	case 4:
      c1 := subs(x=x-ffrom, x^2/(x+1));
   else;
      ERROR(`invalid type argument`);
   fi;
   c1;
end:

# 	Tested with the following:
# > t1 := a2inf(a, 4);
# > subs(x=a, t1);
# > subs(x=a, diff(t1, x));
# > limit(t1/x, x=infinity);
# > limit(diff(t1, x), x=infinity);
# ==================================================================

poly := proc(lval,ldval, rval, rdval, ffrom, tto)
   local p,k, c, sy1,sy2,sy3,sy4,pp, t1, t2;
   p := sum('x**k*c.k', k=0..3);
   pp := diff(p, x);
   sy1 := eval(subs(x=ffrom, p)) = lval;
   sy2 := eval(subs(x=ffrom, pp)) = ldval;
   sy3 := eval(subs(x=tto, p)) = rval;
   sy4 := eval(subs(x=tto, pp)) = rdval;
   t1 := solve({sy1,sy2,sy3,sy4}, {c0,c1,c2,c3});
   eval(subs(t1, p));
end:		# ok [Nov-06-1994]


pulse := proc(y,min,max)
   local x, mi, ma;
   x:=evalf(y);
   mi := evalf(min);
   ma := evalf(max);
   if (type(x, numeric) ) then;
      if ((x<=ma) and (x>=mi)) then;
	 1;
      else;
	 0;
      fi;
   else;
      'pulse'(y,min,max);
   fi;
end:

`diff/pulse` := proc(a,b,c,x);
   0;
end:

sinc := proc(y)
   local x;
   x := (evalf(y));
   if (type(x, numeric) ) then;
      if (abs(x)<0.04) then;
	 evalhf(1-1/6*Pi^2*x^2+1/120*Pi^4*x^4);
      else;
	 evalhf(sin(Pi*x)/(Pi*x));
      fi;
   else;
      'sinc'(y);
   fi;
end:

`diff/sinc` := proc(a,x);
   sinc_p(a)*diff(a,x);
end:

sinc_p := proc(y)
   local x;
   x := (evalf(y));
   if (type(x, numeric) ) then;
      if (abs(x)<0.04) then;
	 evalhf(-1/3*Pi^2*x+1/30*Pi^4*x^3-1/840*Pi^6*x^5);
      else;
	 evalhf(cos(Pi*x)/x-sin(Pi*x)/Pi/x^2);
      fi;
   else;
      'sinc_p'(y);
   fi;
end:

`diff/sinc_p` := proc(a,x);
   sinc_pp(a)*diff(a,x);
end:

sinc_pp := proc(y)
   local x;
   x := (evalf(y));
   if (type(x, numeric) ) then;
      if (abs(x)<0.04) then;
	 evalhf(-1/3*Pi^2+1/10*Pi^4*x^2-1/168*Pi^6*x^4);
      else;
	 evalhf(-sin(Pi*x)*Pi/x-2*cos(Pi*x)/x^2+2*sin(Pi*x)/Pi/x^3);
      fi;
   else;
      'sinc_pp'(y);
   fi;
end: