/*-*-MACSYMA-*-*/
EVAL_WHEN(BATCH,TTYOFF:TRUE)$
/* Or use the BATCHLOAD command to load this with TTYOFF:TRUE */
/*  NOTE:  THE CURRENT VERSION OF VECT IS THE ONE DUE TO STOUTEMYER.
IT WILL BE REPLACED SOON BY AN EXTENDED VERSION WHICH HANDLES BOTH
VECTORS AND DYADICS.
	MICHAEL C. WIRTH (MCW)
	12/18/78
Style changes made in order to TRANSLATE, 3/1/81 George Carrette (GJC)
*/


HERALD_PACKAGE('VECT)$
PUT('VECT,TRUE,'VERSION);

EVAL_WHEN([TRANSLATE,batch,demo,load,loadfile],
MATCHDECLARE([ETRUE,TTRUE,VTRUE], TRUE,
LESSP, BEFORE, SCALARM, VSCALARP),
	TR_BOUND_FUNCTION_APPLYP:FALSE
	/* we do not want FOO(F):=F(1) to mean APPLY(F,[1]) */
	)$

INFIX("~", 134, 133, EXPR, EXPR, EXPR) $
PREFIX("GRAD", 142, EXPR, EXPR) $
PREFIX("DIV", 142, EXPR, EXPR) $
PREFIX("CURL", 142, EXPR, EXPR) $
PREFIX("LAPLACIAN", 142, EXPR, EXPR) $ 

DECLARE([".",ORDER],COMMUTATIVE,
        ORDERN,NARY,
        ["GRAD","DIV","CURL","LAPLACIAN"],OUTATIVE,
        "CURL",NONSCALAR)$

TELLSIMPAFTER(0~ETRUE, 0) $
TELLSIMPAFTER(ETRUE~0, 0) $
TELLSIMPAFTER(ETRUE~ETRUE, 0)$
TELLSIMPAFTER(ETRUE~TTRUE.VTRUE, ETRUE.TTRUE~VTRUE)$
TELLSIMP(ETRUE~LESSP, -LESSP~ETRUE) $
TELLSIMPAFTER(DIV (CURL( ETRUE)), 0) $
TELLSIMPAFTER(CURL (GRAD( ETRUE)), 0) $
TELLSIMPAFTER(VECTORPOTENTIAL(ETRUE,TTRUE),
   (SCALEFACTORS(TTRUE), VPOT(ETRUE))) $
TELLSIMPAFTER(VECTORPOTENTIAL(ETRUE), VPOT(ETRUE)) $
TELLSIMPAFTER(EXPRESS(ETRUE), EXPRESS1(ETRUE)) $
TELLSIMPAFTER(EXPRESS(ETRUE,TTRUE),
   (SCALEFACTORS(TTRUE), EXPRESS(ETRUE))) $
TELLSIMPAFTER(POTENTIAL(ETRUE, TTRUE),
   (SCALEFACTORS(TTRUE),  POTENTIAL1(ETRUE))) $
TELLSIMPAFTER(POTENTIAL(ETRUE), POTENTIAL1(ETRUE)) $


/* Variables and switches */

DEFINE_VARIABLE(COORDINATES, '[X, Y, Z],ANY)$
DEFINE_VARIABLE(DIMENSION,3,FIXNUM)$
DEFINE_VARIABLE(DIMENIMBED,1,FIXNUM)$
DEFINE_VARIABLE(TRYLENGTH,1,FIXNUM)$
DEFINE_VARIABLE(BESTLENGTH,1,FIXNUM)$
DEFINE_VARIABLE(sfprod,1,any)$
DEFINE_VARIABLE(sf,[1,1,1],list)$

DEFINE_VARIABLE(EXPANDALL,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDDOT,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDDOTPLUS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDGRAD,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDPLUS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDALL,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDGRADPLUS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDDIV,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDDIVPLUS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDCURL,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDCURLPLUS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDLAPLACIAN,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDLAPLACIANPLUS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDPROD,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDGRADPROD,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDDIVPROD,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDCURLCURL,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDLAPLACIANTODIVGRAD,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDLAPLACIANPROD,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDCROSS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDCROSSCROSS,false,BOOLEAN)$
DEFINE_VARIABLE(EXPANDCROSSPLUS,false,BOOLEAN)$
DEFINE_VARIABLE(FIRSTCROSSSCALAR,false,BOOLEAN)$


EV_DIFF(X):=APPLY('EV,[X,'DIFF])$


declare(jacobian,special)$
SCALEFACTORS(TRANSFORMATION) := BLOCK(
   IF LISTP(FIRST(TRANSFORMATION)) THEN (
      COORDINATES: REST(TRANSFORMATION),
      TRANSFORMATION: FIRST(TRANSFORMATION))
   ELSE COORDINATES: LISTOFVARS(TRANSFORMATION),
   DIMENSION: LENGTH(COORDINATES),
   DIMENIMBED: LENGTH(TRANSFORMATION),
   FOR ROW:1 THRU DIMENSION DO
      FOR COL:1 THRU DIMENIMBED DO JACOBIAN[ROW,COL]:
         TRIGSIMP(RATSIMP(DIFF(TRANSFORMATION[COL],
            COORDINATES[ROW]))),
   SFPROD:1,
   FOR ROW:1 THRU DIMENSION DO (
      FOR COL:1 THRU ROW-1 DO (
         SF[ROW]: GCOV(ROW,COL),
         IF SF[ROW]#0 THEN PRINT("WARNING: COORDINATE SYSTEM IS NONORTHOGONAL UNLESS FOLLOWING SIMPLIFIES TO ZERO:", SF[ROW])),
      SF[ROW]: RADCAN(SQRT(GCOV(ROW,ROW))),
      SFPROD: SFPROD*SF[ROW])) $

GCOV(II,JJ) := TRIGSIMP(RATSIMP(SUM(
   JACOBIAN[II,KK]*JACOBIAN[JJ,KK], KK, 1, DIMENIMBED))) $

EXPRESS1(EXPN) := BLOCK(
   [ANS],
   IF MAPATOM(EXPN) THEN
      IF NONSCALARP(EXPN) THEN (ANS:[],
         FOR JJ: DIMENSION STEP -1 THRU 1 DO
            ANS: CONS(EXPN[COORDINATES[JJ]], ANS),
         RETURN(ANS))
      ELSE RETURN(EXPN),
   EXPN: MAP('EXPRESS1, EXPN),
   IF MAPATOM(EXPN) OR LISTP(EXPN) THEN RETURN(EXPN),

   IF INPART(EXPN,0)="GRAD" THEN (ANS:[],
      EXPN: INPART(EXPN,1),
      FOR JJ: DIMENSION STEP -1 THRU 1 DO ANS:
         CONS('DIFF(EXPN,COORDINATES[JJ])/SF[JJ], ANS),
      RETURN(ANS)),

   IF PIECE="DIV" THEN (EXPN: INPART(EXPN,1),
      IF NOT LISTP(EXPN) THEN ERROR("DIV CALLED ON SCALAR ARG:",
         EXPN),
      RETURN(SUM('DIFF(SFPROD*EXPN[JJ]/SF[JJ],
         COORDINATES[JJ]), JJ, 1, DIMENSION)/SFPROD)),

   IF PIECE="LAPLACIAN" THEN RETURN(SUM('DIFF(SFPROD*'DIFF(
      INPART(EXPN,1),COORDINATES[JJ])/SF[JJ]**2,
      COORDINATES[JJ]), JJ, 1, DIMENSION) / SFPROD),

   IF PIECE="CURL" THEN (EXPN:INPART(EXPN,1),
      IF LISTP(EXPN) THEN (
         IF LENGTH(EXPN)=2 THEN RETURN(('DIFF(SF[2]*EXPN[2],
            COORDINATES[1])-'DIFF(SF[1]*EXPN[1],
            COORDINATES[2]))/ SF[1]/SF[2]),
         IF DIMENSION=3 THEN RETURN([
             ('DIFF(SF[3]*EXPN[3],COORDINATES[2])-
             'DIFF(SF[2]*EXPN[2],COORDINATES[3]))/
             SF[2]/SF[3],
             ('DIFF(SF[1]*EXPN[1],COORDINATES[3])-
              'DIFF(SF[3]*EXPN[3],COORDINATES[1]))/
             SF[1]/SF[3],
             ('DIFF(SF[2]*EXPN[2],COORDINATES[1]) -
              'DIFF(SF[1]*EXPN[1],COORDINATES[2]))/
             SF[1]/SF[2]])),
      ERROR("CURL USED IN SPACE OF WRONG DIMENSION")),

   IF PIECE="~" THEN (
      ANS: INPART(EXPN,1),  EXPN:INPART(EXPN,2),
      IF LISTP(ANS) AND LISTP(EXPN) AND LENGTH(ANS)=LENGTH(EXPN)
         THEN (IF LENGTH(ANS)=2 THEN RETURN(ANS[1]*EXPN[2]
             -ANS[2]*EXPN[1]),
            IF LENGTH(ANS)=3 THEN RETURN([ANS[2]*EXPN[3]-
               ANS[3]*EXPN[2], ANS[3]*EXPN[1]-ANS[1]*EXPN[3],
               ANS[1]*EXPN[2]-ANS[2]*EXPN[1]])),
      ERROR("~ USED WITH IMPROPER ARGUMENTS:",ANS,EXPN)),

   EXPN) $

TRIGSIMP(EXPN) :=
   RATSIMP(TRIGSIMP1(NUM(EXPN))/TRIGSIMP1(DENOM(EXPN))) $

TRIGSIMP1(EXPN) := BLOCK(
   [LISTOFTRIGSQ, BESTLENGTH, TRYLENGTH],
   LISTOFTRIGSQ: LISTOFTRIGSQ(EXPN),
   BESTLENGTH: 999999,
   IF LISTOFTRIGSQ#[] THEN IMPROVE(EXPN, LISTOFTRIGSQ),
   RETURN(EXPN)) $

IMPROVE(SUBSOFAR, LISTOFTRIGSQ) :=
   IF LISTOFTRIGSQ=[] THEN (
      TRYLENGTH:mode_identity(fixnum, EXPNLENGTH(SUBSOFAR)),
      IF TRYLENGTH<BESTLENGTH THEN (
         EXPN: SUBSOFAR,
         BESTLENGTH: TRYLENGTH))
   ELSE (IMPROVE(SUBSOFAR, REST(LISTOFTRIGSQ)),
      FOR ALT IN FIRST(LISTOFTRIGSQ) DO 
         IMPROVE(RATSUBST(
            IF INPART(ALT,0)='SIN THEN 1-COS(INPART(ALT,1))**2
            ELSE IF PIECE='COS THEN 1-SIN(INPART(ALT,1))**2
            ELSE IF PIECE='SINH THEN COSH(INPART(ALT,1))**2-1
            ELSE 1+SINH(INPART(ALT,1))**2,
               ALT**2, SUBSOFAR), REST(LISTOFTRIGSQ))) $

LISTOFTRIGSQ(EXPN) :=
   IF ATOM(EXPN) THEN []
   ELSE BLOCK([INFLAG, ANS],
      IF INPART(EXPN,0)="^" AND INTEGERP(INPART(EXPN,2))
         AND PIECE>=2 THEN
         IF ATOM(EXPN:INPART(EXPN,1)) THEN RETURN([])
         ELSE IF MEMBER(INPART(EXPN,0),'[SIN,COS,SINH,COSH])
            THEN RETURN([[EXPN]]),
   INFLAG:TRUE,
      ANS:[],
      FOR ARG IN EXPN DO
         ANS: SPECIALUNION(LISTOFTRIGSQ(ARG), ANS),
      RETURN(ANS)) $

SPECIALUNION(LIST1,LIST2) :=
   IF LIST1=[] THEN LIST2
   ELSE IF LIST2=[] THEN LIST1
   ELSE BLOCK([ALTERNATES],
      ALTERNATES: FIRST(LIST1),
      FOR ALT IN ALTERNATES DO LIST2:
         IF INPART(ALT,0)='SIN THEN UPDATE(ALT,'COS,list2)
         ELSE IF PIECE='COS THEN UPDATE(ALT,'SIN,list2)
         ELSE IF PIECE='SINH THEN UPDATE(ALT,'COSH,list2)
         ELSE UPDATE(ALT,'SINH,list2),
      RETURN(SPECIALUNION(REST(LIST1),LIST2))) $

UPDATE(FORM, COMPLEMENT,list2) := BLOCK(
   [ANS],
   COMPLEMENT: APPLY('COMPLEMENT,[INPART(FORM,1)]),
   ANS: FOR ELEMENT IN LIST2 DO
      IF MEMBER(FORM, ELEMENT) THEN RETURN('FOUND)
      ELSE IF MEMBER(COMPLEMENT,ELEMENT) THEN RETURN(
         CONS([FORM,COMPLEMENT], DELETE(ELEMENT,LIST2))),
   IF ANS='FOUND THEN LIST2
   ELSE IF ANS='DONE THEN CONS([FORM],LIST2)
   ELSE ANS) $

/*EXPNLENGTH(EXPR):=LENGTH(?STRING(EXPR))$*/
EXPNLENGTH(EXPR) := 
   IF ATOM(EXPR) THEN 1
   ELSE 1 + ARGSLENGTH(SUBSTINPART("[", EXPR, 0)) $

ARGSLENGTH(ARGS) :=
   IF ARGS=[] THEN 0
   ELSE EXPNLENGTH(FIRST(ARGS)) + ARGSLENGTH(REST(ARGS)) $


DOTASSOC: DOTEXPTSIMP: FALSE$
DOTSCRULES: TRUE $
define_variable(EXPANDFLAGS, '[
   EXPANDALL,
      EXPANDDOT,
         EXPANDDOTPLUS,
      EXPANDCROSS,
         EXPANDCROSSPLUS,
         EXPANDCROSSCROSS,
      EXPANDGRAD,
         EXPANDGRADPLUS,
         EXPANDGRADPROD,
      EXPANDDIV,
         EXPANDDIVPLUS,
         EXPANDDIVPROD,
      EXPANDCURL,
         EXPANDCURLPLUS,
         EXPANDCURLCURL,
      EXPANDLAPLACIAN,
         EXPANDLAPLACIANPLUS,
         EXPANDLAPLACIANPROD,
   EXPANDLAPLACIANTODIVGRAD,
   EXPANDPLUS,
   EXPANDPROD ],list) $


APPLY('DECLARE, ['EXPANDFLAGS, 'EVFLAG]) $

VECTORSIMP(EXPN) := BLOCK(
   [DOTDISTRIB, DOTSCRULES, INFLAG, FIRSTCROSSSCALAR],
   INFLAG: FIRSTCROSSSCALAR: TRUE,
   DOTDISTRIB: EXPANDALL OR EXPANDDOT OR EXPANDDOTPLUS
      OR EXPANDPLUS,
   IF EXPANDALL OR EXPANDGRAD OR EXPANDGRADPLUS OR EXPANDPLUS
      THEN DECLARE("GRAD",ADDITIVE),
   IF EXPANDALL OR EXPANDDIV OR EXPANDDIVPLUS OR EXPANDPLUS
      THEN DECLARE("DIV",ADDITIVE),
   IF EXPANDALL OR EXPANDCURL OR EXPANDCURLPLUS OR EXPANDPLUS
      THEN DECLARE("CURL",ADDITIVE),
   IF EXPANDALL OR EXPANDLAPLACIAN OR EXPANDLAPLACIANPLUS OR
      EXPANDPLUS THEN DECLARE("LAPLACIAN",ADDITIVE),
   EXPN: VSIMP(EXPN),
   IF EXPANDALL THEN EXPN: RATEXPAND(EXPN),
   IF EXPANDALL OR EXPANDGRAD OR EXPANDGRADPLUS OR EXPANDPLUS
      THEN REMOVE("GRAD",ADDITIVE),
   IF EXPANDALL OR EXPANDDIV OR EXPANDDIVPLUS OR EXPANDPLUS
      THEN REMOVE("DIV",ADDITIVE),
   IF EXPANDALL OR EXPANDCURL OR EXPANDCURLPLUS OR EXPANDPLUS
      THEN REMOVE("CURL",ADDITIVE),
   IF EXPANDALL OR EXPANDLAPLACIAN OR EXPANDLAPLACIANPLUS OR
      EXPANDPLUS THEN REMOVE("LAPLACIAN",ADDITIVE),
   EXPN) $


BEFORE(ARG) := INPART(ORDER(ETRUE,ARG),1)#ETRUE$

VSCALARP(ARG) := NOT NONSCALARP(ARG)$

VSIMP(EXPN) :=
   IF MAPATOM(EXPN) THEN EXPN
   ELSE BLOCK([PV, QV, RV, SV],
      EXPN: MAP('VSIMP, EXPN),
      IF MAPATOM(EXPN) THEN RETURN(EXPN),
      IF INPART(EXPN,0)="~" THEN EXPN:REMOVECROSSSC1(EXPN,pv,rv,sv)
      ELSE IF PIECE="GRAD" THEN (
         IF (EXPANDALL OR EXPANDGRAD OR EXPANDGRADPROD OR
            EXPANDPROD) AND
            NOT MAPATOM(PV:INPART(EXPN,1)) AND INPART(PV,0)="*"
            THEN EXPN:APPLY("+", MAPLIST(lambda([u],GRADPROD(u,pv)), PV)))
      ELSE IF PIECE="DIV" THEN(
         IF (EXPANDALL OR EXPANDDIV OR EXPANDDIVPROD OR
            EXPANDPROD) AND NOT
            MAPATOM(PV:INPART(EXPN,1)) AND INPART(PV,0)="*" THEN
            EXPN: APPLY("+", MAPLIST(lambda([u],DIVPROD(u,pv)), PV)))
      ELSE IF PIECE="CURL" THEN (
         IF (EXPANDALL OR EXPANDCURL OR EXPANDCURLCURL) AND NOT
            MAPATOM(PV:INPART(EXPN,1)) AND INPART(PV,0)="CURL"
            THEN EXPN: GRAD( DIV(PV:INPART(PV,1))) - LAPLACIAN( PV))
      ELSE IF PIECE="LAPLACIAN" THEN
         IF EXPANDLAPLACIANTODIVGRAD THEN
            EXPN: DIV (GRAD( INPART(EXPN,1)))
         ELSE IF (EXPANDALL OR EXPANDLAPLACIAN OR
            EXPANDLAPLACIANPROD OR EXPANDPROD) AND NOT MAPATOM
            (PV:INPART(EXPN,1)) AND INPART(PV,0)="*" THEN(
            QV: INPART(PV,1),
            RV: DELETE(QV,PV),
            EXPN: RV*LAPLACIAN(QV) + 2*GRAD (RV) * GRAD( QV) + QV*
               LAPLACIAN(RV)),
      EXPN) $

CROSSSIMP(EX,pv,rv,sv) :=
   IF NOT MAPATOM(EX) AND INPART(EX,0)="~" THEN (
      IF EXPANDALL OR EXPANDCROSS OR EXPANDCROSSCROSS THEN
         EX: TRYCROSSCROSS(EX,pv,rv,sv),
      IF NOT MAPATOM(EX) AND INPART(EX,0)="~" AND(
         EXPANDALL OR EXPANDCROSS OR EXPANDCROSSPLUS OR
         EXPANDPLUS) THEN EX: TRYCROSSPLUS(EX,pv,rv,sv),
      EX)
   ELSE EX $

REMOVECROSSSC(EXPN,pv,rv,sv) :=
   IF NOT MAPATOM(EXPN) AND INPART(EXPN,0)="~" THEN
      REMOVECROSSSC1(EXPN,pv,rv,sv)
   ELSE EXPN $

REMOVECROSSSC1(EXPN,pv,rv,sv) :=  BLOCK(
      [LEFT, RIGHT],
      LEFT: PARTITIONSC(INPART(EXPN,1)),
      RIGHT: PARTITIONSC(INPART(EXPN,2)),
      IF FIRSTCROSSSCALAR AND (LEFT[2]=1 OR RIGHT[2]=1)
         THEN( PRINT("WARNING: DECLARE VECTOR INDETERMINANTS 
NONSCALAR TO AVOID ERRORS & TO GET FULL SIMPLIFICATION"),
            FIRSTCROSSSCALAR:FALSE,
                RETURN(EXPN)),
      LEFT[1]*RIGHT[1]*CROSSSIMP(LEFT[2]~RIGHT[2],pv,rv,sv))$

PARTITIONSC(EX) :=
   IF MAPATOM(EX) THEN
      IF NONSCALARP(EX) THEN [1,EX]
      ELSE [EX,1]
   ELSE IF INPART(EX,0)="*" THEN BLOCK([SC,NONSC],
      SC: NONSC: 1,
      FOR FACT IN EX DO
         IF NONSCALARP(FACT) THEN NONSC:NONSC*FACT
         ELSE SC:SC*FACT,
      [SC,NONSC])
   ELSE [1,EX] $

TRYCROSSPLUS(EXPN,pv,rv,sv) :=(
   PV:INPART(EXPN,1), RV:INPART(EXPN,2),
   IF NOT MAPATOM(PV) AND INPART(PV,0)="+" THEN
      IF NOT MAPATOM(RV) AND INPART(RV,0)="+" THEN
         MAP(lambda([u],TRYCROSSPLUS(u,pv,rv,sv)), 
         MAP(lambda([u],CROSSRV(u,rv,sv)), PV))
      ELSE MAP(lambda([u],CROSSRV(u,rv,sv)), PV)
   ELSE IF NOT MAPATOM(RV) AND INPART(RV,0)="+" THEN
      MAP(lambda([u],PVCROSS(pv,u,sv)), RV)
   ELSE EXPN) $

TRYCROSSCROSS(EXPN,pv,rv,sv) := (
      PV:INPART(EXPN,1),  RV:INPART(EXPN,2),
      IF NOT MAPATOM(RV) AND INPART(RV,0)="~" THEN (
         SV: INPART(RV,2), RV:INPART(RV,1),
         RV*PV.SV - SV*PV.RV)
      ELSE IF NOT MAPATOM(PV) AND INPART(PV,0)="~" THEN(
         SV:INPART(PV,2), PV:INPART(PV,1),
         SV*RV.PV - PV*RV.SV)
      ELSE EXPN) $

PVCROSS(pv,RV,sv) := REMOVECROSSSC(PV~RV,pv,rv,sv) $

CROSSRV(PV,rv,sv) := REMOVECROSSSC(PV~RV,pv,rv,sv) $

GRADPROD(UU,pv) := DELETE(UU,PV)*GRAD(UU) $

DIVPROD(UU,pv) := BLOCK([DOTSCRULES],
   DOTSCRULES: FALSE,
 
   IF NONSCALARP(UU) THEN DELETE(UU,PV)*DIV(UU)
   ELSE DELETE(UU,PV).GRAD(UU) )$


POTENTIAL1(GR) := BLOCK(
   [ORIGIN, GRPERM, JJ, RESULT,%dum],
   IF NOT LISTP(GR) OR LENGTH(GR)#DIMENSION THEN ERROR(
      "1ST ARG OF POTENTIAL MUST BE A LIST OF LENGTH EQUAL TO",
      "THE DIMENSIONALITY OF THE COORDINATE SYSTEM"),
   ORIGIN: ZEROLOC(),
   RESULT: [],
   FOR JJ:DIMENSION STEP -1 THRU 1 DO
      RESULT: CONS(SF[JJ]*GR[JJ], RESULT),
   GRPERM:[],
   FOR EQN IN ORIGIN DO (
      JJ:1,
      WHILE LHS(EQN)#COORDINATES[JJ] DO JJ:JJ+1,
      GRPERM: ENDCONS(RESULT[JJ], GRPERM)),
   RESULT:SUM(MYINT(SUBLESS(JJ,origin,grperm), %DUM, RHS(ORIGIN[JJ]),
      LHS(ORIGIN[JJ])), JJ, 1, DIMENSION),
   GR: GR-EXPRESS1(GRAD( RESULT)),
   GR: EV_DIFF(GR),
   GR: TRIGSIMP(RADCAN(GR)),
   ORIGIN:1,
    /* Variable name should NOT be re-used! */
   WHILE ORIGIN<=DIMENSION AND GR[ORIGIN]=0 DO 
         (MODE_DECLARE(ORIGIN,FIXNUM),ORIGIN:ORIGIN+1),
   IF ORIGIN<=DIMENSION THEN PRINT("UNABLE TO PROVE THAT THE",
      "FOLLOWING DIFFERENCE BETWEEN THE INPUT AND THE GRADIENT",
      "OF THE RETURNED RESULT IS ZERO", GR),
   TRIGSIMP(RADCAN(RESULT))) $

define_variable(POTENTIALZEROLOC, 0, any) $

ZEROLOC() := 
   IF NOT LISTP(POTENTIALZEROLOC) THEN
       MAP(LAMBDA([UU],UU=POTENTIALZEROLOC), COORDINATES)
   ELSE IF DISJUNCT(COORDINATES,MAP('LHS,POTENTIALZEROLOC)) # []
      THEN ERROR("POTENTIALZEROLOC MUST BE A LIST OF LENGTH",
      "EQUALING THE DIMENSIONALITY OF THE COORDINATE SYSTEM",
      "CONTAINING EQUATIONS WITH EACH COORDINATE VARIABLE",
      "ON THE LHS OF EXACTLY 1 EQUATION,",
      "OR ELSE POTENTIALZEROLOC MUST NOT BE A LIST")
   ELSE POTENTIALZEROLOC$



EVAL_WHEN([TRANSLATE,BATCH,DEMO,load,loadfile],

CYC(II) ::= BUILDQ([II], 1 + REMAINDER(II+SHIFT,3)) )$

VPOT(KURL) := BLOCK(
   [ORIGIN, SHIFT],
   MODE_DECLARE(SHIFT,FIXNUM),
   IF NOT LISTP(KURL) OR LENGTH(KURL)#3 THEN ERROR(
      "1ST ARG OF VECTORPOTENTIAL MUST BE A LIST OF LENGTH 3"),
   ORIGIN: ZEROLOC(),
   SHIFT: 1,
   WHILE SHIFT<=3 AND LHS(ORIGIN[1])#COORDINATES[SHIFT] DO 
      SHIFT:SHIFT+1,
   SHIFT: SHIFT+1,
   IF SHIFT>4 OR LHS(ORIGIN[2])#COORDINATES[CYC(2)] OR
      LHS(ORIGIN[3])#COORDINATES[CYC(3)] THEN ERROR(
      "LEFT SIDES OF POTENTIALZEROLOC MUST BE A CYCLIC",
      "PERMUTATION OF COORDINATES"),
   ORIGIN: [(MYINT(SF[CYC(1)]*SF[CYC(3)]*KURL[CYC(2)],
      LHS(ORIGIN[3]),RHS(ORIGIN[3]),LHS(ORIGIN[3])) - MYINT(
      SF[CYC(1)]*SF[CYC(2)]*SUBST(ORIGIN[3],KURL[CYC(3)]),
      LHS(ORIGIN[2]),RHS(ORIGIN[2]),LHS(ORIGIN[2])))/SF[CYC(1)],
      -MYINT(SF[CYC(2)]*SF[CYC(3)]*KURL[CYC(1)],
      LHS(ORIGIN[3]),RHS(ORIGIN[3]),LHS(ORIGIN[3]))/SF[CYC(2)],
      0],
   ORIGIN: [ORIGIN[CYC(CYC(1))], ORIGIN[CYC(CYC(2))],
      ORIGIN[CYC(CYC(3))]],
   KURL: KURL-EXPRESS1(CURL (ORIGIN)),
   KURL: EV_DIFF(KURL),
   KURL: TRIGSIMP(RADCAN(KURL)),
   FOR JJ:1 THRU 3 DO IF KURL[JJ]#0 THEN PRINT(
      "UNABLE TO PROVE THAT THE FOLLOWING DIFFERENCE BETWEEN A",
      "COMPONENT OF THE INPUT AND OF THE CURL OUTPUT IS ZERO",
      KURL[JJ]),
   ORIGIN) $



DISJUNCT(L1,L2) := APPEND(SETDIFF(L1,L2), SETDIFF(L2,L1)) $

SETDIFF(L1,L2) :=
   IF L1=[] THEN []
   ELSE IF MEMBER(FIRST(L1),L2) THEN SETDIFF(REST(L1),L2)
   ELSE CONS(FIRST(L1), SETDIFF(REST(L1),L2)) $

SUBLESS(KK,origin,grperm) := (MODE_DECLARE(KK,FIXNUM),BLOCK([ANS,%dum],
   ANS: RATSUBST(%DUM, LHS(ORIGIN[KK]), GRPERM[KK]),
   FOR L1:1 THRU KK-1 DO
      ANS: RATSUBST(RHS(ORIGIN[L1]), LHS(ORIGIN[L1]),ANS),
   ANS)) $

MYINT(FUN,VAR,LOW,HIGH):=BLOCK([RESULT,ATLOW,ATHIGH],
  RESULT:INTEGRATE(FUN,VAR),
  IF FREEOF(nounify('INTEGRATE),RESULT) THEN (
	ATLOW:EVLIMIT(RESULT,VAR,LOW),
	IF ATLOW=FALSE THEN GO(NOGOOD),
	ATHIGH:EVLIMIT(RESULT,VAR,HIGH),
	IF ATHIGH=FALSE THEN GO(NOGOOD),
	RETURN(RADCAN(ATHIGH-ATLOW))),
  NOGOOD, DEFINT(FUN,VAR,LOW,HIGH))$

EVLIMIT(EXPR,VAR,LIM):=BLOCK([TEMP],
  IF LIM='MINF OR LIM='INF THEN GO(USELIMIT),
  TEMP:ERRCATCH(SUBST(LIM,VAR,EXPR)),
  IF TEMP#[] THEN RETURN(TEMP[1]),
  USELIMIT, TEMP:LIMIT(EXPR,VAR,LIM),
  IF MEMBER(TEMP,'[INF,MINF,UND,IND,INFINITY]) THEN RETURN(FALSE),
  IF FREEOF(nounify('LIMIT),TEMP) THEN TEMP)$

EVAL_WHEN(BATCH,TTYOFF:FALSE)$