Artifact 31c86431a80ce30cf7432fdd865f72a793e76af4e4b6310d9910c30362b1efaa:
- Executable file
r36/src/camal.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 25824) [annotate] [blame] [check-ins using] [more...]
module camal; % Package for calculations in celestial mechanics. % Author: John P. Fitch <jpff@maths.bath.ac.uk> create!-package('(camal fourdom fourplus makefour hsub hdiff), '(contrib camal)); %% This section is to define macros and simple functions to handle the %% data structures for harmonic forms. %% The structure is a vector: %% Coeff | FN | Angle | Next % %% This version only allows 8 angles. Consider extending this later. switch fourier; %% A vector and counter to record link between angle names and index global '(next!-angle!* fourier!-name!*); next!-angle!* := 0; if vectorp fourier!-name!* then << for i :=0:7 do remprop(getv(fourier!-name!*, i),'fourier!-angle) >>; fourier!-name!* := mkvect 7; %% For non Cambridge LISP add smacro procedure putv!.unsafe(x,y,z); putv(x,y,z); smacro procedure getv!.unsafe(x,y); getv(x,y); %% Data abtraction says that we should define macros for access to %% the parts of the Fourier structure smacro procedure fs!:set!-next(f,p); putv!.unsafe(f, 3, p); smacro procedure fs!:next(f); getv!.unsafe(f,3); smacro procedure fs!:set!-coeff(f,p); putv!.unsafe(f, 0, p); smacro procedure fs!:coeff(f); getv!.unsafe(f, 0); smacro procedure fs!:set!-fn(f,p); putv!.unsafe(f, 1, p); smacro procedure fs!:fn(f); getv!.unsafe(f, 1); smacro procedure fs!:set!-angle(f,p); putv!.unsafe(f, 2, p); smacro procedure fs!:angle(f); getv!.unsafe(f, 2); %% Some support functions for angle expressions symbolic procedure fs!:make!-nullangle(); begin scalar ans; ans := mkvect 7; for i:=0:7 do putv!.unsafe(ans,i,0); return ans; end; symbolic procedure fs!:null!-angle!: u; fs!:null!-angle cdr u; symbolic procedure fs!:null!-angle u; begin scalar ans, i, x; x := fs!:angle u; ans := t; i := 0; top: if not(getv!.unsafe(x,i)=0) then return nil; i := i+1; if (i<8) then go to top; return ans; end; endmodule; module fourdom; % Domain definitions for angles and fourier series % Author: John Fitch 1991. global '(domainlist!*); domainlist!*:=union('(!:fs!:),domainlist!*); put('fourier,'tag,'!:fs!:); put('!:fs!:,'dname,'fourier); flag('(!:fs!:),'field); %% Should be ring really put('!:fs!:,'i2d,'i2fourier); put('!:fs!:,'minusp,'fs!:minusp!:); put('!:fs!:,'plus,'fs!:plus!:); put('!:fs!:,'times,'fs!:times!:); put('!:fs!:, 'expt,'fs!:expt!:); put('!:fs!:,'difference,'fs!:difference!:); put('!:fs!:,'quotient,'fs!:quotient!:); put('!:fs!:, 'divide, 'fs!:divide!:); put('!:fs!:, 'gcd, 'fs!:gcd!:); put('!:fs!:,'zerop,'fs!:zerop!:); put('!:fs!:,'onep,'fs!:onep!:); put('!:fs!:,'prepfn,'fs!:prepfn!:); put('!:fs!:,'specprn,'fs!:prin!:); put('!:fs!:,'prifn,'fs!:prin!:); put('!:fs!:,'intequivfn,'fs!:intequiv!:); flag('(!:fs!:),'ratmode); % conversion functions put('!:fs!:,'!:mod!:,mkdmoderr('!:fs!:,'!:mod!:)); % put('!:fs!:,'!:gi!:,mkdmoderr('!:fs!:,'!:gi!:)); % put('!:fs!:,'!:rn!:,mkdmoderr('!:fs!:,'!:rn!:)); put('!:rn!:,'!:fs!:,'!*d2fourier); put('!:ft!:,'!:fs!:,'cdr); put('!:gi!:,'!:fs!:,'!*d2fourier); put('!:gf!:,'!:fs!:,'!*d2fourier); put('expt, '!:fs!:, 'fs!:expt!:); % Conversion functions symbolic procedure i2fourier u; if dmode!*='!:fs!: then !*d2fourier u else u; symbolic procedure !*d2fourier u; if null u then nil else begin scalar fourier; fourier:=mkvect 3; fs!:set!-coeff(fourier,(u . 1)); fs!:set!-fn(fourier,'cos); fs!:set!-angle(fourier,fs!:make!-nullangle()); fs!:set!-next(fourier,nil); return get('fourier,'tag) . fourier end; symbolic procedure !*sq2fourier u; if null car u then nil else begin scalar fourier; fourier:=mkvect 3; fs!:set!-coeff(fourier,u); fs!:set!-fn(fourier,'cos); fs!:set!-angle(fourier,fs!:make!-nullangle()); fs!:set!-next(fourier,nil); return get('fourier,'tag) . fourier end; symbolic procedure fs!:minusp!:(x); fs!:minusp cdr x; symbolic procedure fs!:minusp x; if null x then nil else if null fs!:next x then minusf car fs!:coeff x else fs!:minusp fs!:next x; %% Basic algebraic operations symbolic procedure fs!:times!:(x,y); % This function seems to be called with numeric values as well if null x then nil else if null y then nil else if numberp y then get('fourier,'tag) . fs!:timescoeff(y ./ 1, cdr x) else if numberp x then get('fourier,'tag) . fs!:timescoeff(x ./ 1, cdr y) else if not eqcar(x, get('fourier,'tag)) then get('fourier,'tag) . fs!:timescoeff(x,cdr y) else if not eqcar(y, get('fourier,'tag)) then get('fourier,'tag) . fs!:timescoeff(y,cdr x) else get('fourier,'tag) . fs!:times(cdr x, cdr y); symbolic procedure fs!:timescoeff(x, y); if null y then nil else begin scalar ans, coeff; coeff := multsq(x,fs!:coeff y); if coeff = '(nil . 1) then << print "zero in times"; return fs!:timescoeff(x, fs!:next y) >>; ans := mkvect 3; fs!:set!-coeff(ans,coeff); fs!:set!-fn(ans,fs!:fn y); fs!:set!-angle(ans,fs!:angle y); fs!:set!-next(ans, fs!:timescoeff(x, fs!:next y)); return ans end; symbolic procedure fs!:times(x,y); if null x then nil else if null y then nil else begin scalar ans; ans := fs!:timesterm(x, y); return fs!:plus(ans, fs!:times(fs!:next x, y)); end; symbolic procedure fs!:timesterm(x,y); % Treat x as a term and y as a tree if null y then nil else if null x then nil else begin scalar ans; ans := fs!:timestermterm(x,y); return fs!:plus(ans, fs!:timesterm(x, fs!:next y)); end; symbolic procedure fs!:timestermterm(x,y); % x and y are terms. Generate the two answer terms. begin scalar sum, diff, ans, xv, yv, coeff; sum := mkvect 7; xv := fs!:angle x; yv := fs!:angle y; for i:=0:7 do putv!.unsafe(sum,i, getv!.unsafe(xv,i)+getv!.unsafe(yv,i)); diff := mkvect 7; for i:=0:7 do putv!.unsafe(diff,i, getv!.unsafe(xv,i)-getv!.unsafe(yv,i)); coeff := multsq(fs!:coeff x, fs!:coeff y); coeff := multsq(coeff, '(1 . 2)); if null car coeff then return nil; if fs!:fn x = 'sin then if fs!:fn y = 'sin then % sin x*sin y => [-cos(x+y)+cos(x-y)]/2 return fs!:plus(make!-term('cos, sum, negsq coeff), make!-term('cos,diff, coeff)) else % fs!:fn y = 'cos % sin x * cos y => [sin(x+y)+sin(x-y)]/2 return fs!:plus(make!-term('sin, sum, coeff), make!-term('sin, diff,coeff)) else % fs!:fn x='cos if fs!:fn y = 'sin then % cos x*sin y => [sin(x+y)-sin(x-y)]/2 return fs!:plus(make!-term('sin, sum, coeff), make!-term('sin,diff, negsq coeff)) else % fs!:fn y = 'cos % cos x * cos y => [cos(x+y)+cos(x-y)]/2 return fs!:plus(make!-term('cos, sum, coeff), make!-term('cos, diff,coeff)) end; symbolic procedure fs!:expt!:(x,n); begin scalar ans, xx; ans := cdr !*d2fourier 1; x := cdr x; for i:=1:n do ans := fs!:times(ans,x); return get('fourier,'tag) . ans; end; symbolic procedure make!-term(fn, ang, coeff); begin scalar fourier, sign, i; sign := 0; i:=0; top: if getv!.unsafe(ang,i)<0 then sign := -1 else if getv!.unsafe(ang,i)>0 then sign := 1 else if i=7 then << if fn ='sin then return nil >> else << i := i #+ 1; goto top >>; fourier:=mkvect 3; if sign = 1 or fn = 'cos then fs!:set!-coeff(fourier,coeff) else fs!:set!-coeff(fourier, multsq('(-1 . 1), coeff)); fs!:set!-fn(fourier,fn); if sign = -1 then << sign := mkvect 7; for i:=0:7 do putv!.unsafe(sign,i,-getv!.unsafe(ang,i)); ang := sign >>; fs!:set!-angle(fourier,ang); fs!:set!-next(fourier,nil); return fourier end; symbolic procedure fs!:quotient!:(x,y); if numberp y then fs!:times!:(x, !*sq2fourier (1 ./ y)) else rerror(fourier, 98, "Unimplemented"); symbolic procedure fs!:divide!:(x,y); rerror(fourier, 98, "Unimplemented"); symbolic procedure fs!:gcd!:(x,y); rerror(fourier, 98, "Unimplemented"); symbolic procedure fs!:difference!:(x,y); fs!:plus!:(x, fs!:negate!: y); symbolic procedure fs!:negate!: x; get('fourier,'tag) . fs!:negate cdr x; symbolic procedure fs!:negate x; if null x then nil else begin scalar ans; ans := mkvect 3; fs!:set!-coeff(ans,negsq fs!:coeff x); fs!:set!-fn(ans,fs!:fn x); fs!:set!-angle(ans,fs!:angle x); fs!:set!-next(ans, fs!:negate fs!:next x); return ans end; symbolic procedure fs!:zerop!:(u); null u or (not numberp u and null cdr u or (null fs!:next cdr u and ((numberp v and zerop v) where v=fs!:coeff cdr u))); symbolic procedure fs!:onep!:(u); fs!:onep cdr u; symbolic procedure fs!:onep u; null fs!:next u and onep fs!:coeff u and fs!:null!-angle u and fs!:fn(u) = 'cos; symbolic procedure fs!:prepfn!:(x); x; symbolic procedure simpfs u; u; put('!:fs!:,'simpfn,'simpfs); %% PRINTING FUNCTIONS %% We have all the usual problems of unit coefficients, and zero angles smacro procedure zeroterm x; fs!:coeff x = '(nil . 1); symbolic procedure fs!:prin!:(x); << prin2!* "["; fs!:prin cdr x; prin2!* "]" >>; symbolic procedure fs!:prin x; if null x then prin2!* " 0 " else << while x do << fs!:prin1 x; x := fs!:next x; if x then prin2!* " + " >> >>; symbolic procedure fs!:prin1 x; begin scalar first, u, v; first := t; if not(fs!:coeff x = '(1 . 1)) then << prin2!* "("; sqprint fs!:coeff x; prin2!* ")" >>; if not(fs!:null!-angle x) then << prin2!* fs!:fn x; prin2!* "["; u := fs!:angle x; for i:=0:7 do if not((v := getv!.unsafe(u,i)) = 0) then << if v<0 then << first := t; prin2!* "-"; v := -v >>; if not first then prin2!* "+"; if not(v=1) then prin2!* v; first := nil; prin2!* getv!.unsafe(fourier!-name!*, i) >>; prin2!* "]" >> else if fs!:coeff x = '(1 . 1) then prin2!* "1" end; symbolic procedure fs!:intequiv!:(u); null fs!:next x and fs!:null!-angle x and fs!:fn(x) = 'cos and fixp car fs!:coeff x and cdr fs!:coeff x = 1 where x = cdr u; endmodule; module fourplus; %% ARITHMETIC %% Addition of Fourier expressionsis really a merge operation symbolic procedure fs!:plus!:(x,y); %% Top level addition of two fourier series if fs!:zerop!: y then x else if fs!:zerop!: x then y else get('fourier,'tag) . fs!:plus(copy!-tree cdr x, copy!-tree cdr y); % I cannot rely on the CAMAL selective copy, so I take the coward's way % out. symbolic procedure copy!-tree x; if null x then nil else begin scalar ans; ans := mkvect 3; fs!:set!-coeff(ans,fs!:coeff x); fs!:set!-fn(ans,fs!:fn x); fs!:set!-angle(ans,fs!:angle x); fs!:set!-next(ans, copy!-tree fs!:next x); return ans end; symbolic procedure fs!:plus(x, y); %% The real addition. x is a new tree to which y must be merged. if null y then x else if null x then y else if fs!:fn x = fs!:fn y and angles!-equal(fs!:angle x, fs!:angle y) then begin scalar coef; coef := addsq(fs!:coeff x, fs!:coeff y); % Really I should deal with the zero case here if null car coef then return fs!:plus(fs!:next x, fs!:next y); fs!:set!-coeff(x, coef); fs!:set!-next(x, fs!:plus(fs!:next x, fs!:next y)); return x end else if fs!:angle!-order(x, y) then << fs!:set!-next(x, fs!:plus(fs!:next x, y)); x >> else << fs!:set!-next(y, fs!:plus(fs!:next y,x)); y >>; symbolic procedure angles!-equal(x, y); % Are all angles the same? begin scalar i; i := 0; top: if not(getv!.unsafe(x,i)=getv!.unsafe(y,i)) then return nil; i := i+1; if (i<8) then go to top; return t; end; symbolic procedure fs!:angle!-order(x, y); % Ordering function for angle expressions, also taking account of angle. begin scalar ans, i, xx, yy; i := 0; xx := fs!:angle x; yy := fs!:angle y; top: ans := (getv!.unsafe(xx,i)-getv!.unsafe(yy,i)); if not(ans = 0) then return ans>0; i := i+1; if (i<8) then go to top; return if fs!:fn x = fs!:fn y then nil else if fs!:fn x = 'sin then nil else t; end; endmodule; module makefour; %% User interface; all rather iffy at present symbolic procedure harmonicp u; get(u, 'fourier!-angle); symbolic procedure harmonic u; << for each x in u do if not(get(x, 'fourier!-angle)) then << if (next!-angle!* > 7) then rerror(fourier,3,"Too many angles"); put(x, 'fourier!-angle, next!-angle!*); putv!.unsafe(fourier!-name!*, next!-angle!*, x); next!-angle!* := next!-angle!* #+ 1; >> >>; put('harmonic, 'stat, 'rlis); symbolic procedure simpfourier u; %% Handle the form fourier(...) with treating sin and cos as special begin if not(length u = 1) then rerror(fourier,1,"Argument should be single expression"); return simpfourier1 prepsq simp!* car u;; end; symbolic procedure simpfourier1 u; begin scalar ff; if atom u then << if harmonicp u then rerror(fourier,2,"Secular angle not allowed"); return (!*sq2fourier simp u) . 1; >> else if eqcar(u, '!:fs!:) then return u else if (ff := get(car u, 'simpfour)) then return apply1(ff, cdr u) else << rerror(fourier,4,"Unknown function" . car u); return (!*sq2fourier u) . 1; >> end; put('fourier, 'simpfn, 'simpfourier); symbolic procedure simpfouriersin u; % Creation of a simple angle expression and function begin scalar ans, vv; u := car u; if atom u then if harmonicp u then << ans:=mkvect 3; fs!:set!-coeff(ans,(1 . 1)); fs!:set!-fn(ans,'sin); vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); putv!.unsafe(vv, get(u, 'fourier!-angle), 1); fs!:set!-angle(ans,vv); fs!:set!-next(ans,nil); return (get('fourier,'tag) . ans) . 1 >> else return !*sq2fourier(simp list('sin, u)) . 1; if angle!-expression!-p u then << ans:=mkvect 3; fs!:set!-coeff(ans,(1 . 1)); fs!:set!-fn(ans,'sin); vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); compile!-angle!-expression(u,vv); fs!:set!-angle(ans,vv); fs!:set!-next(ans,nil); return (get('fourier,'tag) . ans) . 1 >>; rerror(fourier,99,"Not finished yet"); end; put('sin, 'simpfour, 'simpfouriersin); symbolic procedure simpfouriercos u; % Creation of a simple angle expression and function begin scalar ans, vv; u := car u; if atom u then if harmonicp u then << ans:=mkvect 3; fs!:set!-coeff(ans,(1 . 1)); fs!:set!-fn(ans,'cos); vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); putv!.unsafe(vv, get(u, 'fourier!-angle), 1); fs!:set!-angle(ans,vv); fs!:set!-next(ans,nil); return (get('fourier,'tag) . ans) . 1 >> else return !*sq2fourier(simp list('cos, u)) . 1; if angle!-expression!-p u then << ans:=mkvect 3; fs!:set!-coeff(ans,(1 . 1)); fs!:set!-fn(ans,'cos); vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); compile!-angle!-expression(u,vv); fs!:set!-angle(ans,vv); fs!:set!-next(ans,nil); return (get('fourier,'tag) . ans) . 1 >>; rerror(fourier,99,"Not finished yet"); end; put('cos, 'simpfour, 'simpfouriercos); %% Is the prefix expression u a sum of angles?? symbolic procedure angle!-expression!-p u; if atom u and harmonicp u then t else if eqcar(u,'plus) or eqcar(u,'difference) then angle!-expression!-p cadr u and angle!-expression!-p caddr u else if eqcar(u,'minus) then angle!-expression!-p cadr u else if eqcar(u,'times) then if numberp cadr u then angle!-expression!-p caddr u else angle!-expression!-p cadr u and numberp caddr u else nil; %% We know that u is a sum of angles, so create vector of coefficients. symbolic procedure compile!-angle!-expression(u,v); if atom u and harmonicp u then putv!.unsafe(v, get(u, 'fourier!-angle), 1+getv!.unsafe(v, get(u, 'fourier!-angle))) else if eqcar(u,'plus) then << u := cdr u; while u do << compile!-angle!-expression(car u,v); u := cdr u >>; v >> else if eqcar(u,'difference) then begin scalar vv; compile!-angle!-expression(cadr u,v); vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); compile!-angle!-expression(caddr u,vv); for i:=0:7 do putv!.unsafe(v,i,getv!.unsafe(v,i) - getv!.unsafe(vv,i)); return v end else if eqcar(u,'minus) then begin scalar vv; vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); compile!-angle!-expression(cadr u,vv); for i:=0:7 do putv!.unsafe(v,i,getv!.unsafe(v,i) - getv!.unsafe(vv,i)); return v; end else if eqcar(u,'times) then if numberp cadr u then begin scalar vv; vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); compile!-angle!-expression(caddr u,vv); for i:=0:7 do putv!.unsafe(v, i, cadr u*getv!.unsafe(vv, i) + getv!.unsafe(v,i)) end else begin scalar vv; vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); compile!-angle!-expression(cadr u,vv); for i:=0:7 do putv!.unsafe(v, i, caddr u * getv!.unsafe(vv, i) + getv!.unsafe(v,i)) end else nil; symbolic procedure simpfouriertimes(u); begin scalar z; z := car simpfourier1 car u; u := cdr u; a: if null u then return z ./ 1; z := fs!:times!:(car simpfourier1 car u,z); u := cdr u; go to a end; put('times, 'simpfour, 'simpfouriertimes); symbolic procedure simpfourierexpt(u); fs!:expt!:(car simpfourier1 car u, cadr u) . 1; put('expt, 'simpfour, 'simpfourierexpt); symbolic procedure simpfourierplus(u); begin scalar z; z := car simpfourier1 car u; u := cdr u; a: if null u then return z ./ 1; z := fs!:plus!:(car simpfourier1 car u,z); u := cdr u; go to a end; put('plus, 'simpfour, 'simpfourierplus); symbolic procedure simpfourierdifference(u); fs!:difference!:(car simpfourier1 car u, car simpfourier1 cadr u) ./ 1; put('difference, 'simpfour, 'simpfourierdifference); symbolic procedure simpfourierminus(u); fs!:negate!:(car simpfourier1 car u) . 1; put('minus, 'simpfour, 'simpfourierminus); symbolic procedure simpfourierquot(u); begin scalar v; v := simp!* cadr u; v := cdr v . car v; return fs!:times!:(car simpfourier1 car u, !*sq2fourier v) ./ 1 end; put('quotient, 'simpfour, 'simpfourierquot); symbolic procedure simphsin u; begin if not(length u = 1) then rerror(fourier,5,"Argument should be single expression"); return simpfouriersin list(u := prepsq simp!* car u) end; put('hsin, 'simpfn, 'simphsin); symbolic procedure simphcos u; begin if not(length u = 1) then rerror(fourier,6,"Argument should be single expression"); return simpfouriercos list(u := prepsq simp!* car u) end; put('hcos, 'simpfn, 'simphcos); endmodule; module hsub; %% Harmonic substitution: the CAMAL HSUB operation, as well as other %% substitutions. fluid '(!*trharm); switch trham; symbolic procedure hsub1(x,u,v,a,n); %% Substitute v+A for u in x to order n begin scalar ans, c, tmp, fs!:zero!-generated; %% fs!:zero!-generated := 0; ans := fs!:subang(x, u, v); % c := ensure!-fourier A; c := car a; if c then c := cdr c; a := c; if !*trham then << print "A"; if null a then print 0 else fs!:prin a >>; for i:=1:n do << if !*trham then << print "i="; print i >>; x := hdiff(x, u); if !*trham then << prin2!* "df(x,u,i)="; fs!:prin x; terpri!* t; prin2!* "A^i ="; fs!:prin c; terpri!* t >>; c := fs!:times(cdr !*sq2fourier (1 ./ i), c); if !*trham then << prin2!* "A^i/fact(i) ="; fs!:prin c; terpri!* t>>; tmp := fs!:times(fs!:subang(x, u, v), c); if !*trham then << prin2!* "f'(0)*A^i/fact i = "; fs!:prin tmp; terpri!* t>>; ans := fs!:plus(ans, tmp); if !*trham then << prin2!* "partial sum ="; fs!:prin ans; terpri!* t>>; if not(i=n) then c := fs!:times(c,a); >>; return ans end; symbolic procedure fs!:subang(x, u, v); if null x then nil else begin scalar vv, n; vv := mkvect 7; n := getv!.unsafe(fs!:angle x, u); for i:=0:7 do if i = u then putv!.unsafe(vv, i, n*getv!.unsafe(v,i)) else putv!.unsafe(vv, i, getv!.unsafe(fs!:angle x,i) + n*getv!.unsafe(v,i)); return fs!:plus(fs!:subang(fs!:next x, u, v), make!-term(fs!:fn x, vv, fs!:coeff x)); end; symbolic procedure fs!:sub(x,u); if null x then nil else begin scalar ans; ans := aeval prepsq fs!:coeff x; if not fixp ans then ans := subsq(cadr ans, u) else ans := fs!:coeff x; if eqcar(numr ans, '!:fs!:) then ans := cdar ans else ans := cdr !*sq2fourier ans; ans := fs!:times(make!-term(fs!:fn x, fs!:angle x, 1 ./ 1), ans); return fs!:plus(fs!:sub(fs!:next x, u), ans); end; symbolic procedure simphsub uu; begin scalar x, u, v, vv, a, n, dmode!*; dmode!* := '!:fs!:; if (length uu = 5) then << x := car uu; uu := cdr uu; u := car uu; uu := cdr uu; v := car uu; uu := cdr uu; a := car uu; uu := cdr uu; n := car uu >> else if (length uu = 3) then << x := car uu; uu := cdr uu; u := car uu; uu := cdr uu; v := car uu; uu := cdr uu; if not harmonicp u then << a := ( ((get('fourier, 'tag) . fs!:sub(cdar simp x, list(u . v))) ./ 1) ) where wtl!*=delasc(u,wtl!*); return a; >>; a := 0; n := 0 >>; if not harmonicp u then rerror(fourier, 7, "Not an angle in HSUB"); x := cdar simp x; if not angle!-expression!-p v then rerror(fourier, 8, "Not an angle expression in HSUB"); vv := mkvect 7; for i:=0:7 do putv!.unsafe(vv,i,0); compile!-angle!-expression(v, vv); a := simp!* a; n := simp!* n; if null car n then n := 0 ./ 1 else if not(fixp car n and cdr n = 1) then rerror(fourier, 9, "Non integer expansion in HSUB"); n := car n; return (get('fourier, 'tag) . hsub1(x,get(u,'fourier!-angle),vv,a,n)) ./ 1; end; put('hsub, 'simpfn, 'simphsub); endmodule; module hdiff; %% Harmonic differentiation and Integration. symbolic procedure hdiff(x, u); if null x then nil else fs!:plus(hdiff(fs!:next x,u), hdiffterm(x,u)); symbolic procedure hdiffterm(x, u); begin scalar n; n := getv!.unsafe(fs!:angle x, u); if n = 0 then return nil; n := multsq( n . 1, fs!:coeff x); if fs!:fn x = 'cos then return make!-term('sin, fs!:angle x, negsq n) else return make!-term('cos, fs!:angle x, n) end; symbolic procedure hdiff1(x, u); if null x then nil else begin scalar ans, aaa; ans := diffsq(fs!:coeff x, u); if ans then << aaa := mkvect 3; fs!:set!-coeff(aaa, ans); fs!:set!-fn(aaa, fs!:fn x); fs!:set!-angle(aaa,fs!:angle x); fs!:set!-next(aaa, hdiff1(fs!:next x, u)); return aaa >> else return hdiff1(fs!:next x, u) end; symbolic procedure simphdiff uu; begin scalar x, u; if not (length uu = 2) then rerror(fourier, 10, "Improper number of arguments to HDIFF"); x := car uu; uu := cdr uu; u := car uu; x := simp x; if not eqcar(car x, '!:fs!:) then x := !*sq2fourier x ./ 1; if not harmonicp u then return (get('fourier, 'tag) . hdiff1(cdar x, u)) ./ 1; x := hdiff(cdar x,get(u,'fourier!-angle)); if null x then return nil ./ 1; return (get('fourier, 'tag) . x) ./ 1 end; put('hdiff, 'simpfn, 'simphdiff); symbolic procedure hint(x, u); if null x then nil %% Bind fs!:zero!-generated ?? else fs!:plus(hint(fs!:next x,u), hintterm(x,u)); symbolic procedure hintterm(x, u); begin scalar n; n := getv!.unsafe(fs!:angle x, u); if n = 0 then return make!-term(fs!:fn x, fs!:angle x, fs!:coeff x); n := multsq( 1 ./ n, fs!:coeff x); if fs!:fn x = 'cos then return make!-term('sin, fs!:angle x, n) else return make!-term('cos, fs!:angle x, negsq n) end; symbolic procedure hint1(x , u); if null x then nil else begin scalar aaa; aaa := mkvect 3; fs!:set!-coeff(aaa, simpint list(prepsq fs!:coeff x, u)); fs!:set!-fn(aaa, fs!:fn x); fs!:set!-angle(aaa,fs!:angle x); fs!:set!-next(aaa, hint1(fs!:next x, u)); return aaa end; symbolic procedure simphint uu; begin scalar x, u; if not (length uu = 2) then rerror(fourier, 11, "Improper number of arguments to HINT"); x := car uu; uu := cdr uu; u := car uu; x := simp x; if not eqcar(car x, '!:fs!:) then x := !*sq2fourier x ./ 1; if not harmonicp u then return (get('fourier, 'tag) . hint1(cdar x, u)) ./ 1; x := hint(cdar x,get(u,'fourier!-angle)); if null x then return nil ./ 1; return (get('fourier, 'tag) . x) ./ 1 end; put('hint, 'simpfn, 'simphint); initdmode 'fourier; endmodule; end;