Astronomical Applications Department, U.S. Naval Observatory - Linear Least Squares (Page 7)

Astronomical Applications Department, U.S. Naval Observatory Linear Least Squares Page 7

# solve the normal equations
#---------------------------
debug_print(procname,"Solving the normal equations...",1);
_EnvAllSolutions := true;
map( allvalues, [solve( convert(eqs,set), convert(params,set) )] );
map( sortsols, %, params );
if nops(%)=1 then
solset := op(%);
else
solset := %;
debug_print(procname,"There are ".(nops(solset))." solutions",1);
fi;
debug_print(procname,"raw solution(s)",4,solset);

#---------------------------------------------------
# replace the denominators with another substitution
#---------------------------------------------------

#first define a workhorse procedure
replace_denominators := proc( expr, N, nmax, S, Delta )

local k, p, d;
global _subslist;

#---------------------------------
# determine the common denominator
#---------------------------------
#grab and filter denominators from solution parameter expressions
map( denom, { seq(op(rhs(expr[k])),k=1..nops(expr)) } );
map( (x)->if type(x,{integer,fraction}) then 0 else x fi, % );
d := % minus {0};
if nops(d) > 1 then
debug_print(procname,"denominator set",0,d);
ERROR("Unable to determine a unique denominator!");
elif d={} then
debug_print(procname,"no denominator for solution",2,expr);
RETURN(expr);
else
d := op(d);
fi;

#--------------------------------------
# make the substitution in the solution
#--------------------------------------
_subslist := [ Delta=eval(d,_subslist), op(_subslist) ];
collect( algsubs(d=Delta,expr), [Delta,seq(S[k],k=0..nmax),N], factor );
end;

#now make the substitutions
if type(solset,list(list)) then #multiple solutions
for k from 1 to nops(solset) do
subsop(k=replace_denominators(solset[k],N,nmax,S,Delta),solset);
solset := subs( Delta=delta[k], % );
_subslist := subs( Delta=delta[k], _subslist );
od;
else #only one solution
solset := replace_denominators(solset,N,nmax,S,Delta);

Page 7