############################################################################### 
#
#pertutils.mpl (c) 1998 by George Davies
#
#written for GRTensorII (c) 1994-98 by Peter Musgrave, Denis Pollney and Kayll Lake
#
#pertutils.mpl is provided free of charge on the condition that it is
#acknowledged that we accept no liability for software performance,
#continued maintenance, or damage to data. The author retains any and
#all rights to the software its documentation and distribution. 
#
###############################################################################
# linpert() linerarizes its argument wrt epsilon.
#
############################################################################### 

linpert := proc (gr_obj)

local new_obj, c0, c1:

readlib(coeftayl):

c0 := coeftayl(gr_obj, epsilon = 0, 0):
c1 := coeftayl(gr_obj, epsilon = 0, 1):

new_obj := c0 + c1*epsilon;

RETURN(new_obj):

end:

############################################################################### 
#
# quadpert() expands its argument to second order in epsilon.
#
############################################################################### 

quadpert := proc (gr_obj)

local new_obj, c0, c1, c2:

readlib(coeftayl):

c0 := coeftayl(gr_obj, epsilon = 0, 0):
c1 := coeftayl(gr_obj, epsilon = 0, 1):
c2 := coeftayl(gr_obj, epsilon = 0, 2):

new_obj := c0 + c1*epsilon + c2*epsilon^2;

RETURN(new_obj):

end:

###############################################################################
# termsimp() is designed to simplify an expression once the collect() routine
# has been applied to it. We will select out each term and simplify and
# replace it. This way maple will not undo the `collect' when it tries to do
# the simplification.
###############################################################################

termsimp := proc(terms)

local num_obs, new_terms, new_term, cntr:

num_obs := nops(terms):

new_terms := factor(simplify(op(1,terms))):

for cntr from 2 to num_obs do
   if (op(cntr,terms) <> `` and type(terms,`+`) ) then
      new_terms := new_terms + factor(simplify(op(cntr,terms)))
   elif (op(cntr,terms) <> `` and type(terms,`*`) ) then
      new_terms := new_terms*factor(simplify(op(cntr,terms)))
   fi:
od:

RETURN(new_terms):

end:

###############################################################################
# mcollect() is a collection routine that collects some compound objects.
# Here we can collect factors like U*V if the they appear in the same term.
###############################################################################

mcollect := proc(terms,list)

local new_terms, new_factor, this_term, tmp_terms, tmp_list,
   num_obs, num_in_list, list_prod, 
   cntr1, cntr2, mcount, c1, tmp:

num_obs := nops(terms):
num_in_list := nops(list):

new_terms := `start1`:
new_factor := `start2`:

tmp_terms := terms:

if not type(terms,`+`) then 
   ERROR(`Make sure the expression has no global factors.`)
fi:

