# daubdemo.map - Demonstration of facilities in the daub.map package.
#
# Usage:
#          read  `daubdemo.map`;
#
# (Copyright) Jacques Gelinas, Ph.D., Maths, CMR, aout 1994, gelinas@cmr.ca
#

#read `daub.m`;
read `daub.mpl`;

print( `daubg(p) - Plot scaling function phid and wavelet psid.` );
print( `daubG(p) - Plot antiderivatives Phid and Psid.` );
print( `daubDg(p) - Plot derivatives Dphid and Dpsid.` );

# Reduce the arguments precision to speed up plots.
daubgg := proc(p,f,g,titl)                 local old, x; global Ndyad;
  old := Ndyad; Ndyad:=12;
  plot( { 'f'(p,x), 'g'(p,x) }, x=0..2*p,
         title=`Daubechies, `.titl.`, p=`.p );
  Ndyad := old; ""
end:

daubg  :=  p ->  daubgg( p, 'phid', 'psid', `functions phi, psi`  ):
daubG  :=  p ->  daubgg( p, Phid, Psid, `antiderivatives Phi, Psi` ):
daubDg :=  p ->  daubgg( p, Dphid,Dpsid,`derivatives Dphi, Dpsi` ):


print( `puzzle(p) - Exact formulas using a few points for moments ?` );
                                pz := proc(x,m) if m=0 then 1 else x^m fi end:
puzzle := proc( p )                         local L, m, x; L:=2*p-1; 
            print( `Error of numerical integration for [0,`.L.`] moments, `.
                         (L*8).` trapezoids, `.Digits.` Digits:`);
seq( mphid(p,m) - rich(0,L,L*2, unapply(pz(x,m) * 'phid'(p,x),x)), m=0..p-1 );
print( [`phid`, op(evalf(["],1)) ]);
seq( mPhid(p,m) - rich(0,L,L*2, unapply(pz(x,m) * 'Phid'(p,x),x)), m=0..p-1 );
print( [`Phid`, op(evalf(["],1)) ]);
seq( mpsid(p,m) - rich(0,L,L*2, unapply(pz(x,m) * 'psid'(p,x),x)), m=0..L );
print( [`psid`, op(evalf(["],1)) ]);
seq( mPsid(p,m) - rich(0,L,L*2, unapply(pz(x,m) * 'Psid'(p,x),x)), m=0..L );
print( [`Psid`, op(evalf(["],1)) ]);
if p > 3 then
seq(mDphid(p,m) - rich(0,L,L*2, unapply(pz(x,m) * 'Dphid'(p,x),x)), m=0..p-1 );
print( [`Dphid`, op(evalf(["],1)) ]);
seq(mDpsid(p,m) - trap(0,L,L*2, unapply(pz(x,m) * 'Dpsid'(p,x),x)), m=0..L );
print( [`Dpsid`, op(evalf(["],1)) ])   fi;
end:

print(`Jacques Gelinas, Ph.D., Maths, CMR, aout 1994, gelinas@cmr.ca`);

puzzle(4);

daubg(2);
daubG(2);
daubDg(4);

# end of daubdemo.map
