#=============================================================================
# Maple startup procedures
#=============================================================================
# Marc A. Murison
# Astronomical Applications Dept.
# U.S. Naval Observatory
# 3450 Massachusetts Ave., NW
# Washignton, DC 20392
# murison@aa.usno.navy.mil
# http://aa.usno.navy.mil/murison/
#=============================================================================

unprotect( norm, dot, mag, vec, mat, cross );

startup[norm] := proc(A) linalg[norm](A,2); end proc;

startup[dot] := proc(u,v) 
    local U, V;
    if type(evalm(u),'array'(1..posint,1..1)) then
        U := convert(evalm(u),vector);
    else
        U := copy(u);
    end if;
    if type(evalm(v),'array'(1..posint,1..1)) then
        V := convert(evalm(v),vector);
    else
        V := copy(v);
    end if;
    linalg[dotprod](U,V,'orthogonal'); 
end proc;

startup[mag] := proc(u) 
	norm(u);
end proc;

startup[vec] := proc() 
    if nargs=1 and type(args[1],list) then
        linalg[vector](args);
    elif nargs=1 and type(args[1],vector) then
        args[1];
    else
        linalg[vector]([args]); 
    end if;
end proc;

startup[mat] :=  proc () 
    local i,j,tmp,_H_;
    if nargs = 1 and type(args[1],vector) then 
        convert( args[1], 'matrix' );
    elif nargs = 1 and type(args[1],list) then 
        tmp := linalg[matrix]([seq([_H_.i],i = 1..nops(args[1]))]);
        for j from 1 to nops(args[1]) do
            tmp[j,1] := subs( _H_.j=args[1][j], tmp[j,1] );
        od;
        eval(tmp);
    elif nargs = 1 and type(args[1],'array'(2)) then
        eval(args[1]);
    else
        linalg[matrix]([seq([args[i]],i = 1 .. nargs)]);
    fi
end proc;

startup[matfunc] := proc( eqn::(matrix=matrix), func::procedure )
    map(func,lhs(eqn),args[3..nargs]) = map(func,rhs(eqn),args[3..nargs]);
end proc;

startup[cross] := proc(u::{vector,list,'array'(2)},v::{vector,list,'array'(2)}) 
    local U, V;
    if type(u,{list,'array'(2)}) then
        U := convert(evalm(u),vector);
    else
        U := copy(u);
    end if;
    if type(v,{list,'array'(2)}) then
        V := convert(evalm(v),vector);
    else
        V := copy(v);
    end if;
    if size(U)=2 then
        U := vec( U[1], U[2], 0 );
    end if;
    if size(V)=2 then
        V := vec( V[1], V[2], 0 );
    end if;
    linalg[crossprod](U,V); 
end proc;

protect( norm, dot, mag, vec, mat, cross );

#
# fix table indexing so that we can make assignments
# to table/array/set/sequence elements (see Heck, p. 312)
#
# WARNING (version V.4): this breaks fsolve()
#
#`index/newtable` := proc( ind, tab, val )
#   if not assigned(tab) then
#       tab := table();
#       tab[op(ind)] := op(val);
#   elif nops([eval(tab)]) > 1 and type(ind,[integer])  #expression sequence
#     and ind[1] >= 1 and ind[1] <= nops([eval(tab)]) then
#       tab := op( subsop(ind[1]=op(val),[eval(tab)]) );
#   elif type(eval(tab),set) and type(ind,[integer])
#     and ind[1] >= 1 and ind[1] <= nops(eval(tab)) then
#       tab := subsop( ind[1]=oop(val), eval(tab) );
#   else
#       ERROR(`: not a table/array/set/sequence`);
#   end if;
#end proc;


#------------------------------------------------------------------------------
# functions to read my procedure files
#------------------------------------------------------------------------------
startup[build_utils]  := proc() read cat(UTILS_PATH,"utils.p");         with(utils);       end proc;
startup[build_ode]    := proc() read cat(ODE_PATH,"diffeqn.p");         with(ode);         end proc;
startup[build_plot]   := proc() read cat(PLOTTING_PATH,"plotting.p");   with(plotting);    end proc;
startup[build_curves] := proc() read cat(CURVES_PATH,"curves.p");       with(curves);      end proc;
startup[build_surfs]  := proc() read cat(SURFACES_PATH,"surfaces3D.p"); with(surfaces3D);  end proc;
startup[build_poly]   := proc() read cat(POLY_PATH,"polynomials.p");    with(polynomials); end proc;
startup[build_start]  := proc() read cat(STARTUP_PATH,"startup.p");     with(startup);     end proc;

save( startup, cat(STARTUP_PATH,"startup.m") );