for cntr1 from 1 to num_in_list do
   if ( op(0,op(1,op(1,list[cntr1]))) = 'diff' ) then
      tmp_terms := subs(list[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
od:

for cntr1 from 1 to num_in_list do
   if ( op(0,op(1,list[cntr1])) = 'diff' ) then
      tmp_terms := subs(list[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
od:

for cntr1 from 1 to num_in_list do
   if ( op(0,list[cntr1]) = 'diff' ) then
      tmp_terms := subs(list[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
od:

for cntr1 from 1 to num_in_list do
   if ( op(0,list[cntr1]) <> 'diff' ) then
      tmp_terms := subs(list[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
   if ( cntr1 = 1) then
      list_prod := tmp_list[1]:
   else
      list_prod := list_prod*tmp_list[cntr1]:
   fi:
od:

if (num_in_list = 1) then
   list_prod := list_prod*list_prod:
fi:

for cntr1 from 1 to num_obs do
      
   this_term := op(cntr1,tmp_terms):
   
   if (this_term <> ``) then

      mcount := 0:

      for cntr2 from 1 to num_in_list do

 	 if (match(this_term=c1*tmp_list[cntr2],tmp_list[cntr2],'tmp') and (num_in_list <> 1)) then

	    unassign('tmp'):

	    mcount := mcount + 1:

         elif (match(this_term=c1*tmp_list[cntr2]^2,tmp_list[cntr2],'tmp') and (num_in_list = 1)) then

            unassign('tmp'):

            mcount := 1:

            break:
	    
         fi:

      od:

      if (mcount = num_in_list) then 

         if (new_factor = `start2`) then
	    new_factor := this_term/list_prod:
         else 
	    new_factor := new_factor + this_term/list_prod:
	 fi:

      else

         if (new_terms = `start1`) then
            new_terms := this_term:
         else
            new_terms := new_terms + this_term:
         fi:

      fi:
      
   fi:

od:

if (new_factor = `start2`) then

   RETURN(terms):

else

   for cntr1 from 1 to num_in_list do

      new_terms := subs(tmp_list[cntr1]=list[cntr1],new_terms):
      list_prod := subs(tmp_list[cntr1]=list[cntr1],list_prod):
      new_factor := subs(tmp_list[cntr1]=list[cntr1],new_factor):

   od:

   if (new_terms = `start1`) then
      new_terms := factor(simplify(new_factor))*list_prod:
   else
      new_terms := factor(simplify(new_factor))*list_prod + new_terms:
   fi:

   RETURN(new_terms):

fi:

end:

###############################################################################
# hcollect() is a hierarcial collection routine that will collect wrt the 
# first supplied list and then collect the prefactors wrt the second list.
###############################################################################

hcollect := proc(terms,list1,list2)

local this_term, tmp_terms, new_terms, new_term, 
      num_obs, cntr1, cntr2, num_in_list1, c1, tmp,
      tmp_list:

num_in_list1 := nops(list1):

new_terms := `start1`:

tmp_terms := collect(terms,list1):

for cntr1 from 1 to num_in_list1 do
   if ( op(0,op(1,op(1,list1[cntr1]))) = 'diff' ) then
      tmp_terms := subs(list1[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
od:

for cntr1 from 1 to num_in_list1 do
   if ( op(0,op(1,list1[cntr1])) = 'diff' ) then
      tmp_terms := subs(list1[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
od:

for cntr1 from 1 to num_in_list1 do
   if ( op(0,list1[cntr1]) = 'diff' ) then
      tmp_terms := subs(list1[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
od:

for cntr1 from 1 to num_in_list1 do
   if ( op(0,list1[cntr1]) <> 'diff' ) then
      tmp_terms := subs(list1[cntr1]=tmp_list[cntr1],tmp_terms):
   fi:
od:

if type(tmp_terms,`*`) then
   num_obs := 1: 
else
   num_obs := nops(tmp_terms):
fi:

for cntr1 from 1 to num_obs do
   
   if (num_obs = 1) then
      this_term := tmp_terms:
   else
      this_term := op(cntr1,tmp_terms):
   fi:
     
   new_term := `start2`:

   if (this_term <> `` and type(this_term,`*`) ) then

      for cntr2 from 1 to nops(list1) do

         if match(this_term=c1*tmp_list[cntr2],tmp_list[cntr2],'tmp') then
	    
	    unassign('tmp'):

	    new_term := termsimp(collect(this_term/tmp_list[cntr2],list2)):
	    new_term := new_term*tmp_list[cntr2]:

	    break:

         fi:
       
      od:

   fi:
     
   if (new_term = `start2`) then
      new_term := this_term:
   fi:

   if (new_terms = `start1`) then
      new_terms := new_term:
   else
      new_terms := new_terms + new_term:
   fi:

od:

for cntr1 from 1 to num_in_list1 do
   new_terms := subs(tmp_list[cntr1]=list1[cntr1],new_terms):
od:

RETURN(new_terms):

end:

###############################################################################
# kfactor() pulls a user specifed factor out of a type `+` expression.
###############################################################################

kfactor := proc(terms,my_factor)

local num_obs, new_terms, new_term, cntr:

if not type(terms,`+`) then 
   ERROR(`Make sure the expression has no global factors.`)
fi:

num_obs := nops(terms):

new_terms := 0:

for cntr from 1 to num_obs do

   if (op(cntr,terms) <> ``) then

      new_terms := new_terms + factor(simplify(op(cntr,terms)/my_factor)):

   fi:
od:

new_terms := my_factor*new_terms:

RETURN(new_terms):

end: