
(* Last updated: August 4, 2007 at 00:55 at home in Christchurch *)
(* File: c0804-07.m *)
(* Code now works under Mathematica versions 4 and 5 *)

(* Bug fixes see: WH 08/04/2007 *)

(* Based in file : c060717h.m *)
(* Previously updated: 07/17/2006 by Hereman at office (Golden) at 15:00 *)

(* Provided the option of using the old algorithm or the new algorithm *)
(* to investigate the compatibility conditions. *)
(* See useoldalgorithm and usenewalgorithm *)

(* Rewrote the normal function to normalize the density. *)
(* Add newanalyzer to dynamically clean up the system. *)
(* New method of selecting the c[i] for the case were there are parameters *)
(* in the system. *)
(* The code now determines the listzeroc of the c[i] that must be zero, *)
(* computes listnonzeroc as the complement relative to the list unknowns, *)
(* and set these c[i] equal to one (one at the time). *)

(* Added a function to check is the next rho is independent from the *)
(* previously computed rho *) 
(* Keep track of all the independent rhos in listallrhos. *)
(* Similarly, keep track of independent fluxes in a listallfluxes *)

(* Changes marked by 07/03/2006 WH, also look for word debugflags *)
(* Purpose: Put print statements in the routines that analyze and solve *)
(* the system for the constants c[i] and find the compatibility conditions *)

(* ********************************************************************** *)

(* Based on condens.m from 06/10/2003 *)
(* 06/10/2003 Added a LogicalExpand to Reduce in solver for the c[i] *)
(* Replaced worklog.m by workcond.m *)
(* Always saves the algebraic system etc in workcond.m now *)
(* 02/23/2002 Hereman *)
   
(*****************************************************************************)
(*                                                                           *)
(*          *** M A T H E M A T I C A   P R O G R A M ***                    *)
(*                                                                           *)
(*      SYMBOLIC COMPUTATION of CONSERVED DENSITIES for SYSTEMS of           *)
(*                  NONLINEAR EVOLUTION EQUATIONS                            *)
(*                                                                           *)
(* program name: condens.m                                                   *)
(*                                                                           *)
(* purpose: computation of the conserved densities and corresponding flux    *)
(*          with possible compatibility conditions, verification of the      *)
(*          conservation law.                                                *)
(*                                                                           *)
(* input to condens.m : system of nonlinear evolution equations of any order,*)
(*                      any degree, polynomial type, in variables x and t,   *)
(*                      only constant parameters, functions in x and t are   *)
(*                      NOT allowed as parameters                            *)
(*                                                                           *)
(*                      u[i]_t = f(u[1],...,u[N],u[1]_{nx},...,u[N]_{nx})    *)
(*                      with i=1,...,N; and n=1,2,...                        *)
(*                                                                           *)
(* output : density and flux of desired rank (if it exists),                 *)
(*          and compatibility conditions on parameters, if applicable        *)
(*                                                                           *)
(* tested on : PCs with Windows XP and workstations with Linux               *)
(*                                                                           *)
(* language : Mathematica v. 5.0, backward compatible with v. 4.0            *)
(*                                                                           *)
(* authors: Unal Goktas and Willy Hereman                                    *)
(*          Department of Mathematical and Computer Sciences                 *)
(*          Colorado School of Mines                                         *)
(*          Golden, CO 80401-1887, USA                                       *)
(*                                                                           *)
(* Last updated: August 4, 2007                                              *)
(* Based on version 3.2: May 4, 1998                                         *)
(*                                                                           *)
(* Copyright 1998-2007                                                       *)
(*                                                                           *)
(*****************************************************************************)

(* WH 08/04/2007 print statement added *)
Print["Loading package condens.m of August 4, 2007."];

Clear["Global`*"]; 

(* WH 08/04/2007 set debug flags *)
(* 07/03/2006 WH debugflags *)
debugstripper = False;
debugmystripper = False;
debugroutine4 = False;
debugnewanalyze = False;
debugsolve = False;
debugnormal = False;
debugevaluate = False;
usenewalgorithm = True;
useoldalgorithm = False;

(* ------------------------------------------------------------------------- *)

(*****************************************************************************)
(* commentinter[]: prints a welcome message                                  *)
(*****************************************************************************)

commentinter[] := Block[{},
Print["*********************************************************"];
Print["        WELCOME TO THE MATHEMATICA PROGRAM               "];
Print["         by UNAL GOKTAS and WILLY HEREMAN                "];
Print["    FOR THE COMPUTATION OF CONSERVED DENSITIES.          "];
Print["      Last updated: August 4, 2007                       "];
Print["   Based on version 3.2 released on May 4, 1998          "];
Print["           Copyright 1998-2007                           "];
Print["*********************************************************"]
]; (* end commentinter *)

(* ------------------------------------------------------------------------- *)

(*****************************************************************************)
(* cls: clears the screen                                                    *)
(*****************************************************************************)

(* 
cls := Print["\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"];
*)
(*****************************************************************************)
(* printpage[ , ]: a subroutine for menu                                     *)
(*****************************************************************************)
 
printpage[n_,page_] := Module[{i,choice,lenpage},
cls; 
lenpage = Length[page[n]];
Print["  *** MENU INTERFACE ***  (page: ",n,")"];
Print["-------------------------------------------"];
For[i=1,i <= lenpage,i++,
   Print[Part[Part[page[n],i],1]]
   ];
Print[" nn) Next Page"];
Print[" tt) Your System"];
Print[" qq) Exit the Program"];
Print["------------------------------------------"];
choice = Input["ENTER YOUR CHOICE: "];
Return[choice] 
]; (* end printpage *)

(*****************************************************************************)
(* menu: creates the menu                                                    *)
(*****************************************************************************)

menu := Module[
{counterpage = 1,menulist,numpages,page,
choice1,control,choice2,choice3,lenmenulist,i},
menulist = {
{"  1) KdV Equation (d_kdv.m)"},
{"  2) Modified KdV Equation (d_mkdv.m)"},
{"  3) Fifth Order KdV Equation-parameterized (d_5kdv.m)"},
{"  4) Fifth Order KdV Equation-Lax Case (d_5lax.m)"},
{"  5) Fifth Order KdV Equation-SK Case (d_5sk.m)"},
{"  6) Fifth Order KdV Equation-KK Case (d_5kk.m)"},
{"  7) Fifth Order KdV Equation-Ito Case (d_5ito.m)"},
{"  8) Seventh Order KdV Equation-parameterized (d_7kdv.m)"},
{"  9) Seventh Order KdV Equation-Lax Case (d_7lax.m)"},
{" 10) Seventh Order KdV Equation-SK-Ito Case (d_7ski.m)"},
{" 11) Seventh Order KdV Equation KK Case (d_7kk.m)"},
{" 12) Schamel Equation (d_scham.m)"},
{" 13) Rosenau-Hyman Equation (d_roshy.m)"},
{" 14) Hirota-Satsuma System-parameterized (d_phrsat.m)"},
{" 15) Hirota-Satsuma System (d_hirsat.m)"},
{" 16) Ito System-parameterized (d_pito.m)"},
{" 17) Ito System (d_ito.m)"},
{" 18) NLS System (real and imaginary parts) (d_nls.m)"},
{" 19) NLS System (equation and conjugate) (d_nlsalt.m)"},
{" 20) DNLS System (d_dnls.m)"},
{" 21) MVDNLS System (d_mvdnls.m)"},
{" 22) Kaup System-parameterized (d_pkaup.m)"},
{" 23) Kaup System (d_kaup.m)"},
{" 24) Kaup-Broer System (d_broer.m)"},
{" 25) Drinfel'd-Sokolov System (d_soko.m)"},
{" 26) Dispersive Long Wave System (d_disper.m)"}, 
{" 27) Non-dispersive Long Wave System (d_nodisp.m)"}, 
{" 28) 3-Component KdV System (d_3ckdv.m)"},
{" 29) 2-Component Nonlinear Schrodinger Equation (d_2cnls.m)"},
{" 30) Boussinesq System (d_bous.m)"}  
}; (* closes menulist *)

lenmenulist = Length[menulist];
numpages = Ceiling[lenmenulist/10];
For[i = 1,i <= numpages,i++,
   page[i] = If[lenmenulist >= (i*10),
               menulist[[Table[k,{k,(i-1)*10+1,i*10}]]],
               menulist[[Table[k,{k,(i-1)*10+1,lenmenulist}]]]]
   ];

choice1 = printpage[counterpage,page];

control := (
Switch[choice1,
tt,Print["Make sure that you have prepared the data file for the system"];
   Print["you want to test (similar to the data files we supplied)."];
   choice2 = Input["If your file is ready, press 1, else 2: "];
   If[choice2 === 1,
    choice3 = InputString["Enter the name of your data file: "];
    Get[choice3],
    Print["Aborting the computations!"];
    Print["Prepare your data file first, then start over."];
    Abort[] 
     ], 
nn,If[counterpage < numpages,counterpage++;
     choice1 = printpage[counterpage,page]; control,
     counterpage = 1; choice1 = printpage[1,page]; control],
qq,Print["Aborting the computations!"];Abort[],
1,<<d_kdv.m,2,<<d_mkdv.m,3,<<d_5kdv.m,4,<<d_5lax.m,5,<<d_5sk.m,
6,<<d_5kk.m,7,<<d_5ito.m,8,<<d_7kdv.m,9,<<d_7lax.m,10,<<d_7ski.m,
11,<<d_7kk.m,12,<<d_scham.m,13,<<d_roshy.m,14,<<d_phrsat.m,15,<<d_hirsat.m,
16,<<d_pito.m,17,<<d_ito.m,18,<<d_nls.m,19,<<d_nlsalt.m,20,<<d_dnls.m,
21,<<d_mvdnls.m,22,<<d_pkaup.m,23,<<d_kaup.m,24,<<d_broer.m,25,<<d_soko.m,
26,<<d_disper.m,27,<<d_nodisp.m,28,<<d_3ckdv.m,29,<<d_2cnls.m,30,<<d_bous.m,_,
Print["Aborting the computations!"];Abort[]
]; (* closes Switch *)
); (* end control *)
control; 
If[Not[StringQ[myfile]],
  myfile = InputString["Enter the name of the output file: "]]
]; (* end menu *)

(*****************************************************************************)
(* EulerD[  ]: applies the Euler operator for variational derivatives.       *)
(*             Written by Yu He in 1992, simplified by W. Hereman            *)
(*****************************************************************************)

