Astronomical Applications Department, U.S. Naval Observatory Linear Least Squares Page 6
indets(value(expand(partials)),function), sum ),
list );
#degree() can't handle summation subscripts, so remove them
YTsums := subs( t=T, y=Y, sums );
_subslist := [];
nmax := 0;
debug_print(procname,"summation terms",4,sums);
debug_print(procname,"YT terms",5,YTsums);
Tcount := 0;
Scount := 0;
for k from 1 to nops(YTsums) do
if ispoly then
select(has,[op(YTsums[k])],T);
if %=[] then
n := 0;
else
n := degree( op(%), T );
fi;
nmax := max(n,nmax);
fi;
if has(YTsums[k],Y) then
if not ispoly then
n := Scount;
nmax := max(n,nmax);
Scount := Scount + 1;
fi;
_subslist := [op(_subslist),sums[k]='S'[n]];
else
if not ispoly then
n := Tcount;
Tcount := Tcount + 1;
fi;
_subslist := [op(_subslist),sums[k]='T'[n]];
fi;
od;
#---------------------------------------------------
# substitute for the summation terms in the partials
#---------------------------------------------------
partials := collect( eval( value(expand(partials)), _subslist ),
[seq(S[k],k=1..nmax),N], factor );
#put substitution list into form suitable for putting the sums back
_subslist := map( (x)->rhs(x)=lhs(x), _subslist );
#--------------------------
# form the normal equations
#--------------------------
eqs := collect( partials, [op(params),sum], factor );
if printlevel >= 1 then
debug_print(procname,"normal equations",1);
for p in eqs do
print(p);
od;
fi;
#---------------------------
Page 6