EulerD[f_, (y_)[x_, r___], w:{x_, r___}] := 
  Module[{Dfuncs, Dtimes, dummyfunc}, 
    Dfuncs = Union[Cases[{f}, Derivative[__][y][__], Infinity]]; 
    Dtimes = (Head[Head[#1]] & ) /@ Dfuncs /. Derivative -> List; 
    Expand[D[f, y[x, r]] + (ReleaseHold[Thread[dummyfunc[(D[f, 
                 #1] & ) /@ Dfuncs, 
             (Hold[Apply[Sequence, #1]] & ) /@ 
              (Thread[{w, #1}] & ) /@ Dtimes]]] /. dummyfunc -> D) . 
        ((-1)^#1 & ) /@ (Apply[Plus, #1] & ) /@ Dtimes] 
    ]; (* end EulerD *)
 
EulerD[f_, v:{(y_)[x_, r___], ___}, w:{x_, r___}] := 
   (EulerD[f, #1, w] & ) /@ v /; 
   If[Apply[And, (MatchQ[#1, _[Apply[Sequence, w]]] & ) /@ v], 
   True, Message[EulerD::argx, w]];

EulerD[f_, (y_)[x_], x_] := EulerD[f, y[x], {x}];
EulerD[f_, v:{(y_)[x_], ___}, x_] := EulerD[f, v, {x}];

(*****************************************************************************)
(* pdeform[]: takes an expression and prints it with subscripted derivatives *)
(*            and removes the functional dependence (i.e. f[x,t] -> f)       *)
(*            Written by Tracy Otto & Tony Miller (CSM-1995).                *)
(*            Modified by U. Goktas.                                         *)
(*****************************************************************************)

pdeform[expres_] := expres /. {
        Derivative[n__][u[k_]][x__] :>
                SequenceForm[u, Subscript[k],Subscript[","],
                        Subscript[
                         SequenceForm @@ Flatten[Table[#[[1]], {#[[2]]}]& /@
                                        Transpose[{{x}, {n}}]]]],
        u[n_][x__] :> SequenceForm[u,Subscript[n]] 
}; (* end pdeform *)

(*****************************************************************************)
(* length[]: a modified length function                                      *)
(*****************************************************************************)

length[a_] := Module[{l},
If[Head[a] === Plus,l = Length[a],
  If[Head[a] === Times,l = Length[a],l = 1 ] 
  ]; Return[l] 
  ]; (* end length *)

(*****************************************************************************)
(* split[]: if free coefficients c[i] remain in the density, the density     *)
(*          is further split in independent pieces corresponding to each     *)
(*          coefficient c[i]                                                 *)
(*****************************************************************************)

split[expr_] := Module[{listofc,lenlistofc,ci,rhoi,rest}, 
    rest=Expand[expr];
    If[FreeQ[rest,c[__]], 
    Print["The density has no free coefficients."], 
    listofc=Union[Cases[rest,c[__],12]];
    lenlistofc=Length[listofc];
    Print["There is/are ", lenlistofc, " free coefficient(s) in the density."];
    Print["These free coefficients are: ", listofc, "."];
    Print["Splitting the density in independent pieces."];
     Do[
       ci=Part[listofc,i];
       rhoi=Expand[Coefficient[rest,ci]]; 
       rest=Expand[rest-ci*rhoi];
       Print["Part of the density with coefficient ", ci, ":"];
       Print[rhoi],{i,1,lenlistofc}]; (* end do *)
     Print["Part of the density that has no free coefficient:"];
     Print[rest]
    ] (* end if *) 
]; (* end module *)

(* new normal function *)
(*****************************************************************************)
(* normal[]: normalizes the density according to the coefficient of the      *)
(*           first term in rho[x,t], returns normalized forms of rho[x,t]    *)
(*           and j[x,t]. Does not divide by the coefficient if the           *)
(*           coefficient is not a number. Returns a list with rho in the     *)
(*           first place, and j in the second place. To be used on the list  *)
(*           that comes out of stripper[rho[x,t]]                            *)
(*****************************************************************************)

normal[mylist_List] := Module[
{expandrho,firsttermrho,firstpartmylist,normalizationcoef,
normalizationcoefcipart}, 
expandrho = Expand[rho[x,t]];
If[debugnormal, 
   Print["At pt. NORMAL01 in normal, expanded form of rho, expandrho: "];
   Print[expandrho]
  ];
If[Head[expandrho] === Plus, (* then 1 *)
   firsttermrho = Part[expandrho,1];
   If[debugnormal, 
    Print["At pt. NORMAL02 in normal, rho is a sum, firsttermrho: "];
     Print[firsttermrho]
   ], (* else *)
   firsttermrho = expandrho;
   If[debugnormal, 
     Print["At pt. NORMAL03 in normal, rho is not a sum, firsttermrho: "];
     Print[firsttermrho]
     ]
  ]; (* end if 1 *)
If[Head[firsttermrho] === Times, (* then 2 *)
   firstpartmylist=Part[mylist,1];
   If[debugnormal, 
     Print["At pt. NORMAL04 in normal, Part[mylist,1], firstpartmylist:"];
     Print[firstpartmylist]
     ];
   normalizationcoef = firsttermrho/firstpartmylist; 
(* WH 08/04/2007 taking the part of normalizationcoef that is independent *)
(* of c[i] *)
   If[debugnormal, 
     Print["At pt. NORMAL05 in normal, normalizationcoef:"];
     Print[normalizationcoef]
     ];
   normalizationcoefcipart = Cases[normalizationcoef,c[_]];
   If[debugnormal, 
      Print["At pt. NORMAL06 in normal, normalizationcoefcipart:"];
      Print[normalizationcoefcipart]
     ];
   If[
     (normalizationcoefcipart =!= {} && Length[normalizationcoefcipart] >= 1),
     (* then *)
     normalizationcoef = 
     Coefficient[normalizationcoef,normalizationcoefcipart[[1]]]
     ];
   If[debugnormal, 
     Print["At pt. NORMAL07 in normal, updated normalizationcoef:"];
     Print[normalizationcoef]
     ];
   If[debugnormal, 
   Print["At pt. NORMAL08 in normal, divide firsttermrho by firstpartmylist,"<>
         " normalizationcoef:"];
     Print[normalizationcoef]
     ];
rho[x,t] = Expand[Factor[rho[x,t]/normalizationcoef]]; 
If[debugnormal, 
   Print["At pt. NORMAL09 in normal, after division by normalizationcoef,"<>
         " rho[x,t]: "];
   Print[rho[x,t]]
  ];
j[x,t] = Expand[Factor[j[x,t]/normalizationcoef]];
If[debugnormal, 
   Print["At pt. NORMAL10 in normal, after division by normalizationcoef,"<>
         " j[x,t]: "];
   Print[j[x,t]]
  ];
If[Not[NumberQ[normalizationcoef]],  
  Print["WARNING!"];
  Print["Both rho[x,t] and j[x,t] were divided by: ",normalizationcoef,"."]
  ] (* no else 2, so leave rho and j unscaled *)
]; (* end if 2 *)
rhojlist = {rho[x,t],j[x,t]};
If[debugnormal, 
   Print["At pt. NORMAL11 in normal, rhojlist :"];
   Print[rhojlist]
   ];
Return[rhojlist]
]; (* end normal *)

(* old normal function *)

(* begin taken out 
(*****************************************************************************)
(* normal[]: normalizes the density according to the coefficient of the      *)
(*           first term in rho[x,t], returns normalized forms of rho[x,t]    *)
(*           and j[x,t]. Does not divide by the coefficient if the           *)
(*           coefficient is not a number. Returns a list with rho in the     *)
(*           first place, and j in the second place. To be used on the list  *)
(*           that comes out of stripper[rho[x,t]]                            *)
(*****************************************************************************)

normal[mylist_List] := Module[
{leadingcoef,facexpr,denrho,rhojlist}, 
(* Cancel common denominators first, keep track of these with parameters *)
facexpr = Factor[rho[x,t]];
If[debugnormal, 
   Print["At pt. NORMAL01 in normal, factored form of rho, facexpr: "];
   Print[facexpr]
  ];
denrho = Denominator[facexpr];
If[debugnormal, 
   Print["At pt. NORMAL02 in normal, denominator of rho, denrho: "];
   Print[denrho]
  ];
If[Not[NumberQ[denrho]], denrho = Factor[denrho]];
(* if 1 *)
If[Not[FreeQ[Head[denrho],Times]],
(* if 2 *)
   If[NumberQ[Part[denrho,1]] && Part[denrho,1] =!= 0, 
     denrho = denrho/Part[denrho,1];
   If[debugnormal, 
      Print["At pt. NORMAL03 in normal, after division by part[denhro,1],"<>
            " denrho: "];
      Print[denrho]
      ]
] (* end if 2 *)
]; (* end if 1 *)
If[Not[NumberQ[denrho]],  
  Print["WARNING!"];
  Print["Both rho[x,t] and j[x,t] were multiplied by: ",denrho,"."]; 
  rho[x,t] = Expand[Factor[rho[x,t]*denrho]]; 
  If[debugnormal, 
     Print["At pt. NORMAL04 in normal, after multiplication by denrho,"<>
           " rho[x,t]: "];
     Print[rho[x,t]]
    ];
  j[x,t] = Expand[Factor[j[x,t]*denrho]]];
  If[debugnormal, 
     Print["At pt. NORMAL05 in normal, after multiplication by denrho,"<>
           " j[x,t]: "];
     Print[j[x,t]]
    ];
leadingcoef = Coefficient[rho[x,t],Part[mylist,1]];
  If[debugnormal, 
     Print["At pt. NORMAL06 in normal, leadingcoef :"];
     Print[leadingcoef]
    ];
If[NumberQ[leadingcoef] && leadingcoef =!= 0,
  rho[x,t] = Expand[rho[x,t]/leadingcoef];
  If[debugnormal, 
     Print["At pt. NORMAL07 in normal, after division by leadingcoef,"<>
           " rho[x,t]: "];
     Print[rho[x,t]]
    ];
  j[x,t] = Expand[j[x,t]/leadingcoef];
  If[debugnormal, 
     Print["At pt. NORMAL08 in normal, after division by leadingcoef,"<>
           " j[x,t]: "];
     Print[j[x,t]]
    ]
  ]; 
  rhojlist = {rho[x,t],j[x,t]};
  If[debugnormal, 
     Print["At pt. NORMAL09 in normal, rhojlist :"];
     Print[rhojlist]
    ];
Return[rhojlist]
]; (* end normal *)
 end taken out *)

(* 01/23/07 WH renamed variable givenexpr to givenexprfac *)
(* 07/11/2006 WH mystripper *)
(*********************************************************************** *)
(* mystripper[]: cancels the numerical factors and members of nzsymbols, *)
(* and anything that has no C in each term of the argument               *)
(*********************************************************************** *)

freeofc[x_]:=FreeQ[x,c];
notfreeofc[x_]:=Not[FreeQ[x,c]];

mystripper[givenexpr_] := Module[
{faclist,factorsfreeofc,factorswithc,factorstokeep,factorstocancel,
factorscancelled,givenexprfac},
givenexprfac = Factor[givenexpr];
  If[debugmystripper,
    Print["At pt. MYS1, factored form givenexprfac: ", givenexprfac]
    ]; 
If[Head[givenexpr] === Times, (* then *)
faclist = FactorList[givenexprfac];
  If[debugmystripper,
    Print["At pt. MYS2, faclist: ", faclist]
    ]; 
factorsfreeofc = Select[faclist,freeofc[#]&];
  If[debugmystripper,
    Print["At pt. MYS3, factorsfreeofc: ", factorsfreeofc]
    ]; 
factorswithc = Complement[faclist, factorsfreeofc];
  If[debugmystripper,
    Print["At pt. MYS4, factorswithc: ", factorswithc]
    ]; 
factorstocancel = Select[factorsfreeofc,(Not[Head[Part[#,1]] === Plus])&];
  If[debugmystripper,
    Print["At pt. MYS5, factorstocancel: ", factorstocancel]
    ]; 
factorstokeep = Select[factorsfreeofc,(Head[Part[#,1]] === Plus)&];
  If[debugmystripper,
    Print["At pt. MYS6, factorstokeep: ", factorstokeep]
    ]; 
strippedequation = 
  Part[Flatten[factorswithc],1]*Product[
      factorstokeep[[k]][[1]]^factorstokeep[[k]][[2]],
      {k,1,Length[factorstokeep]}];
  If[debugmystripper,
    Print["At pt. MYS7, strippedequation: ", strippedequation]
    ]; 
factorscancelled = Factor[givenexpr/strippedequation];
  If[debugmystripper,
    Print["At pt. MYS8, factorscancelled: ", factorscancelled]
    ], (* else head is not times *)
strippedequation = givenexpr
]; (* end if for head testing *)
Return[strippedequation]
]; (* end of mystripper *)

(*****************************************************************************)
(* stripper[]: cancels the numerical factors of each term of the argument    *)
(*             and puts these into a list after cancellation                 *)
(*****************************************************************************)

stripper[givenexpr_] := Module[
{lenexpr1,list = {},expr1,expr2,expr3,expr3s,lenexpr3,i,k,s},
expr1 = Expand[givenexpr];
  If[debugstripper,
    Print["At pt. STR1, expand givenexpr: ", expr1]
    ]; 
If[Head[expr1] === Plus,
   lenexpr1 = Length[expr1],
   lenexpr1 = 1
  ];
(* start for 1 *)
For[k = 1,k <= lenexpr1,k++,
   If[lenexpr1 == 1,
     expr2 = expr1,
     expr2 = Part[expr1,k]
     ];
   expr2 = expr2 /. c[i_] :> 1;
   If[debugstripper,
     Print["At pt. STR2, for k=",k,", setting c's to one, expr2: "];
     Print[expr2]
     ]; 
   expr3 = FactorList[expr2];
   If[debugstripper,
     Print["At pt. STR3, for k=",k,", factorlist of expr2, expr3: "];
     Print[expr3]
     ]; 
   lenexpr3 = Length[expr3];
   If[debugstripper,
     Print["At pt. STR4, for k=",k,", lenexpr3: ", lenexpr3]
     ]; 
   (* start for 2 *)
   For[s = 1,s <= lenexpr3,s++,
      expr3s = Part[expr3,s];
      If[debugstripper,
        Print["At pt. STR5, for k=",k," and part s=",s," of expr3, expr3s:"];
        Print[expr3s]
        ]; 
      If[FreeQ[expr3s,u],
        If[FreeQ[weightpars,Part[expr3s,1]],
          expr2 = expr2/Part[expr3s,1]^Part[expr3s,2];
          If[debugstripper,
            Print["At pt. STR6, for k=",k,", after division, expr2: ", expr2]
            ]
          ],
        s = lenexpr3+1;
        ];
      ]; (* end for 2 *)
   list = Append[list,expr2];
   If[debugstripper,
     Print["At pt. STR7, for k=",k,", updated list: "];
     Print[list]
     ]
   ]; (* end for 1 *)
Return[list]
]; (* end stripper *)

(*****************************************************************************)
(* highestpair[]: a subroutine for highest order function                     *)
(*****************************************************************************)

highestpair[pairlist_] := Module[
{length,k,highestlist},
If[pairlist != {},
  length = Length[pairlist];
  highestlist = {Part[pairlist,1]};
  If[length >= 2,
    For[k = 2,k <= length,k++,
      If[Part[pairlist,k][[1]] > Part[highestlist,Length[highestlist]][[1]],
       highestlist = {Part[pairlist,k]},
       If[Part[pairlist,k][[1]] === Part[highestlist,Length[highestlist]][[1]],
	 highestlist = Union[Append[highestlist,pairlist[[k]]]]
         ];
        ];
      ];
   ];
 Return[highestlist],
 If[debug,
   Print["The expression (given to the function that computes the"]; 
   Print["highest order) does not have any derivatives!"]];
 Return[{{0,0}}];
 ]  
]; (* end highestpair *)
 
(*****************************************************************************)
(* highestorder[]: returns the highest order (or the term with the highest   *)
(*                 order) of the argument                                    *)
(*****************************************************************************)
 
highestorder[gexpr_,debug11_,debug12_,debug21_,debug22_,debug31_,debug32_] := 
Module[
{length1,i,newexpr1,length2,newexpr6,temporder1,temporder2,newexpr5,k,
newexpr2,newexpr3,temporder3,newexpr4,orderlist1 ={},orderlist2 = {},
orderlist3 = {},highestlist1,highestlist2,highestlist3,expr},
expr = Expand[gexpr];
If[Not[FreeQ[expr,Derivative]], 
  length1 = length[expr];
  temporder1 = {}; temporder2 = {}; temporder3 = {}; newexpr5 = 1;
  For[i = 1,i <= length1,i++,
     If[Not[FreeQ[expr,Plus]],temporder1 = {}; temporder2 = {};
            temporder3 = {};newexpr5 = 1];
     If[length1 != 1,newexpr1 = expr[[i]],newexpr1 = expr]; 
       If[Not[FreeQ[newexpr1,Derivative]],  
	 length2 = length[newexpr1];
	 For[k = 1,k <= length2,k++,
	    If[length2 != 1,newexpr2 = newexpr1[[k]],newexpr2 = newexpr1];
	    If[Not[FreeQ[newexpr2,Derivative]],
	      If[FreeQ[newexpr2,Power],
		newexpr3 = FullForm[newexpr2];
                newexpr6 = newexpr3[[1]][[0]][[0]];
		temporder1 = Append[temporder1,newexpr6[[1]]];
		temporder2 = Append[temporder2,newexpr6[[2]]];
		temporder3 = Append[temporder3,newexpr6[[1]]+newexpr6[[2]]];
		newexpr4 = newexpr2;
		newexpr5 = newexpr5*newexpr2, 
		newexpr3 = FullForm[newexpr2];
                newexpr6 = newexpr3[[1]][[1]][[0]][[0]];  
		temporder1 = Append[temporder1,newexpr6[[1]]]; 
		temporder2 = Append[temporder2,newexpr6[[2]]];
		temporder3 = Append[temporder3,newexpr6[[1]]+newexpr6[[2]]];
		newexpr4 = newexpr3[[1]];
		newexpr5 = newexpr5*newexpr3[[1]];
		];  
	      ];  
	    ];  
         temporder1 = Max[temporder1];
         temporder2 = Max[temporder2];
         temporder3 = Max[temporder3];
	 orderlist1 = Append[orderlist1,{temporder1,newexpr4}];
	 orderlist2 = Append[orderlist2,{temporder2,newexpr4}];
	 orderlist3 = Append[orderlist3,{temporder3,newexpr5}];
	 ];  
     ];        
  ];                  
Which[
debug11,highestlist1 = highestpair[orderlist1];
  If[debug12,Return[Table[highestlist1[[k]][[2]],{k,1,Length[highestlist1]}]],
     Return[highestlist1[[1]][[1]]]], 
debug21,highestlist2 = highestpair[orderlist2];
  If[debug22,Return[Table[highestlist2[[k]][[2]],{k,1,Length[highestlist2]}]],
      Return[highestlist2[[1]][[1]]]],
debug31,highestlist3 = highestpair[orderlist3];
  If[debug32,Return[Table[highestlist3[[k]][[2]],{k,1,Length[highestlist3]}]],
      Return[highestlist3[[1]][[1]]]] 
 ] (* end which *)
]; (* end highestorder *)

(*****************************************************************************)
(* setuniformrank[]: a subroutine for the scaling function                   *)
(*****************************************************************************)

setuniformrank[list_] := Module[
{i,lengthlist,syslist={}},
lengthlist = Length[list];
For[i = 1,i < lengthlist,i++,
   syslist = Append[syslist,Part[list,i] == Part[list,lengthlist]]; 
   ];
Return[syslist] 
]; (* end setuniformrank *)

(*****************************************************************************)
(* scaling[]: determines the scaling properties, if the equations are        *)
(*            incompatible it will print the appropriate messages            *)
(*****************************************************************************)

scaling[eqlist_List] := Module[
{pointsym,expr,list,syslist,msyslist = {},i,j,k,scalesol,
lenmsyslist,tempscalesol = {},tempmsyslist,trouble,troubleleft,
troubleright,posleft,posright,troublelist},
If[weightpars =!= {},
  pointsym = Union[Table[weightpars[[i]] -> E^weight[weightpars[[i]]],
                     {i,1,Length[weightpars]}],
       {u[i_Integer][x,t] :> E^weightu[i],
        D[u[i_Integer][x,t],{x,j_Integer}] :> E^(weightu[i]+j),
        D[u[i_Integer][x,t],{t,j_Integer}] :> E^(weightu[i]+j*weight[d/dt]),
        x -> E, t -> E^weight[d/dt]}],
  pointsym = {u[i_Integer][x,t] :> E^weightu[i],
        D[u[i_Integer][x,t],{x,j_Integer}] :> E^(weightu[i]+j),
        D[u[i_Integer][x,t],{t,j_Integer}] :> E^(weightu[i]+j*weight[d/dt]),
        x -> E, t -> E^weight[d/dt]};
  ];
For[i = 1, i <= noeqs, i++,
   expr[i] = Part[eqlist,i];
If[debugstripper, 
   Print["At pt. SCALING100, calling stripper in scaling function!"]
   ];
   list[i] = stripper[expr[i]];
   list[i] = list[i] /. pointsym;
   list[i] = PowerExpand[Log[list[i]]];
   syslist[i] = setuniformrank[list[i]];
   msyslist = Union[msyslist,syslist[i]];
   ];
scalesol = Flatten[Solve[msyslist]];
If[MemberQ[msyslist,False],
  Print["Fatal Error! Weights are incorrect."];
  Print["Check your choice(s) for the weights in the data file."];
  Print["Aborting the computations!"];
  CloseLog[]; Abort[],
  If[scalesol === {} && msyslist =!= {True},
    Print["In the given system there is at least one equation with terms"];
    Print["of unequal rank. Scaling properties can not be determined for"];
    Print["this system. The program will try to find the conflict, and,"];
    Print["if successful, provide suggestions to help resolve the conflict."];
    lenmsyslist = Length[msyslist];
    i = 1;
    While[tempscalesol === {} && i <= lenmsyslist,
         tempmsyslist = Complement[msyslist,{Part[msyslist,i]}];
         tempscalesol = Flatten[Solve[tempmsyslist]];
         i++;
         ];
    i--;
    If[tempscalesol =!= {},
      trouble = Part[msyslist,i];
      troubleleft = Part[trouble,1];
      troubleright = Part[trouble,2];
      For[j = 1, j <= noeqs, j++,
         posleft = Flatten[Position[list[j],troubleleft,{1}]];
         posright = Flatten[Position[list[j],troubleright,{1}]];
         If[posleft =!= {} && posright =!= {},
           troublelist = Union[
             Table[Part[Expand[eq[j][x,t]],posright[[k]]],
                     {k,1,Length[posright]}],
             Table[Part[Expand[eq[j][x,t]],posleft[[k]]],
                     {k,1,Length[posleft]}]
                ];
           Print["The terms"];
           Print[troublelist];
           Print["in equation ",j," are incompatible. Try to introduce an"];
           Print["auxiliary parameter with weight as coefficient of one of"];
           Print["these terms. Aborting the computations! "];
           ];
         ],
    Print["Try to introduce auxiliary parameters with weight as coefficients"];
    Print["into the system. Aborting the computations!"];
    ];
    CloseLog[];
    Abort[],
    Return[scalesol];
    ];
  ] 
]; (* end scaling *)

(*****************************************************************************)
(* subroutine4[]: a subroutine for the function evaluate                     *)
(*****************************************************************************)

(* 07/03/2006 WH added print statements to function subroutine4 *)
(* 07/15/2006 WH checked for presence of c[i] first, to avoid redundant work *)
subroutine4[givenexpr_] := Module[
{expr1,nexpr1,freerule={},part1,lennexpr1,i},
(* 07/03/2006 WH replaced expand by factor since if was factor in next line *)
(* expr1 = Expand[givenexpr]; *)
If[debugroutine4, 
   Print["At pt. SUB01, expresssion for rho, givenexpr: "];
   Print[givenexpr]
  ];
expr1 = Factor[givenexpr]; 
If[debugroutine4, 
   Print["At pt. SUB02, factored from of givenexpr (for rho), expr1: "];
   Print[expr1]
  ];
If[ Cases[expr1,_c,Infinity] === {}, (* then *)
  freerule={};
If[debugroutine4, 
   Print["At pt. SUB03, no c[i] were present, hence, freerule: "];
   Print[freerule]
  ], (* else continue with computation of freerule *)
nexpr1 = Numerator[expr1];
If[debugroutine4, 
   Print["At pt. SUB04, numerator of expr1, nexpr1: "];
   Print[nexpr1]
  ];
lennexpr1 = length[nexpr1];
If[debugroutine4, 
   Print["At pt. SUB05, length of nexpr1, lennexpr1: "];
   Print[lennexpr1]
  ];
For[i = 1,i<= lennexpr1,i++,
   part1 = Part[nexpr1,i];
If[debugroutine4, 
   Print["At pt. SUB06, first part of nexpr1, part1: "];
   Print[part1]
  ];
   If[MemberQ[unknownlist,part1],
     freerule = Union[freerule,{part1 -> 1}] 
     ];
If[debugroutine4, 
   Print["At pt. SUB07, freerule: "];
   Print[freerule]
  ]
   ]; (* end for loop *)
If[freerule =!= {},
  Print["Caution! There is a common (free) factor in the density."];
  Print["Setting ",freerule,"."];
  ]; 
If[debugroutine4, 
   Print["At pt. SUB08, final selected, freerule: "];
   Print[freerule]
  ]
]; (* end if *)
Return[freerule]
]; (* end subroutine4 *)

(*****************************************************************************)
(* subroutine5[]: a subroutine for the function constructformrho             *)
(*****************************************************************************)

subroutine5[nodims_,givenlist_] := Module[
{m,boundary,scale,var,len,tempwei,temppar,list,pair,ctm = 0,k,s},
list[0] = {};
For[m = 1,m <= nodims,m++,
   pair = Part[givenlist,m];
   var = Part[pair,1];
   scale = Part[pair,2];
   If[scale =!= 0,
     If[m === 1,
       boundary = Floor[rhorank/scale];
       list[m] = Table[{i*scale,var^i},{i,0,boundary}];
       ctm = m,
       len = Length[list[m-1]];
       list[m] = {};
       For[k = 1,k <= len,k++,
          tempwei = Part[Part[list[m-1],k],1];
          temppar = Part[Part[list[m-1],k],2];
          boundary = Floor[(rhorank-tempwei)/scale];
          For[s = 0,s <= boundary,s++,
             list[m] = Union[list[m],{{tempwei+s*scale,temppar*var^s}}];
             ctm = m;
             ];
          ];
       ];
     ];
   ];
Return[list[ctm]] 
]; (* end subroutine5 *)

(*****************************************************************************)
(* constructformrho[]: returns the form of rho                               *)
(*****************************************************************************)
 
constructformrho[nodims_,varscalelist_] := Module[
{list,len,tempwei,temppar,inttest,formrholist = {},triviallist = {},
difflist1,difflist2,lendifflist1,i,j}, 
Print["Starting the construction of the form of the density."];
list = subroutine5[nodims,varscalelist];
len = Length[list];
For[i = 1,i <= len,i++,
   tempwei = Part[Part[list,i],1];
   temppar = Part[Part[list,i],2];
   inttest = rhorank-tempwei;
   If[IntegerQ[inttest] && Not[FreeQ[temppar,u]],
     If[inttest === 0, 
       formrholist = Union[formrholist,{temppar}],
     If[debugstripper, 
        Print["At pt. CONSTRFORMRHO200, calling stripper in constructformrho!"]
        ];
       difflist1 = stripper[D[temppar,{x,inttest-1}]];
       lendifflist1 = Length[difflist1];
       For[j = 1, j <= lendifflist1,j++,
        If[debugstripper, 
        Print["At pt. CONSTRFORMRHO300, calling stripper in constructformrho!"]
        ];
          difflist2 = stripper[D[Part[difflist1,j],x]];
          triviallist = Union[triviallist,{Last[difflist2]}];
          formrholist = Union[formrholist,Complement[difflist2,triviallist]]; 
          ];
       ];
     ];
   ];

lenformrho = Length[formrholist];
formrho[x,t] = 0;
For[i = 1,i <= lenformrho,i++,
   formrho[x,t] = formrho[x,t]+c[i]*Part[formrholist,i];
   ];
Print["For RANK = ",rhorank,", this is the form of density rho[x,t]: "];
Print[pdeform[Expand[formrho[x,t]]]];
If[formrho[x,t] === 0,
  Print["The only density is the trivial density rho[x,t] = 0."];
  Print["Aborting the computations! "]; 
  CloseLog[]; Abort[] ] 
]; (* end constructformrho *)               

(*****************************************************************************)
(* evaluate[]: computes the density and flux, with verification              *)
(*****************************************************************************)

evaluate[expr1_,expr2_,solc_,rules_] := Module[
{test,factortest,newrhot,newformj,freerule,rhojlist,printflag,
leadingtermofrho}, 
rho[x,t] = expr1 /. solc;
rho[x,t] = Expand[rho[x,t]];
newrhot[x,t] = expr2 /. Flatten[{rules,solc}]; 
newrhot[x,t] = Expand[newrhot[x,t]];
newformj[x,t] = (-1)*Integrate[newrhot[x,t],x];
j[x,t] = Expand[newformj[x,t]];
(* start if 100 *)
If[FreeQ[j[x,t],Integrate] && rho[x,t] =!= 0, (* then if 100 *)
  If[rho[x,t] =!= 0,freerule = subroutine4[rho[x,t]],freerule = {}];
  rho[x,t] = rho[x,t] /. freerule;
  j[x,t] = j[x,t] /. freerule;
  test = Expand[D[rho[x,t],t]+D[j[x,t],x] /. rules];
  factortest = Expand[test];
  If[factortest =!= 0, factortest = Factor[factortest]];
  If[factortest =!= 0, (* if factortest not zero, then *)
    Print["Automatic verification of the conservation law FAILED!"],
    (* else factortest was zero *)
    If[debugstripper, 
      Print["At pt. EVALUATE400, calling stripper in evaluate function!"]
      ];
  leadingtermofrho = stripper[Part[Expand[rho[x,t]],1]]; 
(* first application of normal function *)
(* WH 08/04/2007 normal function may return c[i] in denominators *)
  rhojlist = normal[leadingtermofrho];
  rho[x,t] = Part[rhojlist,1];
  rho[x,t] = Factor[rho[x,t]];
  rho[x,t] = Expand[rho[x,t]]; 
(* 07/14/2006 WH capturing the rho in list *)
(* Extra test to verify if the new rho is a multiple of first *)
(* computed rho *)
(* TO BE DONE *)
(* Extra test to verify if the new rho is independent from all previously *)
(* computed rhos *)
If[listallrhos =!= {}, (* if x *)
   ratiorhos = Factor[rho[x,t]/listallrhos[[1]]]; 
   If[debugevaluate,
     Print["At pt. EVAL1, ratiorhos = "];
     Print[ratiorhos]
    ], (* else x *)
   ratiorhos = u
  ]; (* end if x *)
  If[rho[x,t] =!= 0 && Not[FreeQ[ratiorhos,u]], (* then *)
(* WH 08/04/2007 to avoid duplicates, added a Union *)
     listallrhos = Append[listallrhos,rho[x,t]];
     listallrhos = Union[listallrhos];
     printflag = True, (* else *)
     printflag = False
    ];
  If[debugevaluate,
    Print["At pt. EVAL2, listallrhos = "];
    Print[listallrhos]
    ];
  j[x,t] = Part[rhojlist,2];
  j[x,t] = Factor[j[x,t]];
  j[x,t] = Expand[j[x,t]];
(* 07/14/2006 WH capturing the fluxes in list *)
  If[j[x,t] =!= 0 && Not[FreeQ[ratiorhos,u]], 
    listallfluxes = Append[listallfluxes,j[x,t]];
(* WH 08/04/2007 to avoid duplicates, added a Union *)
    listallfluxes = Union[listallfluxes];
     printflag = True, (* else *)
     printflag = False
    ];
  If[debugevaluate,
    Print["At pt. EVAL3, listallfluxes = "];
    Print[listallfluxes]
    ];
If[printflag === True, (* if printflag is true, then *)
  Print["*******************************************************"];
  Print["The normalized density rho[x,t] is:"];
  Print[pdeform[rho[x,t]]];
  split[pdeform[rho[x,t]]];
  If[Length[j[x,t]] <= 20 || debug,
    Print["*******************************************************"];
    Print["The corresponding flux j[x,t] is:"];
    Print[pdeform[j[x,t]]], (* else length *)
    Print["*******************************************************"];
    Print["Flux j[x,t] has been computed. It has ",Length[j[x,t]]," terms."];
    Print["Since there are more than 20 terms, j[x,t] will not be shown."];
    Print["To see the flux at the end, type: pdeform[j[x,t]] or j[x,t]"]
    ];
  Print["*******************************************************"];
  Print["Result of explicit verification: (rho_t + J_x) = ",factortest];
  Print["*******************************************************"], 
  (* else when printflag is false *)
  Print["The rho[x,t] and j[x,t] are linearly dependent on the ones"<>
        " computed above!"]
  ], (* end if printflag *)
  (* else if 100 *)
  Print["*******************************************************"];
  Print["There is NO DENSITY of this form!"]; 
  Print["*******************************************************"]
  ] (* end if factortest *)
  ] (* end if 100 *)
]; (* end evaluate *) 

(*****************************************************************************)
(* constructeqlist[]: forms the system of equations for coefficients c[i]    *)
(*****************************************************************************)

constructeqlist[expr_] := Module[
{expr1,rules1,rules2,eqlist = {},tempexpr,expr2,templist,lentemplist,s,
tempcoef,parttemplist,expr3},
rules1 = D[u[i_][x,t],{x,j_}] :> 0;
rules2 = u[i_][x,t] :> 0;
expr1 = Expand[expr];
tempexpr = expr1 //. rules1;
tempexpr = tempexpr //. rules2;
eqlist = Union[eqlist,{tempexpr}]; 
expr2 = expr1-tempexpr;
If[expr2 =!= 0,
   If[debugstripper, 
      Print["At pt. CONSTREQLIST500, calling stripper in constructeqlist!"]
     ];
  templist = Union[stripper[expr2]];
  lentemplist = Length[templist];
  expr3 = expr2;
  For[s = 1,s <= lentemplist,s++,
     parttemplist = Part[templist,s];
     tempcoef = Coefficient[Expand[expr2],parttemplist];
     tempcoef = tempcoef //. rules1;
     tempcoef = tempcoef //. rules2;
     If[debug,expr3 = expr3 - (tempcoef*parttemplist)];
     eqlist = Union[eqlist,{tempcoef}];
     ];
  If[debug,
    expr3 = Factor[expr3];
    Print["Remainder of the expression during the construction of"];
    Print["the system for the coefficients c[i]: ",expr3]];
  If[debug && expr3 =!= 0,
    Print["Error in the construction of the system for the c[i]!"];
    Print["Aborting the computations!"];
    CloseLog[]; Abort[] ];
  ]; (* end if *)
Return[eqlist] 
]; (* end constructeqlist *)

(*****************************************************************************)
(* mysimplify[]: simplifies the system of equations for the c[i]             *)
(*****************************************************************************)
 
mysimplify[list_] := Module[
{newlist,lennewlist,k,fterm,simplelist = {}},
Print["Starting the simplification of the system."];
newlist = Complement[list,{0}];
lennewlist = Length[newlist];
For[k = 1 ,k <= lennewlist,k++,
   fterm = Factor[newlist[[k]]];
   If[Part[fterm,0] === Times && NumberQ[Part[fterm,1]],  
     fterm = fterm/Part[fterm,1] ];
   If[Intersection[Expand[simplelist],Expand[{fterm*(-1)}]] === {},
     simplelist = Union[simplelist,{fterm}]];
   ]; (* end for *)
Return[simplelist] 
]; (* end mysimplify *)

(*****************************************************************************)
(* newanalyzer[]: determines the coefficients c[i] that must be zero,        *)
(* Dynamically simplifies system and makes list of c[i] that must be zero    *)
(* used in the search for compatibility conditions                           *)
(*****************************************************************************)
(* analyzelist will be list of two pieces, *)
(* first piece: list of simplified lhs of equations *)
(* second piece: list of the c[i] that must be zero *)

newanalyzer[doublelist_List] := 
Module[{currentsystem,currentzeroc,strippedsystem,updatedzeroc,analyzelist,
updatedzerocrules,updatedsystem}, 
Print["Starting the repeated analysis and simplification of the system."];
currentsystem = Part[doublelist,1];
If[debugnewanalyze, 
   Print["At pt. ANN01 in newanalyzer, currentsystem: "];
   Print[currentsystem]
  ];
currentzeroc = Part[doublelist,2];
If[debugnewanalyze, 
   Print["At pt. ANN02 in newanalyzer, currentzeroc: "];
   Print[currentzeroc]
  ];
strippedsystem = Union[Map[mystripper,currentsystem]];
If[debugnewanalyze, 
   Print["At pt. ANN03 in newanalyzer, after mapping mystripper,"<>
         " strippedsystem: "];
   Print[strippedsystem]
  ];
updatedzeroc = Union[Join[currentzeroc,Cases[strippedsystem,_c,1]]];
If[debugnewanalyze, 
   Print["At pt. ANN04 in newanalyzer, updatedzeroc: "];
   Print[updatedzeroc]
  ];
updatedzerocrules = Union[Map[Rule[#,0]&, updatedzeroc]];
If[debugnewanalyze, 
   Print["At pt. ANN05 in newanalyzer, updatedzerocrules: "];
   Print[updatedzerocrules]
  ];
updatedsystem = Complement[strippedsystem /. updatedzerocrules,{0}];
If[debugnewanalyze, 
   Print["At pt. ANN06 in newanalyzer, after application of"<>
         " updatedzerocrules, updatedsystem: "];
   Print[updatedsystem]
  ];
analyzelist = {updatedsystem, updatedzeroc};
If[debugnewanalyze, 
   Print["At pt. ANN07 in newanalyzer, analyzelist: "];
   Print[analyzelist]
  ];
If[debugnewanalyze, 
  Print["At pt. ANN08 in newanalyzer, "];
  Print["  !!!!!!!!! END OF ONE RUN OF NEWANALYZER !!!!!!!!!!  "];
  ];
Return[analyzelist]
]; (* end newanalyzer *)

(*****************************************************************************)
(* picklhs[]: makes a system, setting LHS == 0                               *)
(*****************************************************************************)

picklhs[list_,k_] := Part[list,k] == 0;

(*****************************************************************************)
(* coefmat[]: creates the coefficient matrix of the system for the c[i]      *)
(*****************************************************************************)

coefmat[system_List,unknowns_List] := Module[
{i,mat = {},lensys}, 
If[system === {} || unknowns === {}, Print["Fatal Error"],
lensys = Length[system];
For[i = 1,i <= lensys,i++, 
   mat = Append[mat,Table[Coefficient[Expand[Part[Part[system,i],1]],
   Part[unknowns,k]],{k,1,Length[unknowns]}]]]];
Return[mat] 
]; (* end coefmat *)

(*****************************************************************************)
(* main[]: solves the system for the c[i], computes the compatibility        *)
(*         conditions (if applicable)                                        *)
(*****************************************************************************)

main[utlist_] := Module[
{rules = {},trules,ressubst,formj,tformj,dtformj,comcond,k,comcondfactab,
myeqlist,originalmaineqlist,maineqlist,maineqlistlhs,lengthpar,syscondlist1={},
syscond,inputval,inputrule,sol1,lengthsol1,nrules1,par,newmaineqlist,
solc,nrules2,nsolc,funclist,j,seqlist,lenseqlist,analyzelist,analyzeliststart,
analyzelistend,listzeroc,listzerocrules,listnonzeroc,myruletrick,inputpart,
inputlist,i,iii,complexanalysis,comcondfac},
rhseq[i_] := D[u[i][x,t],t] /. Part[utlist,i];
rule[i_,n_] := D[D[u[i][x,t],{x,n}],{t,1}] -> D[rhseq[i],{x,n}];
For[i = 1,i <= noeqs,i++,
   trules[i] = Table[rule[i,n],{n,0,highestorderrho}];
   rules = Union[rules,trules[i]]
   ];
If[debug,
  Print["These are the replacement rules from the given system: ",rules]];
(* Computation of j *)
rhot[x,t] = D[formrho[x,t],t];
ressubst = rhot[x,t] /. rules;
ressubst = Expand[-ressubst]; 
If[debug,Print["The Euler operator will be applied to: ",
  pdeform[ressubst]]]; 
funclist = Table[u[i][x,t],{i,1,noeqs}];
Print["Applying the Euler operator (the variational derivative)."];
formj = EulerD[ressubst,funclist,{x,t}];
Print["Finished with variational differentiation."];
tformj = formj;
(* begin if point 1 *)
If[Length[Cases[tformj,x_0]] === noeqs, (* then point 1 *)
  Print["End of the computations."];
  evaluate[formrho[x,t],rhot[x,t],{},rules], (* else point 1 *)
  dtformj = tformj; 
  dtformj = Expand[dtformj];
  Print["Starting the construction of the system for the c[i]."];
  myeqlist = {};
  For[j = 1,j <= noeqs,j ++,
     If[Part[dtformj,j] =!= 0,
       myeqlist = Union[myeqlist,constructeqlist[Part[dtformj,j]]]];
     ];
 seqlist = mysimplify[myeqlist];
 lenseqlist = Length[seqlist];
 originalmaineqlist = Flatten[Table[picklhs[seqlist,i],{i,1,lenseqlist}]];
 system = MapAll[Factor,originalmaineqlist];
If[ (* start if 0 *)
(Not[debug] && Length[originalmaineqlist]>50) && (parameters==={} && 
 weightpars==={}), 
 Print["Since there are no parameters, and the linear system"];
 Print["for the c[i] has ",Length[maineqlist],
       " equations, which is more than 50,"];
 Print["the system will not be shown."];
 Print["To see the system at the end of the computations, type: system"],
 Print["This is the system for the c[i]:"]; 
 Print[system];
 Print["List of unknown coefficients c[i]:"];
 Print[unknownlist]
 ]; (* end if 0 *)

(* 02/23/2002 Hereman: workcond.m instead of worklog.m *)
(* complexanalysis > 20 replaced by complexanalysis > 0 to force it *)
(* to always save the information in a file workcond.m *)

   lengthpar = Length[parameters];
   Print["The unsimplified system for the coefficients c[i] has ",
          lengthpar," parameter(s)."];
   complexanalysis=lengthpar+Length[originalmaineqlist]+Length[unknownlist];
   If[complexanalysis > 0, 
   coefmatrix = coefmat[system,unknownlist];
   Print["The system of PDEs, rank and form of rho, the system for"];
   Print["c[i], its coefficient matrix, and the lists of unknowns and"]; 
   Print["parameters, will all be saved in the file workcond.m for"];
   Print["further analysis. For analysis use the Mathematica functions"]; 
   Print["Reduce, Solve, etc. Load the file workcond.m with the command"]; 
   Print["<<workcond.m ."];
   Save["workcond.m",myfile,eqlist,rhorank,formrho,system,coefmatrix,
   unknownlist,parameters];
   Clear[eqlist,system,coefmatrix];
   ]; (* end if complexityanalysis *)

(* 07/13/2006 WH completely reworked *)
(* start of the simplification and analysis of the system *)
If[parameters =!= {}, (* then 8 start for parameter test *)
 analyzeliststart = {seqlist,{}};
 If[debugnewanalyze, 
   Print["At pt. ANN00, entering newanalyzer, with analyzeliststart : "];
   Print[analyzeliststart]
   ];
   analyzelistend = FixedPoint[newanalyzer,analyzeliststart,250];
 If[debugnewanalyze, 
   Print["At pt. ANN100, leaving newanalyzer, with analyzelistend : "];
   Print[analyzelistend]
   ];
maineqlistlhs = Part[analyzelistend,1];
 If[debugnewanalyze, 
   Print["At pt. ANN200, leaving newanalyzer, with maineqlistlhs : "];
   Print[maineqlistlhs]
   ];
maineqlist = Table[Part[maineqlistlhs,k] == 0, {k,1,Length[maineqlistlhs]} ]; 
 If[debugnewanalyze, 
   Print["At pt. ANN300, leaving newanalyzer, with maineqlist : "];
   Print[maineqlist]
   ];
listzeroc = Part[analyzelistend,2];
 If[debugnewanalyze, 
   Print["At pt. ANN400, leaving newanalyzer, with listzeroc : "];
   Print[listzeroc]
   ];
listzerocrules = Union[Map[Rule[#,0]&,listzeroc]];
 If[debugnewanalyze, 
   Print["At pt. ANN500, leaving newanalyzer, with listzerocrules : "];
   Print[listzerocrules]
   ];
listnonzeroc = Complement[unknownlist,listzeroc];
 If[debugnewanalyze, 
   Print["At pt. ANN600, leaving newanalyzer, with listnonzeroc : "];
   Print[listnonzeroc]
   ];
(* 07/03/2006 WH added print statements to solver *)
 For[i=1,i <= lengthpar,i++,
    syscondlist1 = Append[syscondlist1,parameters[[i]] != 0];
    ];
   syscond = Apply[And,syscondlist1];
   Print["Starting the search for compatibility conditions."];
   Print["The software will assume that ",syscond,"."]; 
   Print["Cases where one or more of these parameters are zero must be"<>
         " run separately!"];
   Print["The program will set the nonzero coefficients equal to 1,"<> 
         " one at the time."]; 
(* remove all c[i] that must be zero *)
   inputlist = Complement[unknownlist,listzeroc];
 If[debugsolve, 
   Print["At pt. SO1, c[i] to be set to one (one at the time), inputlist : "];
   Print[inputlist]
   ];
   If[inputlist === {}, (* then 6 no else 6 *)
     Print["All of the coefficients c[i] in the density have to vanish."];
     Print["The search for compatibility conditions ends."]; 
     ];
   While[inputlist =!= {},
    inputpart = Part[inputlist,1];
    inputval = inputpart == 1;
    inputlist = Complement[inputlist,{inputpart}];
    inputrule = ToRules[inputval];
    Print["-------------------------------------------------------"];
    Print["* Setting ",inputrule[[1]],":"];
    Print["Computation of the compatibility conditions."];
 If[debugsolve, 
   Print["At pt. SO2a, entering Eliminate after applying inputrule,"<>
         " maineqlist: "];
   Print[maineqlist /. inputrule]
   ];
 If[debugsolve, 
   Print["At pt. SO2b, entering Eliminate, listnonzeroc : "];
   Print[listnonzeroc]
   ];
    comcond = Eliminate[maineqlist /. inputrule, listnonzeroc];
 If[debugsolve, 
   Print["At pt. SO3a, leaving Eliminate, comcond : "];
   Print[comcond]
   ];
(* 07/03/2006 WH move terms to left hand side of equation, then factor! *)
myruletrick = {lhs__ == rhs__ -> lhs - rhs == 0};
comcond = comcond /. myruletrick;
 If[debugsolve, 
   Print["At pt. SO3b, moving all terms to left hand side, comcond : "];
   Print[comcond]
   ];
    comcondfac = MapAll[Factor,comcond];
 If[debugsolve, 
   Print["At pt. SO4, after mapall factor on comcond, comcondfac : "];
   Print[comcondfac]
   ];
    (* Remove possible duplicates *)
    If[Head[comcondfac] === Or, 
      comcondfactab = Table[Part[comcondfac,k], {k,1,Length[comcondfac]}];
      comcondfac = Union[comcondfactab]; 
      If[debugsolve, 
        Print["At pt. SO5, after removing possible duplicates, comcondfac : "];
        Print[comcondfac]
        ]
      ]; (* end if head *)
    Print["This is the compatibility condition:"];
    Print[comcondfac];
(* 06/10/2004 LogicalExpand added to make it work under version Math 5.0 *)
 If[debugsolve, 
   Print["At pt. SO6, entering Reduce, parameters: "];
   Print[parameters]
   ];
 If[debugsolve, 
   Print["At pt. SO7, entering Reduce, comcond && syscond: "];
   Print[comcond && syscond]
   ];
 If[debugsolve, 
   Print["At pt. SO7bis, entering Reduce, comcond: "];
   Print[comcond]
   ];
(* HERE HERE HERE *)
(* Bug fix WH 08/04/2007 different version to make code work under *)
(* earlier versions of Mathematica *)
(* For Mathematica version 5.0 and higher *)
If[$VersionNumber >= 5.,
    flip[equ_] := equ /. { Equal[a_,b_] -> Equal[b,a] };
    sol1 = LogicalExpand[Reduce[comcond && syscond, parameters]];
    sol2 = LogicalExpand[Reduce[comcond && syscond, Reverse[parameters]]]; 
    sol2 = Map[flip, sol2];
    rootfindersol2 = Cases[sol2, Power[_,1/2], Infinity];
    If[rootfindersol2 === {}, finalsol = sol2];
    rootfindersol1 = Cases[sol1, Power[_,1/2], Infinity];
    If[rootfindersol1 === {}, finalsol = sol1];
    sol1 = finalsol
  ];

(* For Mathematica versions lower than 5.0 *)
If[$VersionNumber < 5.,
    flip[equ_] := equ /. { Equal[a_,b_] -> Equal[b,a] };
    sol1 = Reduce[comcond && syscond, parameters];
    sol2 = Reduce[comcond && syscond, Reverse[parameters]]; 
    sol2 = Map[flip, sol2];
    rootfindersol2 = Cases[sol2, Power[_,1/2], Infinity];
    If[rootfindersol2 === {}, finalsol = sol2];
    rootfindersol1 = Cases[sol1, Power[_,1/2], Infinity];
    If[rootfindersol1 === {}, finalsol = sol1];
    sol1 = finalsol
  ];
(* END HERE HERE HERE *)

If[debugsolve, 
   Print["At pt. SO8, leaving Reduce,  sol1: "];
   Print[sol1]
   ];
    If[comcondfac === True,   
      Print["The compatibility condition is satisfied without constraints"];
      Print["on the parameters."] 
      ];
    If[sol1 =!= False, (* then 4 *)
      If[!FreeQ[sol1,Or], (* then 5 case where sol1 has or in it *)
        lengthsol1 = Length[sol1];  
        (* for iii *)
        For[iii = 1,iii <= lengthsol1,iii++,
           If[debugsolve, 
              Print["At pt. SO9, taking Part[sol1,iii] : "];
              Print[Part[sol1,iii]]
              ];
           nrules1 = ToRules[Part[sol1,iii]];
           If[debugsolve, 
              Print["At pt. SO10, nrules1 : "];
              Print[nrules1]
              ];
           If[comcondfac=!=True, 
             For[k=1,k <= lengthpar,k++,
           If[debugsolve, 
              Print["At pt. SO11, before applying nrules1, parameters[[k]]: "];
              Print[parameters[[k]]]
              ];
                par[k]=parameters[[k]] /. nrules1;   
           If[debugsolve, 
              Print["At pt. SO12, after applying nrules1, par[k]: "];
              Print[par[k]]
              ];
                Print["For ",parameters[[k]]," = ",par[k] ];
                ] 
             ];
           If[debugsolve, 
              Print["At pt. SO13a, before applying nrules1, maineqlist: "];
              Print[maineqlist]
              ];
           newmaineqlist = maineqlist /. nrules1;
           If[debugsolve, 
              Print["At pt. SO13b, after applying nrules1, newmaineqlist: "];
              Print[newmaineqlist]
              ];
           newmaineqlist = newmaineqlist /. inputrule;
           If[debugsolve, 
             Print["At pt. SO14, after applying inputrule, going into solver"<>
                   " newmaineqlist: "];
             Print[newmaineqlist]
             ];
           If[debugsolve, 
              Print["At pt. SO15, unknowns going into solver, listnonzeroc: "];
              Print[listnonzeroc]
              ];
           solc = Flatten[Solve[newmaineqlist, listnonzeroc]];
           If[debugsolve, 
              Print["At pt. SO16, after solving for c[i], solc: "];
              Print[solc]
              ];
           solc = Union[solc,inputrule];
           If[debugsolve, 
              Print["At pt. SO17a, after union with inputrule, solc: "];
              Print[solc]
              ]; 
           solc = Union[solc,listzerocrules];
           If[debugsolve, 
              Print["At pt. SO17b, after union with listzerocrules, solc: "];
              Print[solc]
              ]; 
           Print["The solution of the system is:"];
           Print[solc];
           If[debugsolve, 
              Print["At pt. SO18, rules: "];
              Print[rules]
              ];
           If[debugsolve, 
              Print["At pt. SO19, nrules1: "];
              Print[nrules1]
              ];
           nrules2 = rules /. nrules1;
           If[debugsolve, 
              Print["At pt. SO20, rules /. nrules1 gives nrules2: "];
              Print[nrules2]
              ];
           evaluate[formrho[x,t],rhot[x,t],solc,nrules2];
If[useoldalgorithm, (* if useoldalgorithm then case 1 *)
(* begin taken out 1 *)
(* case with parameters, and sol1 has or in it, *)
(* begin old routine for compatibility analysis *)
           If[debugsolve, 
              Print["At pt. OLDALGO SO21a, before INTERSECTION,"<> 
                    " original inputlist: "];
              Print[inputlist]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21b, solc: "];
              Print[solc]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21c, Position[solc,_->0]: "];
              Print[Position[solc,_->0]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21d, Flatten[Position[solc,_->0]]: "];
              Print[Flatten[Position[solc,_->0]]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21e,"<>
                    " solc[[ Flatten[Position[solc,_->0]] ]]:"];
              Print[solc[[Flatten[Position[solc,_->0]]]]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21f, after rule"<>
                    " x:(c[a_]->0)->c[a] gives:"];
              Print[solc[[Flatten[Position[solc,_->0]]]] /.
                    x:(c[a_]->0)->c[a] ]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21g, First[#]& /@ solc :"];
              Print[First[#]& /@ solc]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21h,"<>
                    " Complement[inputlist,First[#]& /@ solc] :"];
              Print[Complement[inputlist,First[#]& /@ solc]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO21i, union of modified solc and"<>
                    " complement :"];
              Print[Union[
                   solc[[ Flatten[Position[solc,_->0]] ]] /. 
                          x:(c[a_]->0)->c[a],
                   Complement[inputlist,First[#]& /@ solc] 
                   ]
              ]];
           inputlist = 
             Intersection[
               inputlist,
               Union[
                   solc[[ Flatten[Position[solc,_->0]] ]] /. 
                          x:(c[a_]->0)->c[a],
                   Complement[inputlist,First[#]& /@ solc]
               ] 
             ]; (* end intersection *)
           If[debugsolve, 
              Print["At pt. OLDALGO SO22, after INTERSECTION,"<>
                    " updated inputlist: "];
              Print[inputlist]
              ], (* else for usenewalgorithm case 1 *)
(* case with parameters, and sol1 has or in it, *)
(* end old routine for compatibility analysis *)
(* end taken out 1 *)
(* case with parameters, and sol1 has or in it, *)
(* begin new routine for compatibility analysis *)
         If[inputlist =!= {}, (* then *)
           inputlist = Drop[inputlist,1];
           If[debugsolve, 
              Print["At pt. NEWALGO SO23, case with parameters, after"<>
                    " removing the first element, updated inputlist: "];
              Print[inputlist]
              ], (* else *)
           inputlist = {}
           ] (* end if *)
]; (* end if useoldalgorithm case 1 *)
(* case with parameters, and sol1 has or in it, *)
(* end new routine for compatibility analyis *)
           Clear[nrules1,nrules2,solc,newmaineqlist];
           ], (* for iii *) (* else 5 where sol1 has no or in it *)
        nrules1 = ToRules[sol1];
        If[debugsolve, 
           Print["At pt. OLDALGO or NEWALGO SO24, nrules1: "];
           Print[nrules1]
           ];
        If[comcondfac=!=True, (* then 6 *)
           For[k = 1,k <= lengthpar,k++,
           If[debugsolve, 
              Print["At pt. OLDALGO or NEWALGO SO25,"<>
                    " before applying nrules1, parameters[[k]]: "];
              Print[parameters[[k]]]
              ];
              par[k]=parameters[[k]] /. nrules1;
           If[debugsolve, 
              Print["At pt. OLDALGO or NEWALGO SO26,"<>
                    " after applying nrules1, par[k]: "];
              Print[par[k]]
              ];
              Print["For ",parameters[[k]]," = ", par[k]];
              ]
           ]; (* end if 6 *)
           If[debugsolve, 
              Print["At pt. OLDALGO or NEWALGO SO27,"<>
                    " before applying nrules1, maineqlist: "];
              Print[maineqlist]
              ];
        newmaineqlist = maineqlist /. nrules1;
           If[debugsolve, 
              Print["At pt. OLDALGO or NEWALGO SO28,"<>
                    " after applying nrules1, newmaineqlist: "];
              Print[newmaineqlist]
              ];
        newmaineqlist = newmaineqlist /. inputrule;
           If[debugsolve, 
             Print["At pt. OLDALGO or NEWALGO SO29,"<>
                   " after applying inputrule, going into solver"<>
                   " newmaineqlist: "];
             Print[newmaineqlist]
             ];
           If[debugsolve, 
              Print["At pt. OLDALGO or NEWALGO SO30,"<>
                    " unknowns going into solver, listnonzeroc: "];
              Print[listnonzeroc]
              ];
        solc = Flatten[Solve[newmaineqlist, listnonzeroc]];
           If[debugsolve, 
              Print["At pt. OLDALGO or NEWALGO SO31,"<>
                    " after solving for c[i], solc: "];
              Print[solc]
              ];
           solc = Union[solc,inputrule];
              If[debugsolve, 
                Print["At pt. OLDALGO or NEWALGO SO32a,"<>
                      " after union with inputrule, solc: "];
                Print[solc]
              ];
           solc = Union[solc,listzerocrules];
           If[debugsolve, 
              Print["At pt. OLDALGO or NEWALGO SO32b,"<>
                    " after union with listzerocrules, solc: "];
              Print[solc]
              ]; 
        Print["The solution of the system is "];
        Print[solc];
        If[debugsolve, 
           Print["At pt. OLDALGO or NEWALGO SO33,"<>
                 " system replacement substitution rules: "];
           Print[rules]
           ];
        If[debugsolve, 
           Print["At pt. OLDALGO or NEWALGO SO34, nrules1: "];
           Print[nrules1]
           ];
        nrules2 = rules /. nrules1;
        If[debugsolve, 
           Print["At pt. OLDALGO or NEWALGO SO35,"<>
                 " rules /. nrules1 gives nrules2: "];
           Print[nrules2]
           ];
        evaluate[formrho[x,t],rhot[x,t],solc,nrules2];
If[useoldalgorithm, (* if useoldalgorithm then case 2 *)
(* begin taken out 2 *)
(* case with parameters, and sol1 has no or in it, *)
(* begin old routine for compatibility analysis *)
           If[debugsolve, 
              Print["At pt. OLDALGO SO36a, before INTERSECTION,"<>
                    " original inputlist: "];
              Print[inputlist]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36b, solc: "];
              Print[solc]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36c, Position[solc,_->0]: "];
              Print[Position[solc,_->0]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36d, Flatten[Position[solc,_->0]]: "];
              Print[Flatten[Position[solc,_->0]]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36e,"<>
                    " solc[[ Flatten[Position[solc,_->0]] ]]:"];
              Print[solc[[Flatten[Position[solc,_->0]]]]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36f, after rule"<>
                    " x:(c[a_]->0)->c[a] gives:"];
              Print[solc[[Flatten[Position[solc,_->0]]]] /.
                    x:(c[a_]->0)->c[a] ]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36g, First[#]& /@ solc :"];
              Print["All c[i] that were fixed after solving:"];
              Print[First[#]& /@ solc]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36h,"<>
                    " Complement[inputlist,First[#]& /@ solc] :"];
              Print["All c[i] that remained undetermined after solving:"];
              Print[Complement[inputlist,First[#]& /@ solc]]
              ];
           If[debugsolve, 
              Print["At pt. OLDALGO SO36i, union of modified solc and"<>
                    " complement :"];
              Print["All c[i] that are undetermined or zero:"];
              Print[Union[
                   solc[[ Flatten[Position[solc,_->0]] ]] /. 
                          x:(c[a_]->0)->c[a],
                   Complement[inputlist,First[#]& /@ solc] 
                   ]
               ]
              ];
        inputlist = Intersection[inputlist,
             Union[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a],
                   Complement[inputlist,First[#]& /@ solc]] 
             ]; 
           If[debugsolve, 
              Print["At pt. OLDALGO SO37, after INTERSECION,"<>
                    " updated inputlist: "];
              Print["All c[i] that were zero after solving, plus all "];
              Print["undetermined c[i], minus c[i] that must always be zero:"];
              Print[inputlist];
              Print["Investigate what happens when these c[i] are set to one"]
              ], (* else if usenewalgorithm case 2 *)
(* case with parameters, and sol1 had no or in it, *)
(* end old routine for compatibility analysis *)
(* end taken out 2 *)
(* case with parameters, and sol1 has no or in it, *)
(* begin new routine for compatibility analysis *)
         If[inputlist =!= {}, (* then *)
           inputlist = Drop[inputlist,1];
           If[debugsolve, 
              Print["At pt. NEWALGO SO38, case without parameters,"<>
                    " after removing the first element, updated inputlist: "];
              Print[inputlist]
              ], (* else *)
           inputlist = {}
           ]
]; (* end if useoldalgorithm case 2 *)
(* case with parameters, and sol1 had no or in it, *)
(* end new routine for compatibility analysis *)
        Clear[nrules1,nrules2,solc,newmaineqlist];
        ], (* end if 5 case where sol1 had or in it *) 
           (* else 4 sol1 was false *)
     Print["The system becomes inconsistent, or the compatibility conditions"];
     Print["require that one or more of the parameters are zero."];
     Print["Not acceptable!"];
     ] (* end if 4 *) 
     ], (* closes While *)  (* else 8 start case without parameters *)
      maineqlist = originalmaineqlist;
      If[debugsolve, 
        Print["At pt. OLDALGO or NEWALGO SO39,"<>
              " before solving (case without parameters), maineqlist: "];
        Print[maineqlist]
        ];
      If[debugsolve, 
        Print["At pt. OLDALGO or NEWALGO SO40,"<>
              " before solving (case without parameters), unknownlist: "];
        Print[unknownlist]
        ];
   nsolc = Flatten[Solve[maineqlist,unknownlist]];
      If[debugsolve, 
        Print["At pt. OLDALGO or NEWALGO SO41, after solving"<>
              " (case without parameters), nsolc: "];
        Print[nsolc]
        ];
   Print["Solution of the system: "];
   Print[nsolc];
   evaluate[formrho[x,t],rhot[x,t],nsolc,rules];
   Clear[solc1,solc2,solc3,nsolc]
   ] (* end 8 end if parameter parameter test *) 
  ] (* end point 1 *) 
]; (* end module main *)
 
(* ------------------------------------------------------------------------- *)

(* set-up for collecting the data in a log file, transcript of computations *)
logfile="";
OpenLog[filename_String] := (logfile = OpenWrite[filename];
  If[logfile === $Failed, Return[]];
  AppendTo[$Echo,logfile]; AppendTo[$Output,logfile];);
CloseLog[] := ($Echo = Complement[$Echo,{logfile}];
  $Output = Complement[$Output,{logfile}];
  Close[logfile];);

(* End of the auxiliary functions and procedures *)
(* Start of the executable code lines *)

menu;
OpenLog[myfile];
commentinter[]; 

(* HERE INITIALLY *)
(* 07/14/2006 WH setting up empty lists for listallrhos, listallfluxes *)
listallrhos = {};
listallfluxes = {};

If[formrho[x,t] === 0, formrho[x,t] = {}];
If[Not[ListQ[formrho[x,t]]], formrho[x,t] = {formrho[x,t]}];
Print["Working with the data file for the ",name,"."]; 
For[i = 1,i <= noeqs,i++,
   eq[i][x,t] = Expand[eq[i][x,t]];
   Print["Equation ",i," of the system with ",noeqs," equation(s):"];
   Print["  ",pdeform[eq[i][x,t]]," = 0"];
   ];
eq[x,t] = Sum[eq[i][x,t],{i,1,noeqs}];
highestordereq = highestorder[eq[x,t],False,False,False,False,True,False];
Print["Highest order of derivative in the system is ",highestordereq,"."];
If[formrho[x,t] === {}, (* start point rhogiven *)
  eqlist = Flatten[Table[eq[i][x,t],{i,1,noeqs}]];
  scalerules = scaling[eqlist];
  counternegweight = 0;
  listfreeweights = {};
  varscalelist = {};
  allchoices = {};
  Print["Program determines the weights of the variables (and parameters)."];
  For[i = 1,i <= noeqs,i++,
     weightu[i] = weightu[i] /. scalerules;
     If[i === 1, Print["For the given system:"] ];
     Print["* weight of u[",i,"] is ",weightu[i],"."]; 
     varscalelist = Union[varscalelist,{{u[i][x,t],weightu[i]}}];
     If[Not[NumberQ[weightu[i]]],
       If[Length[weightu[i]] === 1,
         listfreeweights = Union[listfreeweights,{weightu[i]}]],
       If[weightu[i] < 0,counternegweight++];
       ];
     ];
  lenweightpars = Length[weightpars];
  For[i = 1,i <= lenweightpars,i++,
     weight[Part[weightpars,i]] = 
             weight[Part[weightpars,i]] /. scalerules; 
     Print["* weight of ",Part[weightpars,i]," is ",
           weight[Part[weightpars,i]],"."];
     varscalelist = Union[varscalelist,
                    {{Part[weightpars,i],weight[Part[weightpars,i]]}}];
     If[Not[NumberQ[weight[Part[weightpars,i]]]],
       If[Length[weight[Part[weightpars,i]]] === 1,
         listfreeweights = 
            Union[listfreeweights,{weight[Part[weightpars,i]]}]],
       If[weight[Part[weightpars,i]] < 0,counternegweight++];
       ];
     ]; 
  weight[d/dt] = weight[d/dt] /. scalerules;
  Print["* weight of d/dt is ",weight[d/dt],"."];
  tvarscalelist = Union[varscalelist,{{"d/dt",weight[d/dt]}}];
  If[Not[NumberQ[weight[d/dt]]],
    If[Length[weight[d/dt]] === 1,
      listfreeweights = Union[listfreeweights,{weight[d/dt]}]],
    If[weight[d/dt] < 0,
      Print["CAUTION! Weight of d/dt is negative. Computations"];
      Print["will continue, but the program is not designed for this"]'
      Print["case. Results may not be completely reliable."];
      ];
    ];
  counterfreeweight = Length[listfreeweights];
  If[counternegweight =!= 0, (* start of if 1 *)
    Print["One or more of the weights are negative."];
    Print["Negative weights are not allowed."];
    Print["Aborting the computations!"];
    CloseLog[]; Abort[], 
    If[counterfreeweight =!= 0, (* start if 2 *)
      If[counterfreeweight > 1, (* start if 3 *)
        Print["Two or more of weights have freedom."];
        Print["Enter your values for the free weights or type Abort[];"];
        Print["Your values for the free weights can be entered by typing"];
        Print["`weightu[label] = val' or `weight[variable] = val' and"];
        Print["putting `;' between your entries."];
        Input[": "], (* else of counterfreeweight > 1 *)
        Print["Program will try to determine CHOICES for ",
              listfreeweights[[1]],":"];
        tempvarscalelist = Sort[Table[Part[tvarscalelist,k][[2]], 
                                 {k,1,Length[tvarscalelist]}]];
        tempvarscalelist = 
               Complement[tempvarscalelist,Select[tempvarscalelist,NumberQ]];
        lentempvarscalelist = Length[tempvarscalelist];
        Do[ (* begin do loop *)
          Print["* CHOICE ",k];
          Print["Solving the equation: ",Part[tempvarscalelist,k]," = 1."];
          attempt = Flatten[Solve[Part[tempvarscalelist,k] == 1]];
          If[attempt =!= {},
            solfreeweight = Part[listfreeweights,1] /. attempt;
            While[solfreeweight <= 0, 
                 solfreeweight = solfreeweight + 1;
                 Print["Since the weight was zero or negative,"];
                 Print["it was incremented with 1."]
                 ];
            Print[Part[listfreeweights,1]," = ",solfreeweight];
            weightrule = Part[listfreeweights,1] -> solfreeweight;
            tempvarscalelistval = tempvarscalelist /. weightrule;
            If[MemberQ[NonNegative[tempvarscalelistval],False],
              Print["Choice is rejected!"],
              varscalelistval = varscalelist /. weightrule;
              tvarscalelistval = tvarscalelist /. weightrule;    
              Print["List of all the variables (and parameters) with their"]; 
              Print["weights:"];
              Print[tvarscalelistval];
              testvarscalelist = Table[Part[varscalelistval,k][[2]],
                        {k,1,Length[varscalelist]}];
              If[Union[Positive[testvarscalelist]] === {False},
                Print["Since all the weights of u[i] are zero"];
                Print["this choice is rejected!"], (* else *)
                allchoices = Union[allchoices,{solfreeweight}];
                ];
              ];
            ],{k,1,lentempvarscalelist} 
        ];  (* end do loop *)
        Print["Simple POSITIVE choices are considered."];
        intchoices = Select[allchoices,IntegerQ];
        fracchoices = Complement[allchoices,intchoices];
        If[intchoices =!= {},
          (* Picking the minimum integer choice *)
          intchoices = Min[intchoices]; 
          choicerule = Part[listfreeweights,1] -> intchoices, (* else *)
          If[fracchoices =!= {},
           (* Picking the minimum fractional choice *)
           fracchoices = Min[fracchoices]; 
           choicerule = Part[listfreeweights,1] -> fracchoices, (* else *)
           Print["Enter your choice for the value of ",Part[listfreeweights,1],
                 "by typing its value. NO semi-colon at the end!"];
           choicerule = Part[listfreeweights,1]-> Input[": "];
           ]; 
          ]; 
        varscalelistchoice = varscalelist /. choicerule;
        tvarscalelistchoice = tvarscalelist /. choicerule;
        If[Length[allchoices] > 1,
          Print["List of all positive choices considered: "];
          Print[Union[Table[Part[listfreeweights,1] -> Part[allchoices,k], 
                             {k,1,Length[allchoices]}]]];
          Print["Some of the weights of u[i], but not all, could be zero."];
          Print["In the data file you could enter your choices in the format"];
          Print["`weightu[label] = value;' or `weight[variable] = value;'."];
          Print["Program continues with the choice: ", choicerule];
          Print["corresponding to the weights:"];
          Print[tvarscalelistchoice]
          ];
        varscalelist = varscalelistchoice
      ]  (* end if 3 *)
   ] (* end if 2 *)
  ]; (* end if 1 *)
  Print["The rank of rho should be an integer multiple of the lowest weight"];
  Print["of the DEPENDENT variable(s). Fractional weights are allowed."];
  If[
    formrho[x,t] === {},
    If[Not[NumberQ[rhorank]], rhorank = Input["Enter the rank of rho: "]]
    ];
  If[formrho[x,t] === {},
    Print["Computation of the density (and flux) of RANK = ",rhorank]; 
    ];  
  nodims = noeqs + Length[weightpars];
  varscalelist = Reverse[Sort[varscalelist,OrderedQ[{#1[[2]],#2[[2]]}]&]];
  constructformrho[nodims,varscalelist], (* else point for rhogiven *)
  formrho[x,t] = Part[formrho[x,t],1];
  formrho[x,t] = Expand[formrho[x,t]];
  If[Head[formrho[x,t]] === Times, 
    lenformrho = 1,lenformrho = Length[formrho[x,t]];
    ];
  Print["The program will only test the given form of the density rho."];
  Print["No determination of scaling properties."];
  Print["You have given this form for rho :"];
  Print[pdeform[formrho[x,t]]]
  ]; (* end point rhogiven *)

highestorderrho=highestorder[formrho[x,t],False,False,False,False,True,False];
Print["The highest-order of the terms in rho is ",
highestorderrho,"."];
eqlist = Flatten[Table[eq[i][x,t],{i,1,noeqs}]];
unknownlist = Table[c[i],{i,1,lenformrho}];
parameters = Union[parameters,weightpars];
weightpars = {};
 
Clear[i,k,lenformrho,counternegweight,counterfreeweight,scalerules,
     lenweightpars,varscalelist,tempvarscalelist,tvarscalelist,
     attempt,listfreeweights,allchoices,varscalelistval,temparscalelistval,
     tvarscalelistval,intchoices,fracchoices,choicerule,weightrule,
     nodims,highestordereq,weightu,lentempvarscalelist,solfreeweight,
     varscalelistchoice,tvarscalelistchoice,testvarscalelist,name,weight];

(* Solving for u[i]_t and calculation of the x-derivatives *)
mainsolforlist = Table[D[u[i][x,t],t],{i,1,noeqs}]; 
If[debug,
  Print["*************************************************"];
  Print["Trying to solve the equation(s) for ",mainsolforlist,","];
  Print["respectively."]; 
  Print["*************************************************"] 
  ];
utlist = {};
For[i = 1,i <= noeqs,i++,
   ut[i] = Solve[eq[i][x,t] == 0,mainsolforlist[[i]]][[1]];
   utlist = Flatten[Append[utlist,ut[i]]];
   ];
Clear[ut,mainsolforlist,i];
main[utlist];

(* 07/14/2006 WH added several statements to the summary *)
Print["********************* BEGIN SUMMARY **********************"];
Print["Total CPU time used in the current session is ",Round[TimeUsed[]],
      " seconds."];
If[Length[listallrhos] === 1, (* then *)
  Print["The software computed 1 density with its corresponding flux."], 
  (* else *)
  Print["The software computed ",Length[listallrhos]," different densities"<>
        " with their corresponding fluxes."]
  ];
Print["To see a list of all the densities type: listallrhos," ];
Print["or subscriptform[listallrhos]."];
Print["To see the last density type: rho[x,t], or subscriptform[rho[x,t]]."];
Print["To split the last density in independent pieces type: "];
Print["splitrho[rho[x,t]], or splitrho[subscriptform[rho[x,t]]]."];
Print["To see a list of all the fluxes type: listallfluxes,"];
Print["or subscriptform[listallfluxes]."];
Print["To see the last flux type: j[x,t], or subscriptform[j[x,t]]."];
Print["To split the last flux in independent pieces type: splitflux[j[x,t]],"];
Print["or splitflux[subscriptform[j[x,t]]]."];
Print["************************ END SUMMARY ***************************"];
CloseLog[];
(* WH 08/04/2007 print statement added *)
Print["Package condens.m of August 4, 2007 successfully loaded."];
(* ******************************* END ************************************* *)
