(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialiation Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)







If[$VersionNumber<3,
  Print["This package will not work correctly with versions of Mathematica \
earlier than Mathematica 3.0."]]

Print["HFT; Version 5.01, 12 June 1999"]

Print["This Mathematica package, called HFT.m, is designed for computing with \
harmonic functions. Documentation for the use of this package and information \
about the algorithms used in it is available in the Mathematica notebook \
called HFT.nb."]

Print["For addional information about harmonic functions, see the book \
Harmonic Function Theory, by Sheldon Axler, Paul Bourdon, and Wade Ramey, \
published by Springer."]

Print["This package is copyrighted by Sheldon Axler but is distributed \
electonically without charge. The most recent version of this HFT.m package \
along with its documentation notebook HFT.nb is available over the world wide \
web at: http://math.sfsu.edu/axler"]

Print["Comments, suggestions, and bug reports should be sent by electronic \
mail to: axler@sfsu.edu"]

Print["* The computer is unpacking the HFT.m package. This may take a minute \
or so."]

BeginPackage["HFT`"]



AnnularDirichlet::"usage"=
  "AnnularDirichlet[ f, g, r, R, x ] gives the function of x "<>
    "harmonic on an annular region that equals f on the inner "<>
    "sphere of radius r and equals g on the outer sphere of radius R."

AntiLaplacian::"usage"=
  "AntiLaplacian[f, x] gives a polynomial "<>"whose Laplacian equals f.  "<>
    "Here f must be a polynomial function of x."

Ball::"usage"=
  "Ball is one of the values that may be assigned to the option "<>
    "Orthonormal."

BasisH::"usage"=
  "BasisH[m, x] gives a basis for the space of harmonic "<>
    "polynomials homogeneous of degree m in the variable x."

BergmanKernel::"usage"=
  "BergmanKernel[x, y] gives the Bergman "<>
    "reproducing kernel for the unit ball."

BergmanKernelH::"usage"=
  "BergmanKernel[z, w] gives the Bergman "<>
    "reproducing kernel for the upper half-space."

BergmanProjection::"usage"=
  "BergmanProjection[ f, x] gives the "<>
    "orthogonal projection of a polynomial f, as a function of x, "<>
    "onto the Bergman space of harmonic functions on the ball."

BiDirichlet::"usage"=
  "BiDirichlet[f , x] solves the "<>
    "BiDirichlet problem with boundary data f, "<>"as a function of x."

Delta::"usage"=
  "Delta[j] is the vector that equals 1 in the j-th "<>
    "coordinate and 0 in the other coordinates."

Dimension::"usage"="Dimension[x] is the dimension of the vector x."

DimensionH::"usage"=
  "DimensionH[m, n] gives the vector space dimension of the "<>
    "space of spherical harmonics of degree m in n variables."

Dirichlet::"usage"=
  "Dirichlet[f, x] gives the harmonic function "<>
    "that equals f on the unit sphere.  Here f must be a "<>
    "polynomial function of x."

Divergence::"usage"=
  "Divergence[f, x] gives the divergence of f "<>"with respect to x."

DoubleBracketingBar::"usage"="DoubleBracketingBar is equal to the Norm."

ExpandNorm::"usage"=
  "ExpandNorm[f] gives f with all terms of the form |x + y| "<>
    "replaced by Sqrt[ |x|^2 + |y|^2 + 2 x.y ]."

ExteriorDirichlet::"usage"=
  "ExteriorDirichlet[f , x] solves the "<>
    "exterior Dirichlet problem with boundary data f, "<>"as a function of x."

ExteriorNeumann::"usage"=
  "ExteriorNeumann[f , x] solves the "<>
    "exterior Neumann problem with boundary data f, "<>"as a function of x."

Grad::"usage"="Grad[f, x] gives the gradient of f with "<>"respect to x."

HarmonicConjugate::"usage"=
  "HarmonicConjugate[u, {x,y}] gives "<>
    "the harmonic conjugate of u with respect to the variables x, y."

HarmonicDecomposition::"usage"=
  "HarmonicDecomposition[u, x] "<>
    "gives the harmonic decomposition of u with respect to "<>
    "the variable x."

HilbertSchmidt::"usage"=
  "HilbertSchmidt[A] gives the "<>"HilbertSchmidt norm of a matrix A."

Homogeneous::"usage"=
  "Homogeneous[f, d, x] gives the homogenous "<>
    "part of f of degree d, with respect to x."

IntegrateBall::"usage"=
  "IntegrateBall[f, x] gives the integral "<>
    "of f, as a function of x, over the unit ball."

IntegrateSphere::"usage"=
  "IntegrateSphere[f, x] gives the "<>
    "integral of f, as a function of x, over the unit sphere "<>
    "with respect to normalized surface area measure."

Inversion::"usage"="Inversion[x] gives the inversion of x."

J::"usage"="J[f, x] gives the Jacobian of f with "<>"respect to x."

Kelvin::"usage"=
  "Kelvin[u, x] gives the Kelvin transform of u, "<>
    "thought of as a function of x."

KelvinM::"usage"=
  "KelvinM is the modified Kelvin transform, "<>
    "as defined in Chapter 7 of Harmonic Function Theory."

Laplacian::"usage"=
  "Laplacian[f, x] gives the Laplacian of f "<>"with respect to x."

Multiple::"usage"=
  "Multiple is an option for AntiLaplacian.  "<>
    "The default value is None.  The value Norm^2 produces the "<>
    "unique antiLaplacian that is a polynomial multiple of "<>
    "Norm[x]^2, where x is the variable."

Neumann::"usage"=
  "Neumann[f, g, x] solves the "<>
    "Neumann problem of finding a function of x whose "<>
    "outward normal derivative (on the unit sphere) is f "<>
    "and whose Laplacian is g."

Norm::"usage"="Norm[x] gives the Euclidean norm of x."

NormalD::"usage"=
  "NormalD[f, z] gives the outward normal derivative of f, "<>
    "as a function of z, with respect to the unit ball."

Orthonormal::"usage"=
  "Orthonormal is an option for BasisH.  "<>
    "The default value is None.  The value Ball produces a "<>
    "basis that is orthonormal with respect to volume measure "<>
    "on the ball.  The value Sphere produces a "<>
    "basis that is orthonormal with respect to surface area "<>
    "measure on the sphere."

Partial::"usage" = 
  "Partial[f, \!\(x\_j\)] gives the partial derivative "<>
    "of f with respect to \!\(x\_j\)."

Phi::"usage"=
  "Phi[z] is the modified inversion defined in "<>
    "Chapter 7 of Harmonic Function Theory."

PoissonKernel::"usage"=
  "PoissonKernel[x, z] gives the Poisson kernel for the unit ball."

PoissonKernelH::"usage"=
  "PoissonKernelH[x, y, t] gives the "<>
    "Poisson kernel for the upper half-space."

S::"usage"="S is the south pole."

Schwarz::"usage"=
  "Schwarz[x] gives the maximum of |u[x]|, where "<>
    "u ranges over all harmonic functions on the unit ball with "<>
    "u[0] = 0 and |u| < 1."

SetDimension::"usage"="SetDimension[x, n] sets the Dimension of x to n."

Singularity::"usage"=
  "Singularity is an option for AntiLaplacian.  "<>
    "The default value is None."

Sphere::"usage"=
  "Sphere is one of the values that may be assigned to the option "<>
    "Orthonormal."

SurfaceArea::"usage"=
  "SurfaceArea[n] gives the surface area of "<>
    "the unit ball in n dimensional real Euclidean space."

Taylor::"usage"=
  "Taylor[f, d, x] gives the Taylor expansion of f "<>
    "to degree d, with respect to x."

Togetherness::"usage"=
  "TurnOff[ Togetherness ] turns off the feature "<>
    "that simplifies some output.  TurnOn[ Togetherness ] turns "<>
    "this feature back on."

Tr::"usage"="Tr[A] is the trace of a matrix A."

TurnOff::"usage"="TurnOff[F] turns feature F off."

TurnOn::"usage"="TurnOn[F] turns feature F on."

Volume::"usage"=
  "Volume[n] gives the volume of the unit ball in "<>
    "n-dimensional real Euclidean space."

ZonalHarmonic::"usage"=
  "ZonalHarmonic[m, x, z] gives the "<>
    "zonal harmonic of degree m with pole x."



Begin["`private`"]

Unprotect[Im,IntegerQ,Limit,Transpose]

protectedWords=Unprotect[Abs,Dot,IdentityMatrix,Integer,Power,Times]



SetAttributes[postTogether,HoldAll]

postTogether[Apart[f_]]:=Apart[f]

postTogether[Together[f_]]:=Together[f]

postTogether[Simplify[f_]]:=Simplify[f]

postTogether[Expand[f_]]:=Expand[f]

postTogether[a_=b_]:=a=postTogether[b]

shorter[f_] := With[ {g = Together[f]},
	If[ StringLength[ ToString[ Format[ g, InputForm ] ] ] < 
	StringLength[ ToString[ Format[ f, InputForm ] ] ], g, f ] ]

shorter/:f_Symbol[a___,HoldPattern[shorter][x_],b___]:=
  f[a,x,b]/;
    MemberQ[Attributes[f],HoldAll]||MemberQ[Attributes[f],HoldFirst]||
      MemberQ[Attributes[f],HoldRest]

postTogether[f_]:=Map[shorter,f,{0,-2}]//.HoldPattern[shorter][j_]\[Rule]j

usersPost=If[ToString[$Post]==="$Post",Identity,$Post]

Togetherness::Off=
  "Togetherness has been turned off.  "<>
    "The command  TurnOn[Togetherness]  will turn it back on."

Togetherness::On=
  "Togetherness has been turned on.  "<>
    "The command  TurnOff[Togetherness]  will turn it back off."

Togetherness::AlreadyOn="Togetherness is already turned on."

Togetherness::AlreadyOff="Togetherness is already turned off."

TurnOn[Togetherness]^:= 
  If[ togetherness === True, 
    Message[Togetherness::AlreadyOn],$Post=
      Composition[postTogether,usersPost];togetherness=True; 
    Message[Togetherness::On]]

$Post=Composition[postTogether,usersPost];

togetherness=True;

TurnOff[Togetherness]^:= 
  If[ togetherness === False, 
    Message[Togetherness::AlreadyOff],$Post=usersPost;togetherness=False;
    Message[Togetherness::Off]]





Off[IdentityMatrix::"intnm"]

matrix[IdentityMatrix[_]]^=True;

IdentityMatrix/:a_ .IdentityMatrix[_]:=a

IdentityMatrix/:IdentityMatrix[_].v_:=v



Off[Transpose::"nmtx"]

matrix[Transpose[_]]=True

Transpose[Transpose[x_]]:=x

Transpose[{c_ x_?vector}]:=c Transpose[{x}]

Transpose[c_ x_?matrix]:=c Transpose[x]

Transpose[IdentityMatrix[n_]]:=IdentityMatrix[n]

Transpose[{0}]:={0}



\!\(Tr[x_] := 
    \[Sum]\+\(j = 1\)\%\(Length[x]\)x\[LeftDoubleBracket]j, j
          \[RightDoubleBracket] /; 
      MatrixQ[x] && 
        Length[x] == Length[x\[LeftDoubleBracket]1\[RightDoubleBracket]]\)

Tr[IdentityMatrix[n_]]:=n

Tr[Transpose[A_]]:=Tr[A]

Tr[c_ A_]:=c Tr[A]/;matrix[A]



Norm=DoubleBracketingBar; Protect[Norm]

\!\(\[LeftDoubleBracketingBar]{x_?vector, y_}\[RightDoubleBracketingBar] := 
    \((unmakeVector[y]; 
      \@\(\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2 + y\^2\))
      \)\)

\[LeftDoubleBracketingBar]x_\[RightDoubleBracketingBar]:=
  0/;(makeVector[x];False)

\!\(\[LeftDoubleBracketingBar]x_List\[RightDoubleBracketingBar] := 
    \@\(x . x\)\)

\[LeftDoubleBracketingBar]0\[RightDoubleBracketingBar]=0;

\[LeftDoubleBracketingBar]c_?NumberQ x_\[RightDoubleBracketingBar]:=(
    makeVector[x];
    Abs[c] \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])

\[LeftDoubleBracketingBar]c_ x_?vector\[RightDoubleBracketingBar]:=(
    Unprotect[Im];Im[c]=0;
    Abs[c] \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])



\!\(ExpandNorm[f_] := 
    f //. \[LeftDoubleBracketingBar]g_ + h_\[RightDoubleBracketingBar] 
        \[RuleDelayed] 
        \@\(\[LeftDoubleBracketingBar]g\[RightDoubleBracketingBar]\^2 + 
            \[LeftDoubleBracketingBar]h\[RightDoubleBracketingBar]\^2 + 
            2\ g . h\)\)



\!\(Abs[\[LeftDoubleBracketingBar]x_\[RightDoubleBracketingBar]\^p_. ] := 
    \((makeVector[x]; 
      \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^p)\)\)

\!\(Abs /: Abs[x_\_j_]\^p_?EvenQ := x\_j\^p /; vector[x]\)

\!\(Abs /: Abs[c_]\^p_?EvenQ := 
    c\^p /; Im[c] == 0 || Im[1\/c] == 0 || Im[c\^\(p/2\)] == 0 || 
        Im[1\/c\^\(p/2\)] == 0\)



vector[HoldPattern[+a__]]:=
  If[TrueQ[Or@@Thread[vector[{a}]]],makeVector[+a];True]

vector[HoldPattern[Times[a__]]]:=Or@@Thread[vector[{a}]]

vector[A_?matrix.x_?vector]:=True

vector[x_?vector.Y_?matrix]:=True

vector[a_[".",_]]:=True

vector[a_[_,"."]]:=True

0[j_]:=0

\!\(vector[\(Partial\__[b_]\)[x_]] := True /; vector[b[x]]\)

vector[Grad[_][_]]=True



vectors={};

vectorValuedFunctions={};

\!\(makeVector[x_Symbol] := 
    If[MemberQ[vectors, x], Null, x[j_] = x\_j; AppendTo[vectors, x]; 
      vector[x] ^= True; ]\)

makeVector[c_?NumberQ x_]:=makeVector[x]

makeVector[HoldPattern[+a__]]:=makeVector/@{a}

\!\(\*
  RowBox[{\(makeVector[f_Symbol[_]]\), ":=", 
    RowBox[{
      RowBox[{"If", "[", 
        RowBox[{\(MemberQ[vectorValuedFunctions, f]\), ",", "Null", ",", 
          RowBox[{
          \(\(f[x_]\)[j_] = f\_j[x]\), ";", "\n", "\t\t", 
            \(AppendTo[vectorValuedFunctions, f]\), ";", 
            \(vector[f[_]] ^= True\), ";", 
            RowBox[{
              RowBox[{"vector", "[", 
                RowBox[{
                  SuperscriptBox["f", 
                    TagBox[\((_)\),
                      Derivative],
                    MultilineFunction->None], "[", "_", "]"}], "]"}], "=", 
              "True"}], ";", \(vector[\(Partial\__[f]\)[_]] = True\), 
            ";"}]}], "]"}], "/;", \(f =!= List\)}]}]\)

\!\(makeVector[\(Partial\__[u_]\)[_]] := makeVector[u[_]]\)

Attributes[unmakeVector]={Listable};

unmakeVector[x_Symbol]:=
  If[MemberQ[vectors,x],(x[j$_])=.;
    vectors=Delete[vectors,
        Position[vectors,x]\[LeftDoubleBracket]1,1\[RightDoubleBracket]];
    vector[x]^=False;]



matrix[HoldPattern[Times[a__]]]:=Or@@Thread[matrix[{a}]]

matrix[HoldPattern[+a__]]:=Or@@Thread[matrix[{a}]]

matrix[a_?matrix.x_?matrix]:=True

matrix[J[_][_]]=True;

matrix[{_}]=True;

0[_,_]:=0

matrix[Transpose[x_]]:=matrix[x]



matrices={};

\!\(makeMatrix[x_Symbol] := 
    If[TrueQ[matrix[x]], Null, x[j_, k_] = x\_\(j, k\); 
      AppendTo[matrices, x]; matrix[x] ^= True; vector[x] ^= False]\)



Attributes[Dot]={OneIdentity}

b_ . 0 = 0;

0 . b_ = 0;

HoldPattern[+x__].z_:=Plus@@(#1.z&)/@{x}

a_ .HoldPattern[+x__]:=Plus@@(a.#1&)/@{x}

(c_ x_).y_:=c x.y/;vector[x]||matrix[x]||NumberQ[c]

{c_ x_}.y_:=c {x}.y/;vector[x]||matrix[x]||NumberQ[c]

a_ .(c_ y_):=c a.y/;vector[y]||matrix[y]||NumberQ[c]

a_ .{c_ y_}:=c a.{y}/;vector[y]||matrix[y]||NumberQ[c]

\!\(x_\  . x_ := 
    \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2 /; 
      Head[x] == Symbol\  && vector[x]\)

\!\(f_[x_]\  . \ f_[x_]\  := \ 
    \[LeftDoubleBracketingBar]f[x]\[RightDoubleBracketingBar]\^2\  /; \ 
      Head[f] == Symbol\  && vector[f[x]]\)

\!\(\(\(Derivative[s_]\)[f_]\)[x_]\  . \ \(\(Derivative[s_]\)[f_]\)[x_]\  := 
    \ \[LeftDoubleBracketingBar]\(\(Derivative[s]\)[f]\)[x]\ 
          \[RightDoubleBracketingBar]\^2 /; \ 
      Head[f]\  == \ Symbol\  && vector[f[x]]\)

\!\(\(Partial\_j_[f_]\)[x_]\  . \ \(Partial\_j_[f_]\)[x_]\  := \ 
    \[LeftDoubleBracketingBar]\ 
          \(Partial\_j[f]\)[x]\[RightDoubleBracketingBar]\^2 /; \ 
      Head[f]\  == \ Symbol\  && vector[f[x]]\)

Delta[j_].Delta[k_]:=Delta[j,k]

\!\(Delta[j_] . x_Symbol := \((makeVector[x]; x\_j)\) /; matrix[x] =!= True\)

\!\(a_Symbol . Delta[j_] := \((makeVector[a]; a\_j)\) /; matrix[a] =!= True\)

Delta[j_].f_Symbol[x_]:=(makeVector[f[_]];f[x][j])/;matrix[f[x]]=!=True/;
    f=!=List

b_Symbol[x_].Delta[j_]:=(makeVector[b[_]];b[x][j])/;matrix[b[x]]=!=True/;
    b=!=List

\!\(Delta[j_] . \(Grad[u_]\)[x_] := \(Partial\_{j}[u]\)[x]\)

\!\(\(J[b_]\)[x_] . Delta[j_] := \(Partial\_{j}[b]\)[x]\)

Delta[j_].J[b_][x_]:=Grad[b.Delta[j]][x]

\!\(Delta[j_] . Y_?matrix := Y\_\(j, "\<.\>"\)\)

\!\(A_?matrix . Delta[j_] := A\_\("\<.\>", j\)\)

\!\(Delta[k_] . Y\_\(j_, "\<.\>"\) := Y\_\(j, k\)\)

\!\(A\_\(j_, "\<.\>"\)\  . Delta[k_] := A\_\(j, k\)\)

\!\(Delta[j_] . Y\_\("\<.\>", k_\) := Y\_\(j, k\)\)

\!\(A\_\("\<.\>", k_\)\  . Delta[j_] := A\_\(j, k\)\)

a_?VectorQ.Delta[j_]:=a\[LeftDoubleBracket]j\[RightDoubleBracket]

Delta[j_].v_?VectorQ:=v\[LeftDoubleBracket]j\[RightDoubleBracket]

A_?MatrixQ.Delta[j_]:=Transpose[A]\[LeftDoubleBracket]j\[RightDoubleBracket]

Delta[j_].Y_?MatrixQ:=A\[LeftDoubleBracket]j\[RightDoubleBracket]

x_ .A_?identityMatrixMultiple:=
  A\[LeftDoubleBracket]1,1\[RightDoubleBracket] x/;primitiveVector[x]

A_?identityMatrixMultiple.x_:=
  A\[LeftDoubleBracket]1,1\[RightDoubleBracket] x/;primitiveVector[x]

identityMatrixMultiple[A_]:=
  MatrixQ[A]&&Length[A]>1&&
    Length[A]===Length[A\[LeftDoubleBracket]1\[RightDoubleBracket]]&&
    And@@(A\[LeftDoubleBracket]#1,#1\[RightDoubleBracket]===
              A\[LeftDoubleBracket]1,1\[RightDoubleBracket]&)/@
        Range[Length[A]]&&
    And@@Flatten[
        Table[A\[LeftDoubleBracket]m,n\[RightDoubleBracket]===0||m===n,{m,
            Length[A]},{n,Length[A]}]]

\!\(a_List . x_Symbol := 
    \((makeVector[x]; a . Table[x\_j, {j, Length[a]}])\) /; 
      matrix[x] =!= True\)

\!\(x_Symbol . y_List := 
    \((makeVector[x]; Table[x\_j, {j, Length[y]}] . y)\) /; 
      matrix[x] =!= True\)

\!\(a_List . B_Symbol := 
    a . Table[\ B\_\(j, k\), \ {j, \ \(Dimension[B]\)[\([1]\)]}, 
          \ {k, \ \ \(Dimension[B]\)[\([2]\)]}] /; matrix[B]\)

\!\(B_Symbol . c_List := 
    Table[\ B\_\(j, k\), \ {j, \ \(Dimension[B]\)[\([1]\)]}, 
          \ {k, \ \ \(Dimension[B]\)[\([2]\)]}] . c /; matrix[B]\)

x_?primitiveVector.y_?primitiveVector:=
  y.x/;Order[ToString[y],ToString[x]]==1&&matrix[y]=!=True&&matrix[x]=!=True

primitiveVector[x_Symbol]:=True/;Context[x]=!="System`"

primitiveVector[Grad[_][_]]=True

primitiveVector[Delta[_]]=True



\!\(degree[f_, z_] := 
    degree1[Expand[
        f /. \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar] \[Rule] 
            z\_normReplacement], z]\)

degree1[0,_]:=-\[Infinity]

degree1[HoldPattern[+a__],z_]:=Max[Thread[degree1[{a},z]]]

\!\(degree1[f_, z_] := Exponent[f /. z\__ \[Rule] z, z]\)



\!\(monomialQ[f_, x_Symbol] := 
    PolynomialQ[f, Cases[Variables[x], x\__]] && FreeQ[f, Plus] && 
      \((Head[f] === Times && 
            \(\[InvisibleSpace]\(! 
              \((Or@@\(\((
                      FreeQ[f\[LeftDoubleBracket]#1\[RightDoubleBracket], x]&)
                      \)/@Range[Length[f]]\))\)\)\) || Head[f] === Power || 
          MatchQ[f, \ x\__])\)\)



polynomialDecomposition[0,x_]:=Null

\!\(polynomialDecomposition[p_, x_] := 
    With[{p1 = 
            Expand[p /. 
                \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] 
                  \[Rule] x\_normReplacement]}, 
        With[{t = 
              \(({#1, degree[#1, x]}&)\)/@
                If[Head[p1] === Plus, List@@p1, {p1}]}, 
          With[{s = 
                Table[If[
                    MemberQ[t, {_, j}], {
                      Plus@@\(Transpose[Cases[t, {_, j}]]
                            \)\[LeftDoubleBracket]1\[RightDoubleBracket], j}, 
                    Null], {j, 0, 
                    Max[\(Transpose[t]\)\[LeftDoubleBracket]2
                        \[RightDoubleBracket]]}]}, 
            Delete[s, Position[s, Null]]]]] /. 
      x\_normReplacement \[Rule] 
        \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\)



vector[Delta[_]]^:=True

Delta[j_,k_]:=0/;j\[NotEqual]k

Delta[j_,k_]:=1/;j==k

Delta[j_][k_]:=Delta[j,k]

\[LeftDoubleBracketingBar]Delta[_]\[RightDoubleBracketingBar]^=1;



\!\(Partial[f_, {x_\__, 0}] := \((makeVector[x]; f)\)\)

Partial[f_,{y_Symbol,0}]:=(unmakeVector[y];f)

Partial[f_,{z_,n_Integer}]:=
  Partial[postTogether[Partial[f,{z,n-1}]],z]/;Positive[n]

Partial[f_,z1_,z__]:=Partial[Partial[f,z1],z]

\!\(Partial[f_, x_\__] := \((makeVector[x]; 0)\) /; FreeQ[f, x]\)

Partial[f_,y_Symbol]:=(unmakeVector[y];0/;FreeQ[f,y])

\!\(Partial[x_\_j_, x_\_k_] := \((makeVector[x]; Delta[j, k])\)\)

Partial[y_,y_]:=1

\!\(Partial[x_, x_\_j_] := \((makeVector[x]; Delta[j])\)\)

Partial[HoldPattern[+a__],y_]:=Plus@@Thread[Partial[{a},y]]

Partial[f_ g_,y_]:=f Partial[g,y]+Partial[f,y] g

Partial[f_ .g_,y_]:=(If[TrueQ[matrix[f]],Null,makeVector[f]];
    If[TrueQ[matrix[g]],Null,makeVector[g]];f.Partial[g,y]+Partial[f,y].g)

\!\(Partial[f_[v_List], t_] := 
    Plus@@\(\((
          \(Partial\_{#1}[f]\)[v]\ *\ 
              Partial[v\[LeftDoubleBracket]#1\[RightDoubleBracket], t]&)\)/@
        Range[Length[v]]\)\)

\!\(Partial[f_\^p_, y_] := 
    p\ f\^\(p - 1\)\ Partial[f, y] + f\^p\ Log[f]\ Partial[p, y]\)

\!\(\*
  RowBox[{\(Partial[f_[g_], x_\_j_]\), ":=", 
    RowBox[{"(", 
      RowBox[{\(makeVector[x]\), ";", 
        RowBox[{"If", "[", 
          RowBox[{
          \(TrueQ[vector[g]]\), ",", 
            \(If[TrueQ[vector[f[g]]], \(J[f]\)[g] . Partial[g, x\_j], 
              \(Grad[f]\)[g] . Partial[g, x\_j]]\), ",", 
            RowBox[{
              RowBox[{
                SuperscriptBox["f", "\[Prime]",
                  MultilineFunction->None], "[", "g", "]"}], " ", 
              \(Partial[g, x\_j]\)}]}], "]"}]}], ")"}]}]\)

\!\(\*
  RowBox[{\(Partial[f_[g_], y_Symbol]\), ":=", 
    RowBox[{"(", 
      RowBox[{\(unmakeVector[y]\), ";", 
        RowBox[{"If", "[", 
          RowBox[{
          \(TrueQ[vector[g]]\), ",", 
            \(If[TrueQ[vector[f[g]]], \(J[f]\)[g] . Partial[g, y], 
              \(Grad[f]\)[g] . Partial[g, y]]\), ",", 
            RowBox[{
              RowBox[{
                SuperscriptBox["f", "\[Prime]",
                  MultilineFunction->None], "[", "g", "]"}], " ", 
              \(Partial[g, y]\)}]}], "]"}]}], ")"}]}]\)

\!\(Partial\_j_List[Partial\_k_List[f_]] := Partial\_\(Sort[Join[k, j]]\)[f]\)



NormalD[f_,z_Symbol]:=
  Grad[f,z].z/.\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\[Rule]1

\!\(NormalD[f_, v_List] := 
    \((unmakeVector[v]; 
      Together[Grad[f, v] . v] /. 
        \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]\^2 \[Rule] 1)
      \)\)



Grad[f_,x_Symbol]:=(makeVector[x];0)/;FreeQ[f,x]

\!\(Grad[x_\_j_, x_Symbol] := \((makeVector[x]; Delta[j])\)\)

Grad[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[Grad[{a},x]]

Grad[f_ g_,x_Symbol]:=f Grad[g,x]+Grad[f,x] g

Grad[f_ .g_,x_Symbol]:=f.J[g,x]+g.J[f,x]

\!\(Grad[f_\^p_, x_Symbol] := 
    p\ f\^\(p - 1\)\ Grad[f, x] + f\^p\ Log[f]\ Grad[p, x]\)

\!\(\*
  RowBox[{\(Grad[f_[y_], x_Symbol]\), ":=", 
    RowBox[{"(", 
      RowBox[{\(makeVector[x]\), ";", 
        RowBox[{"If", "[", 
          RowBox[{
          \(TrueQ[vector[y]]\), ",", \(\(Grad[f]\)[y] . J[y, x]\), ",", 
            RowBox[{
              RowBox[{
                SuperscriptBox["f", "\[Prime]",
                  MultilineFunction->None], "[", "y", "]"}], " ", 
              \(Grad[y, x]\)}]}], "]"}]}], ")"}]}]\)

\!\(Grad[DoubleBracketingBar] := #1
        \/\[LeftDoubleBracketingBar]#1\[RightDoubleBracketingBar]&\)

Grad[Dimension]:=0&

Grad[f_,v_List]:=(unmakeVector[v];(Partial[f,#1]&)/@v)



Divergence[f_,x_Symbol]:=(makeVector[x];0)/;FreeQ[f,x]

Divergence[x_,x_Symbol]:=(makeVector[x];Dimension[x])

Divergence[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[Divergence[{a},x]]

Divergence[f_ g_,x_Symbol]:=
  f Divergence[g,x]+Grad[f,x].g/;(makeVector[x];vector[g])

Divergence[A_?matrix.f_,x_Symbol]:=Tr[A.J[f,x]]/;FreeQ[A,x]

Divergence[f_ .T_?matrix,x_Symbol]:=Tr[Transpose[T].J[f,x]]/;FreeQ[T,x]

Divergence[u_List,v_List]:=(unmakeVector[v];
    Plus@@Thread[partial[u,v]]/.partial\[Rule]Partial)



J[f_,x_Symbol]:=(makeVector[f];makeVector[x];0)/;FreeQ[f,x]

J[x_,x_Symbol]:=(makeVector[x];IdentityMatrix[Dimension[x]])

J[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[J[{a},x]]

J[f_ g_?vector,x_Symbol]:=f J[g,x]+Transpose[{Grad[f,x]}].{g}

J[f_Symbol[x_],x_Symbol]:=(makeVector[x];makeVector[f[_]];J[f][x])/;f=!=List

\!\(\*
  RowBox[{\(J[f_[y_], x_Symbol]\), ":=", 
    RowBox[{"If", "[", 
      RowBox[{\(TrueQ[vector[y]]\), ",", \(\(J[f]\)[y] . J[y, x]\), ",", 
        RowBox[{
          RowBox[{"Transpose", "[", 
            RowBox[{"{", 
              RowBox[{
                SuperscriptBox["f", "\[Prime]",
                  MultilineFunction->None], "[", "y", "]"}], "}"}], "]"}], 
          ".", \({Grad[y, x]}\)}]}], "]"}]}]\)

J[A_?matrix.f_,x_Symbol]:=(makeVector[f];A.J[f,x])/;FreeQ[A,x]

J[f_ .T_?matrix,x_Symbol]:=(makeVector[f];Transpose[T].J[f,x])/;FreeQ[T,x]

J[f_List,v_List]:=(unmakeVector[v];
    Table[Partial[f\[LeftDoubleBracket]m\[RightDoubleBracket],
        v\[LeftDoubleBracket]n\[RightDoubleBracket]],{m,Length[f]},{n,
        Length[v]}])



\!\(HilbertSchmidt[A_?MatrixQ] := Plus@@\(Plus@@\(A\^2\)\)\)



Laplacian[f_,v_List]:=(unmakeVector[v];Plus@@(Partial[f,{#1,2}]&)/@v)

Laplacian[f_,x_Symbol]:=(makeVector[x];0)/;FreeQ[f,x]

\!\(Laplacian[x_\_j_, x_Symbol] := \((makeVector[x]; 0)\)\)

Laplacian[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[Laplacian[{a},x]]

Laplacian[f_ g_,x_Symbol]:=
  f Laplacian[g,x]+Laplacian[f,x] g+2 Grad[f,x].Grad[g,x]

\!\(Laplacian[f_\^p_, x_Symbol] := 
    p\ f\^\(p - 1\)\ Laplacian[f, x] + 
      p\ \((p - 1)\)\ f\^\(p - 2\)\ Grad[f, x]\  . \ Grad[f, x]\  + 
      2\ p\ f\^\(p - 1\)\ Log[f]\ Grad[f, x] . Grad[p, x] + 
      f\^p\ Log[f]\ Laplacian[p, x] + 
      2\ f\^\(p - 1\)\ Grad[f, x] . Grad[p, x] + 
      f\^p\ Log[f]\^2\ Grad[p, x]\  . \ Grad[p, x]\)

\!\(Laplacian[f_[a_. \ x_ + b_. ], x_Symbol] := 
    \((makeVector[x]; If[b == 0, Null, makeVector[b]]; 
        a\^2\ \(Laplacian[f]\)[a\ x + b])\) /; FreeQ[a, x] && FreeQ[b, x]\)

Laplacian[a_ .x_,x_Symbol]:=(makeVector[x];0)/;FreeQ[a,x]

Laplacian[x_ .y_,x_Symbol]:=(makeVector[x];0)/;FreeQ[y,x]

Laplacian[(A_ .x_).y_,x_Symbol]:=(makeVector[x];0)/;FreeQ[A,x]&&FreeQ[y,x]

Laplacian[(x_ .Y_).z_,x_Symbol]:=(makeVector[x];0)/;FreeQ[Y,x]&&FreeQ[z,x]

Laplacian[a_ .(B_ .x_),x_Symbol]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[B,x]

Laplacian[a_ .(x_ .Y_),x_Symbol]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[Y,x]

Laplacian[(A_ .x_).x_,x_Symbol]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x]

Laplacian[(x_ .A_).x_,x_Symbol]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x]

Laplacian[x_ .(A_ .x_),x_Symbol]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x]

Laplacian[x_ .(x_ .A_),x_Symbol]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x]

Laplacian[b_ .Grad[u_][b_],b_Symbol]:=(makeVector[b];
    2 Laplacian[u[b],b]+b.Grad[Laplacian[u[b],b],b])

Laplacian[Grad[u_][x_].x_,x_Symbol]:=(makeVector[x];
    2 Laplacian[u[x],x]+x.Grad[Laplacian[u[x],x],x])

\!\(\*
  RowBox[{\(Laplacian[f_[g_], x_Symbol]\), ":=", 
    RowBox[{
      RowBox[{
        RowBox[{
          RowBox[{
            SuperscriptBox["f", "\[DoublePrime]",
              MultilineFunction->None], "[", "g", "]"}], "  ", 
          \(Grad[g, x]\  . \ Grad[g, x]\)}], " ", "+", " ", 
        RowBox[{
          RowBox[{
            SuperscriptBox["f", "\[Prime]",
              MultilineFunction->None], "[", "g", "]"}], " ", 
          \(Laplacian[g, x]\)}]}], "/;", 
      \((makeVector[x]; vector[g] =!= True)\)}]}]\)

\!\(Laplacian[DoubleBracketingBar] := 
    \(Dimension[#1] - 1\)\/\[LeftDoubleBracketingBar]#1
          \[RightDoubleBracketingBar]&\)

\!\(Laplacian[\[LeftDoubleBracketingBar]A_\  . x_\[RightDoubleBracketingBar], 
      x_] := \((makeVector[x]; 
        \(\[LeftDoubleBracketingBar]A . x\[RightDoubleBracketingBar]\^2\ 
              HilbertSchmidt[A]\^2 - 
            \[LeftDoubleBracketingBar]A . x . A\[RightDoubleBracketingBar]\^2
              \)\/\[LeftDoubleBracketingBar]A . x
                \[RightDoubleBracketingBar]\^3)\) /; FreeQ[A, x]\)

\!\(Laplacian[\[LeftDoubleBracketingBar]x_\  . Z_\[RightDoubleBracketingBar], 
      x_] := \((makeVector[x]; 
        \(\[LeftDoubleBracketingBar]x . Z\[RightDoubleBracketingBar]\^2\ 
              HilbertSchmidt[Z]\^2 - 
            \[LeftDoubleBracketingBar]x . Z . 
                  Transpose[Z]\[RightDoubleBracketingBar]\^2
              \)\/\[LeftDoubleBracketingBar]x . Z
                \[RightDoubleBracketingBar]\^3)\) /; FreeQ[Z, x]\)

Laplacian[Dimension]:=0&

Laplacian[f_,x_Symbol,y_Symbol]:=(unmakeVector[y];
    Laplacian[f,x]+Partial[f,{y,2}])

\!\(\(Laplacian\^n_\)[f_, x___] := 
    \(Laplacian\^\(n - 1\)\)[Laplacian[f, x], x] /; 
      Head[n] == Integer && Positive[n]\)

\!\(Laplacian[\(Laplacian\^n_. \)[f_]] := \(Laplacian\^\(n + 1\)\)[f]\)

\!\(\(Laplacian\^m_. \)[\(Laplacian\^n_. \)[u_]] := 
    \(Laplacian\^\(m + n\)\)[u]\)



Dimension[a_ x_]:=Dimension[x]/;vector[x]||matrix[x]

Dimension[a_?NumberQ x_]:=(If[matrix[x]=!=True,makeVector[x]];Dimension[x])

Dimension[a_. S+x_]:=Dimension[x]

Dimension[x_+w_]:=
  If[StringLength[ToString[Dimension[x]]]<
      StringLength[ToString[Dimension[w]]],Dimension[x],Dimension[w]]

Dimension[v_?VectorQ]:=Length[v]

Dimension[A_?MatrixQ]:={Length[A],
    Length[A\[LeftDoubleBracket]1\[RightDoubleBracket]]}

IntegerQ[Dimension[_]]^:=True



SetDimension::Matrix = 
  "`1` will be considered to be a matrix of size `2`-by-`3`."

SetDimension[x_Symbol,{m_,n_}]:=(Dimension[x]={m,n};makeMatrix[x];
    If[Not[IntegerQ[m]],Unprotect[IntegerQ];IntegerQ[m]=True];
    If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];
			With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, 
      SetOptions[$FrontEnd, 
        MessageOptions\[Rule]{Global`WarningAction\[Rule]{}}];
      Message[SetDimension::Matrix, x, m, n]; 
      SetOptions[$FrontEnd, messageOptions];])

SetDimension::Vector = 
  "`1` will be considered to be a vector in `2`-dimensional real Euclidean \
space."

setDimension[x_Symbol,n_]:=(Dimension[x]=n;makeVector[x];
    If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];)

SetDimension[x_Symbol,n_]:=(setDimension[x, n];
    With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, 
      SetOptions[$FrontEnd, 
        MessageOptions\[Rule]{Global`WarningAction\[Rule]{}}];
      Message[SetDimension::Vector, x, n]; 
      SetOptions[$FrontEnd, messageOptions];])

SetDimension::VectorValuedFunction = 
  "`1` will be considered to be a function taking values in `2`-dimensional \
Euclidean space."

SetDimension[f_Symbol[_],
    n_]:=(Dimension[f[_]]=n;makeVector[f[_]];
      If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];
		With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, 
        SetOptions[$FrontEnd, 
          MessageOptions\[Rule]{Global`WarningAction\[Rule]{}}];
        Message[SetDimension::VectorValuedFunction, f, n]; 
        SetOptions[$FrontEnd, messageOptions];])/;f=!=List



\!\(Homogeneous[f_, d_, x_Symbol, c_Symbol] := 
    \((makeVector[x]; makeVector[c]; setDimension[c, Dimension[x]]; 
        Homogeneous[f, d, Union[Cases[Level[f, \(-1\)], x\__]], 
            Union[Cases[Level[f, \(-1\)], x\__]] /. x \[Rule] c] /. 
          Table[\ c\_j, 
              \ {j, \ Length[Union[Cases[Level[f, \(-1\)], x\__]]]}] \[Rule] 
            c)\) /; \(! \((MemberQ[Level[f, {\(-1\)}], x])\)\) && 
        Head[d] == Integer\)

\!\(Homogeneous[f_, d_, x_Symbol, c_Symbol] := 
    \((makeVector[x]; makeVector[c]; setDimension[c, Dimension[x]]; 
        Homogeneous[f, d, x, Table[c\_j, \ {j, \ Dimension[x]}]] /. 
          Table[c\_j, \ {j, \ Dimension[x]}] \[Rule] c)\) /; 
      MemberQ[Level[f, {\(-1\)}], x] && Head[Dimension[x]] == Integer && 
        Head[d] == Integer\)

\!\(Homogeneous[f_, d_, x_Symbol, c_List] := 
    \((makeVector[x]; unmakeVector[c]; 
        Module[{z}, 
          \(\(Homogeneous[f, d, Table[x\_j, \ {j, Length[c]}], c] /. 
                x\_j_ \[Rule] z\_j\) /. x \[Rule] c\) /. z\_j_ \[Rule] x\_j])
        \) /; Head[d] == Integer\)

\!\(Homogeneous[f_, d_, v_List, c_Symbol] := 
    \((makeVector[c]; unmakeVector[v]; setDimension[c, Length[v]]; 
        Homogeneous[f, d, v, Table[c\_j, \ {j, \ Dimension[c]}]])\) /; 
      Head[d] == Integer\)

\!\(Homogeneous[f_, d_, v_List, c_List] := 
    \((unmakeVector[v]; unmakeVector[c]; 
        If[togetherness && 
            \(! \((\[LeftDoubleBracketingBar]c\[RightDoubleBracketingBar] === 
                  0)\)\), TurnOff[Togetherness]]; 
        With[{mI = multiIndices[d, Length[v]]}, 
          Plus@@\(\((
                \((\((Partial@@
                              Prepend[
                                Transpose[{v, 
                                    mI\[LeftDoubleBracket]#1
                                      \[RightDoubleBracket]}], f] /. 
                            Thread[v \[Rule] c])\)\ 
                        Times@@
                          \(\((v - c)\)\^mI\[LeftDoubleBracket]#1
                                \[RightDoubleBracket]\))\)/
                    multiFactorial[
                      mI\[LeftDoubleBracket]#1\[RightDoubleBracket]]&)\)/@
              Range[Length[mI]]\)])\) /; Head[d] == Integer\)

Homogeneous[f_,d_,v_List]:=
  Homogeneous[f,d,v,Table[0,{Length[v]}]]/;Head[d]==Integer

\!\(Homogeneous[f_, d_, x_Symbol] := 
    \((Homogeneous[f, d, Union[Cases[Level[f, \(-1\)], x\__]]] /. 
          Table[0, {Length[Union[Cases[Level[f, \(-1\)], x\__]]]}] \[Rule] 0)
        \) /; \(! \((MemberQ[Level[f, {\(-1\)}], x])\)\) && 
        Head[d] == Integer\)

Homogeneous[f_,d_,
    x_Symbol]:=(
      Homogeneous[f,d,x,
          Table[0,{Dimension[x]}]]/.Table[0,{Dimension[x]}]\[Rule]0)/;
    MemberQ[Level[f,{-1}],x]&&Head[Dimension[x]]==Integer&&Head[d]==Integer

multiFactorial[a_List]:=Times@@(a!)

multiIndices[k_,1]:={{k}}

multiIndices[0,n_]:={Table[0,{n}]}

multiIndices[k_Integer,n_Integer]:=
  multiIndices[k,n]=
    Join@@Table[(Prepend[#1,k-j]&)/@multiIndices[j,n-1],{j,0,k}]



Taylor[f_,d_,vc___]:=
  Plus@@(Homogeneous[f,#1,vc]&)/@Range[0,d]/;Head[d]==Integer



\!\(Volume[n_Integer] := \[Pi]\^\(n/2\)\/Gamma[n\/2 + 1] /; Positive[n]\)



\!\(SurfaceArea[n_Integer] := 
    \(2\ \[Pi]\^\(n/2\)\)\/Gamma[n\/2] /; Positive[n]\)



\!\(IntegrateBall[f_, v_List] := 
    \((unmakeVector[v]; 
      Module[{z}, setDimension[z, Length[v]]; 
        IntegrateBall[f /. Thread[v \[Rule] Table[z\_j, \ {j, \ Length[v]}]], 
          z]])\)\)

IntegrateBall[f_,x_Symbol]:=(makeVector[x];
		f Volume[Dimension[x]])/;FreeQ[f,x]

\!\(IntegrateBall[f_, x_Symbol] := 
    \[Integral]\_\(-1\)\%1
          \((f /. {
                \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2 
                  \[Rule] x\^2, 
                \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] 
                  \[Rule] Abs[x], x\_1 \[Rule] x})\) \[DifferentialD]x /; 
      Dimension[x] == 1\)

\!\(IntegrateBall[f_, x_Symbol] := 
    \((makeVector[x]; 
      Module[{r}, 
        If[Head[Dimension[x]] === Integer, Null, 
          Limit[c_. \ r\^\(Dimension[x] + m_. \), r \[Rule] 0] := 
            0 /; \((m > 0 || m == 0)\) && FreeQ[c, r]]; 
        Dimension[x]\ Volume[Dimension[x]]\ 
          \(\[Integral]\_0
              \%1 Expand[
                r\^\(Dimension[x] - 1\)\ 
                  IntegrateSphere[
                    f /. {\[LeftDoubleBracketingBar]x
                            \[RightDoubleBracketingBar] \[Rule] r, 
                        x\_j_ \[Rule] r\ x\_j, x \[Rule] r\ x}, x]] 
              \[DifferentialD]r\)])\)\)



\!\(IntegrateSphere[f_, v_List] := 
    \((unmakeVector[v]; 
      Module[{z}, setDimension[z, Length[v]]; 
        IntegrateSphere[
          f /. Thread[v \[Rule] Table[\ z\_j, \ {j, Length[v]}]], z]])\)\)

\!\(IntegrateSphere[f_, x_Symbol] := 
    \((makeVector[x]; 
      integrateSphere1[
        Expand[f /. {
              \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] \[Rule] 
                1, x . z_ \[RuleDelayed] 
                Table[x\_j, {j, \ Dimension[x]}] . 
                  Table[z\_j, {j, \ Dimension[x]}], 
              z_\  . x \[RuleDelayed] 
                Table[z\_j, {j, \ Dimension[x]}] . 
                  Table[x\_j, {j, \ Dimension[x]}]}], x])\)\)

integrateSphere1[HoldPattern[+a__],x_]:=Plus@@Thread[integrateSphere1[{a},x]]

integrateSphere1[f_,x_]:=f/;FreeQ[f,x]

\!\(integrateSphere1[c_. \ x_\__, x_] := 0\)

integrateSphere1[c_ f_,x_]:=c integrateSphere1[f,x]/;FreeQ[c,x]

\!\(integrateSphere1[x_\__\^j_, x_] := integrateSphere3[{j}, Dimension[x]]\)

integrateSphere1[f_,x_]:=integrateSphere2[f,x]/;monomialSpecialQ[f,x]

\!\(integrateSphere2[f_, x_] := 
    integrateSphere3[List@@f /. x\__\^b_ \[Rule] b, Dimension[x]]\)

integrateSphere3[b_,n_]:=0/;Or@@OddQ[b]

\!\(integrateSphere3[b_, n_] := 
    Times@@\((\(\((b - 1)\)!!\))\)\/Times@@Range[n, n + Plus@@b - 2, 2]\)

\!\(monomialSpecialQ[f_, x_Symbol] := 
    Head[f] === Times && 
      And@@\(\((
            MatchQ[f\[LeftDoubleBracket]#1\[RightDoubleBracket], x\__\^_]&)
            \)/@Range[Length[f]]\)\)





Dirichlet[p_,x_Symbol]:=
  Plus@@(#1\[LeftDoubleBracket]1\[RightDoubleBracket]&)/@
      HarmonicDecomposition[
        Limit[p/.
            \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\[Rule]r,
          r\[Rule]1],x]

Dirichlet[f_,g_,z_Symbol]:=(makeVector[z];
    With[{t=AntiLaplacian[g,z]},t+Dirichlet[f-t,z]])

Dirichlet[f_,v_List]:=explicit[Dirichlet,f,v]

Dirichlet[f_,g_,v_List]:=explicit[Dirichlet,f,g,v]



d[0,0]=1

\!\(lap[f_, v_] := Plus@@\(\((\[PartialD]\_{#1, 2}f&)\)/@v\)\)

HarmonicDecomposition[0,x_]:={}

\!\(HarmonicDecomposition[p_, x_Symbol] := 
    \((makeVector[x]; 
      With[{deg = 
            degree[p /. 
                \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] 
                  \[Rule] norm, x], n = Dimension[x]}, powers = {}; 
        hd[_] := 0; 
        With[{partition = 
              Partition[
                Flatten[
                  Expand[\((
                          \((m = #1\[LeftDoubleBracket]2
                                  \[RightDoubleBracket]; k = Floor[m\/2]; 
                              If[m > 1, 
                                d[1, 1] = 1\/\(2\ \((2\ m + n - 4)\)\)]; 
                              v = Select[
                                  Variables[#1\[LeftDoubleBracket]1
                                      \[RightDoubleBracket]], 
                                  \(! \((FreeQ[#1, x])\)\)&]; 
                              laplac[0] = #1\[LeftDoubleBracket]1
                                  \[RightDoubleBracket]; 
                              laplac[1] = 
                                lap[#1\[LeftDoubleBracket]1
                                    \[RightDoubleBracket], v]; 
                              Do[d[j, j] = 
                                  \(d[j - 1, j - 1]\ 
                                      \((2\ m + n - 2\ j)\)\)\/\(2\ j\ 
                                      \((2\ m + n - 4\ j + 2)\)\ 
                                      \((2\ m + n - 4\ j)\)\); 
                                laplac[j] = lap[laplac[j - 1], v], {j, 2, 
                                  k}]; Do[
                                d[i, j] = 
                                  d[i, j - 1]\/\(2\ \((i - j)\)\ 
                                      \((2\ \((m - i - j - 1)\) + n)\)\), {i, 
                                  0, k - 1}, {j, i + 1, k}]; 
                              Table[{
                                  \[Sum]\+\(j = i\)\%k d[i, j]\ 
                                      \[LeftDoubleBracketingBar]x
                                        \[RightDoubleBracketingBar]\^\(2\ 
                                        \((j - i)\)\)\ laplac[j], 2\ i}, {i, 
                                  0, k}])\)&)\)/@
                        polynomialDecomposition[
                          p /. \[LeftDoubleBracketingBar]x
                                \[RightDoubleBracketingBar] \[Rule] norm, 
                          x]] //. {a_.  + norm\^e_. \ b_. , j_} 
                      \[Rule] {{a, j}, {b, j + e}}], 2]}, 
          powers = 
            powers \[Union] 
              \(Transpose[partition]\)\[LeftDoubleBracket]2
                \[RightDoubleBracket]; 
          \((\((hd[#1\[LeftDoubleBracket]2\[RightDoubleBracket]] = 
                    hd[#1\[LeftDoubleBracket]2\[RightDoubleBracket]] + #1
                        \[LeftDoubleBracket]1\[RightDoubleBracket])\)&)\)/@
            partition]; 
        aa = Select[
            \(({hd[#1], #1}&)\)/@
              powers, #1\[LeftDoubleBracket]1\[RightDoubleBracket] =!= 0&]; 
        Clear[hd]; aa])\)\)

HarmonicDecomposition[f_,v_List]:=explicit[HarmonicDecomposition,f,v]



Options[AntiLaplacian]={Multiple\[Rule]None,Singularity\[Rule]None}

AntiLaplacian[0,x_,options___?OptionQ]:=0

\!\(AntiLaplacian[p_, x_Symbol, options___?OptionQ] := 
    \((makeVector[x]; 
      If[\((\(Multiple /. {options}\) /. Options[AntiLaplacian])\) === 
            Norm\^2 || 
          \(! \((FreeQ[p, 
                \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]])\)\), 
        antiLaplacianN[Expand[p], x, options], antiLaplacian[Expand[p], x]])
      \)\)

antiLaplacianN[HoldPattern[+a__],z_,options___?OptionQ]:=
  Plus@@Thread[Unevaluated[antiLaplacianN[{a},z,options]]]

antiLaplacianN[c_ g_,z_,options___?OptionQ]:=
  c antiLaplacianN[g,z,options]/;FreeQ[c,z]

antiLaplacianN[g_,z_,options___?OptionQ]:=
  antiLap1[If[Head[g]===Times,
      Select[g,FreeQ[#1,
            \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]]&],
      If[FreeQ[g,\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]],g,
        1]],If[Head[g]===Times,
      Select[g,!(
              FreeQ[#1,
                \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]])&],
      If[FreeQ[g,\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]],1,
        g]],z,options]

antiLap1[p_,f_,z_,options___?OptionQ]:=
  Plus@@(antiLap2[#1\[LeftDoubleBracket]1\[RightDoubleBracket],#1
              \[LeftDoubleBracket]2\[RightDoubleBracket],f,z,options]&)/@
      HarmonicDecomposition[p,z]

antiLap2[p_,j_,f_,z_,options___?OptionQ]:=
  Plus@@(antiLap3[#1\[LeftDoubleBracket]1\[RightDoubleBracket],#1
              \[LeftDoubleBracket]2\[RightDoubleBracket],j,f,z,options]&)/@
      polynomialDecomposition[p,z]

\!\(antiLap3[p_, m_, j_, f_, z_, options___?OptionQ] := 
    With[{n = Dimension[z]}, 
      p\ If[\((\(Singularity /. {options}\) /. Options[AntiLaplacian])\) =!= 
              None, \[Integral]t\^\(1 - 2\ m - n\)\ 
                \((\[Integral]\((
                        t\^\(\(-1\) + 2\ m + n + j\)\ f /. 
                          \[LeftDoubleBracketingBar]z
                              \[RightDoubleBracketingBar] \[Rule] t)\) 
                      \[DifferentialD]t)\) \[DifferentialD]t, 
            With[{int = 
                  \[Integral]\_0\%t
                      \((s\^\(\(-1\) + 2\ m + n + j\)\ f /. 
                          \[LeftDoubleBracketingBar]z
                              \[RightDoubleBracketingBar] \[Rule] s)\) 
                      \[DifferentialD]s}, 
              \[Integral]t\^\(1 - 2\ m - n\)\ 
                  If[FreeQ[int, Indeterminate] && FreeQ[int, \[Infinity]] && 
                      FreeQ[int, Limit] && FreeQ[int, If] && 
                      FreeQ[int, Integrate] && 
                      \((\(Singularity /. {options}\) /. 
                            Options[AntiLaplacian])\) === None, int, 
                    \[Integral]\((
                          t\^\(\(-1\) + 2\ m + n + j\)\ f /. 
                            \[LeftDoubleBracketingBar]z
                                \[RightDoubleBracketingBar] \[Rule] t)\) 
                        \[DifferentialD]t] \[DifferentialD]t]] /. 
        t \[Rule] \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]]\)

\!\(AntiLaplacian[f_, v_List, options___?OptionQ] := 
    \((unmakeVector[v]; 
      Module[{z}, setDimension[z, Length[v]]; 
        AntiLaplacian[
            \(f /. \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]\^2 
                  \[Rule] 
                  \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2
                    \) /. Thread[v \[Rule] Table[z\_j, \ {j, \ Length[v]}]], 
            z, options] /. 
          Prepend[Thread[Table[z\_j, \ {j, \ Length[v]}] \[Rule] v], 
              \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]] \[Rule] 
            \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]])\)\)

antiLaplacian[HoldPattern[+a__],z_]:=Plus@@Thread[antiLaplacian[{a},z]]

\!\(antiLaplacian[g_, z_] := 1\/2\ g\ z\_1\^2 /; FreeQ[g, z]\)

\!\(antiLaplacian[g_, z_] := 
    antiLaplacian3[
        Transpose[
          Sort[\(({Exponent[g, z\_#1], #1}&)\)/@
              \(\((#1\[LeftDoubleBracket]2\[RightDoubleBracket]&)\)/@
                Variables[g]\)]], g, z] /; monomialQ[g, z]\)

antiLaplacian[c_ g_,z_]:=c antiLaplacian[g,z]/;FreeQ[c,z]

\!\(antiLaplacian3[{e_, p_}, g_, z_] := 
    z\_\(p\[LeftDoubleBracket]\(-1\)\[RightDoubleBracket]\)\^\(e
              \[LeftDoubleBracket]\(-1\)\[RightDoubleBracket] + 2\)\/\(\((
            e\[LeftDoubleBracket]\(-1\)\[RightDoubleBracket] + 2)\)\ 
          \((e\[LeftDoubleBracket]\(-1\)\[RightDoubleBracket] + 1)\)\) /; 
      Length[e] == 1\)

\!\(antiLaplacian3[{e_, p_}, g_, z_] := 
    \(g\ z\_\(Last[p]\)\^2\)\/\(\((Last[e] + 1)\)\ \((Last[e] + 2)\)\) - 
        Plus@@\(\((
              \(\((e\[LeftDoubleBracket]#1\[RightDoubleBracket] - 1)\)\ 
                    e\[LeftDoubleBracket]#1\[RightDoubleBracket]\ 
                    antiLaplacian[
                      \(z\_\(Last[p]\)\^2\ g
                          \)\/z\_\(p\[LeftDoubleBracket]#1
                              \[RightDoubleBracket]\)\^2, z]\)\/\(\((
                      Last[e] + 1)\)\ \((Last[e] + 2)\)\)&)\)/@
            Range[\(Position[Positive[e], True]\)\[LeftDoubleBracket]1, 1
                \[RightDoubleBracket], Length[e] - 1]\) /; Length[e] > 1\)



ExteriorDirichlet[f_,z_]:=Factor[Kelvin[Dirichlet[f,z],z]]

ExteriorDirichlet[f_,g_,z_Symbol]:=
  Factor[With[{t=AntiLaplacian[g,z,Singularity\[Rule]0]},
      t+ExteriorDirichlet[f-t,z]]]

ExteriorDirichlet[f_,v_List]:=explicit[ExteriorDirichlet,f,v]

ExteriorDirichlet[f_,g_,v_List]:=explicit[ExteriorDirichlet,f,g,v]



annular[0,r_,s_,z_]=0

\!\(annular[f_, r_, s_, z_] := 
    Plus@@\(\((
            If[#1\[LeftDoubleBracket]2\[RightDoubleBracket] === #1
                    \[LeftDoubleBracket]3\[RightDoubleBracket], 
                \(\((Log[\[LeftDoubleBracketingBar]z
                            \[RightDoubleBracketingBar]] - Log[s])\)\ 
                    r\^#1\[LeftDoubleBracket]2\[RightDoubleBracket]\ #1
                      \[LeftDoubleBracket]1\[RightDoubleBracket]\)\/\(Log[
                      r] - Log[s]\), 
                \(\((\[LeftDoubleBracketingBar]z
                            \[RightDoubleBracketingBar]\^\(2
                              \ #1\[LeftDoubleBracket]2
                                \[RightDoubleBracket] - 
                            2\ #1\[LeftDoubleBracket]3
                                \[RightDoubleBracket]\) - 
                        s\^\(2\ #1\[LeftDoubleBracket]2
                                \[RightDoubleBracket] - 
                            2\ #1\[LeftDoubleBracket]3\[RightDoubleBracket]\))
                      \)\ r\^#1\[LeftDoubleBracket]2\[RightDoubleBracket]\ #1
                      \[LeftDoubleBracket]1\[RightDoubleBracket]\)\/\(r
                      \^\(2\ #1\[LeftDoubleBracket]2\[RightDoubleBracket] - 
                        2\ #1\[LeftDoubleBracket]3\[RightDoubleBracket]\) - 
                    s\^\(2\ #1\[LeftDoubleBracket]2\[RightDoubleBracket] - 
                        2\ #1\[LeftDoubleBracket]3\[RightDoubleBracket]\)\)]&)
            \)/@Flatten[
            \((append[
                    HarmonicDecomposition[#1\[LeftDoubleBracket]1
                        \[RightDoubleBracket], z], #1\[LeftDoubleBracket]2
                      \[RightDoubleBracket]]&)\)/@
              polynomialDecomposition[f, z], 1]\) /; Dimension[z] == 2\)

\!\(annular[f_, r_, s_, z_] := 
    Plus@@\(\((
          \(\((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^\(2 - 
                        Dimension[z] + 
                        2\ #1\[LeftDoubleBracket]2\[RightDoubleBracket] - 
                        2\ #1\[LeftDoubleBracket]3\[RightDoubleBracket]\) - 
                    s\^\(2 - Dimension[z] + 
                        2\ #1\[LeftDoubleBracket]2\[RightDoubleBracket] - 
                        2\ #1\[LeftDoubleBracket]3\[RightDoubleBracket]\))\)\ 
                r\^#1\[LeftDoubleBracket]2\[RightDoubleBracket]\ #1
                  \[LeftDoubleBracket]1\[RightDoubleBracket]\)\/\(r
                  \^\(2 - Dimension[z] + 
                    2\ #1\[LeftDoubleBracket]2\[RightDoubleBracket] - 
                    2\ #1\[LeftDoubleBracket]3\[RightDoubleBracket]\) - 
                s\^\(2 - Dimension[z] + 
                    2\ #1\[LeftDoubleBracket]2\[RightDoubleBracket] - 
                    2\ #1\[LeftDoubleBracket]3\[RightDoubleBracket]\)\)&)\)/@
        Flatten[\((
              append[HarmonicDecomposition[#1\[LeftDoubleBracket]1
                      \[RightDoubleBracket], z], #1\[LeftDoubleBracket]2
                    \[RightDoubleBracket]]&)\)/@
            polynomialDecomposition[f, z], 1]\)\)

append[b_,d_]:=(Append[#1,d]&)/@b

AnnularDirichlet[f_,g_,r_,s_,z_Symbol]:=annular[f,r,s,z]+annular[g,s,r,z]

\!\(AnnularDirichlet[f_, g_, r0_, r1_, v_List] := 
    \((unmakeVector[v]; 
      Module[{x}, setDimension[x, Length[v]]; 
        \(AnnularDirichlet[
              f /. Thread[v \[Rule] Table[x\_j, \ {j, \ Length[v]}]], 
              g /. Thread[v \[Rule] Table[x\_j, \ {j, \ Length[v]}]], r0, r1, 
              x] /. Thread[Table[x\_j, \ {j, \ Length[v]}] \[Rule] v]\) /. 
          \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] \[Rule] 
            \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]])\)\)

AnnularDirichlet[f0_,f1_,g_,r0_,r1_,z_Symbol]:=
  With[{t=AntiLaplacian[g,z]},
    t+AnnularDirichlet[
        f0-t/.\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\[Rule]r0,
        f1-t/.\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\[Rule]r1,
        r0,r1,z]]

\!\(AnnularDirichlet[f0_, f1_, g_, r0_, r1_, v_List] := 
    \((unmakeVector[v]; 
      Module[{x}, setDimension[x, Length[v]]; 
        \(AnnularDirichlet[
              f0 /. Thread[v \[Rule] Table[x\_j, \ {j, \ Length[v]}]], 
              f1 /. Thread[v \[Rule] Table[x\_j, \ {j, \ Length[v]}]], r0, 
              r1, g /. Thread[v \[Rule] Table[x\_j, \ {j, \ Length[v]}]], 
              x] /. Thread[Table[x\_j, \ {j, \ Length[v]}] \[Rule] v]\) /. 
          \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] \[Rule] 
            \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]])\)\)



Neumann[f_,v_List]:=explicit[Neumann,f,v]

Neumann[f_,g_,v_List]:=explicit[Neumann,f,g,v]

Neumann[0,x_Symbol]:=0

\!\(Neumann[f_, x_Symbol] := 
    Plus@@\(\((
          If[#1\[LeftDoubleBracket]2\[RightDoubleBracket] === 0, 
              0, #1\[LeftDoubleBracket]1\[RightDoubleBracket]\/#1
                  \[LeftDoubleBracket]2\[RightDoubleBracket]]&)\)/@
        polynomialDecomposition[Dirichlet[f, x], x]\)\)

\!\(Neumann[f_, g_, x_Symbol] := 
    With[{v = AntiLaplacian[g, x]}, 
      v + Neumann[f - NormalD[v, x], x] - 
        Limit[v /. {x\__ \[Rule] 0, 
              \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] \[Rule] 
                t}, t \[Rule] 0, Direction \[Rule] \(-1\)]]\)



ExteriorNeumann[f_,v_List]:=explicit[ExteriorNeumann,f,v]

ExteriorNeumann[f_,g_,v_List]:=explicit[ExteriorNeumann,f,g,v]

ExteriorNeumann[0,x_Symbol]:=0

\!\(ExteriorNeumann[f_, x_Symbol] := 
    Plus@@\(\((
            \(\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^\(2 - 
                      Dimension[x] - 
                      2\ #1\[LeftDoubleBracket]2\[RightDoubleBracket]\)\ #1
                    \[LeftDoubleBracket]1\[RightDoubleBracket]\)\/\(Dimension[
                    x] + #1\[LeftDoubleBracket]2\[RightDoubleBracket] - 2\)&)
            \)/@polynomialDecomposition[Dirichlet[f, x], x]\) /; 
      Dimension[x] =!= 2\)

\!\(ExteriorNeumann[f_, x_Symbol] := 
    Plus@@\(\((
            \[LeftDoubleBracketingBar]x
                    \[RightDoubleBracketingBar]\^\(\(-2\)\ #1
                      \[LeftDoubleBracket]2\[RightDoubleBracket]\)\ #1
                  \[LeftDoubleBracket]1\[RightDoubleBracket]\ 
                If[#1\[LeftDoubleBracket]2\[RightDoubleBracket] === 0, 0, 
                  1\/#1\[LeftDoubleBracket]2\[RightDoubleBracket]]&)\)/@
          polynomialDecomposition[Dirichlet[f, x], x]\) /; 
      Dimension[x] === 2\)

ExteriorNeumann[f_,g_,x_Symbol]:=
  With[{v=AntiLaplacian[g,x,Singularity\[Rule]0]},
    v+ExteriorNeumann[f+NormalD[v,x],x]]



BiDirichlet[f_,v_List]:=explicit[BiDirichlet,f,v]

BiDirichlet[f_,g_,v_List]:=explicit[BiDirichlet,f,g,v]

\!\(BiDirichlet[f_, z_Symbol] := 
    Together[With[{p = Dirichlet[f, z]}, 
        p + 1\/2\ 
            \((1 - \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2)
              \)\ z . Grad[p, z]]]\)

BiDirichlet[f_,g_,z_Symbol]:=normalIntegral[Dirichlet[g,z],z]+BiDirichlet[f,z]

\!\(normalIntegral[p_, z_] := 
    1\/2\ \((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2 - 1)
        \)\ p\)





ZonalHarmonic[0,x_,z_]:=1

\!\(ZonalHarmonic[m_Integer, x_Symbol, z_Symbol] := 
    \((makeVector[x]; makeVector[z]; 
      Plus@@\(\((
            \(1\/\(\(#1!\)\ \(\((m - 2\ #1)\)!\)\ 2\^#1\)\) 
                \((\((\(-1\))\)\^#1\ \((Dimension[x + z] + 2\ m - 2)\)\ 
                    Times@@Range[Dimension[x + z], 
                        Dimension[x + z] + 2\ m - 2\ #1 - 4, 2]\ 
                    \((x . z)\)\^\(m - 2\ #1\)\ 
                    \[LeftDoubleBracketingBar]x
                        \[RightDoubleBracketingBar]\^\(2\ #1\)\ 
                    \[LeftDoubleBracketingBar]z
                        \[RightDoubleBracketingBar]\^\(2\ #1\))\)&)\)/@
          Range[0, m\/2]\))\)\)

\!\(ZonalHarmonic[m_Integer, x_List, z_List] := 
    \((unmakeVector[x]; unmakeVector[z]; 
      Expand[Plus@@
          \(\((\((\((\(-1\))\)\^#1\ \((Length[x] + 2\ m - 2)\)\ 
                      Times@@Range[Length[x], Length[x] + 2\ m - 2\ #1 - 4, 2]
                          \ \((x . z)\)\^\(m - 2\ #1\)\ 
                      \[LeftDoubleBracketingBar]x
                          \[RightDoubleBracketingBar]\^\(2\ #1\)\ 
                      \[LeftDoubleBracketingBar]z
                          \[RightDoubleBracketingBar]\^\(2\ #1\))\)/
                  \((\(#1!\)\ \(\((m - 2\ #1)\)!\)\ 2\^#1)\)&)\)/@
            Range[0, m\/2]\)])\)\)



DimensionH[0,2]=1;

DimensionH[m_,n_]:=Binomial[n+m-1,n-1]-Binomial[n+m-3,n-1]



Options[BasisH]={Orthonormal\[Rule]None};

\!\(BasisH[0, x_, options___?OptionQ] := 
    With[{orthonormal = \(Orthonormal /. {options}\) /. Options[BasisH]}, 
      If[orthonormal === None, {1}, 
        Which[orthonormal === Sphere, {1}, 
          orthonormal === Ball, {1\/\@Volume[Dimension[x]]}]]]\)

\!\(BasisH[m_Integer, x_Symbol, options___?OptionQ] := 
    postTogether[
        With[{orthonormal = \(Orthonormal /. {options}\) /. Options[BasisH]}, 
          Module[{xx}, 
            With[{ttt = 
                  Expand[Expand[
                        \((Kelvin[#1, xx]&)\)/@
                          \(\((Partial@@
                                  Prepend[Transpose[{{xx\_1, xx\_2}, #1}], 
                                    Log[\[LeftDoubleBracketingBar]xx
                                        \[RightDoubleBracketingBar]]]&)\)/@
                            Join[\((Prepend[#1, 0]&)\)/@multiIndices[m, 1], 
                              \((Prepend[#1, 1]&)\)/@multiIndices[m - 1, 1]]
                              \)]] /. xx \[Rule] x}, 
              If[orthonormal === None, 
                noLeadingMinus/@
                  \(Transpose[FactorTermsList/@ttt]\)\[LeftDoubleBracket]2
                    \[RightDoubleBracket], 
                orthonormalize[ttt, 
                  \(Which[orthonormal === Sphere, IntegrateSphere, 
                        orthonormal === Ball, IntegrateBall]\)[#1\ #2, 
                      x]&]]]]]] /; Dimension[x] == 2\)

\!\(BasisH[m_Integer, x_Symbol, options___?OptionQ] := 
    postTogether[
        With[{orthonormal = \(Orthonormal /. {options}\) /. Options[BasisH]}, 
          Module[{xx}, 
            With[{ttt = 
                  Expand[Expand[
                        \((Kelvin[#1, xx]&)\)/@
                          \(\((Partial@@
                                  Prepend[
                                    Transpose[{
                                        Table[\ xx\_j, 
                                        \ {j, \ Dimension[x]}], #1}], 
                                    \[LeftDoubleBracketingBar]xx
                                        \[RightDoubleBracketingBar]^
                                      \((2 - Dimension[x])\)]&)\)/@
                            Join[\((Prepend[#1, 0]&)\)/@
                                multiIndices[m, Dimension[x] - 1], 
                              \((Prepend[#1, 1]&)\)/@
                                multiIndices[m - 1, Dimension[x] - 1]]\)]] /. 
                    xx \[Rule] x}, 
              If[orthonormal === None, 
                noLeadingMinus/@
                  \(Transpose[FactorTermsList/@ttt]\)\[LeftDoubleBracket]2
                    \[RightDoubleBracket], 
                orthonormalize[ttt, 
                  \(Which[orthonormal === Sphere, IntegrateSphere, 
                        orthonormal === Ball, IntegrateBall]\)[#1\ #2, 
                      x]&]]]]]] /; Head[Dimension[x]] == Integer\)

\!\(BasisH[m_Integer, v_List, options___?OptionQ] := 
    \((unmakeVector[v]; 
      Module[{z}, setDimension[z, Length[v]]; 
        Expand[BasisH[m, z, options] /. 
              Thread[Table[\ z\_j, \ {j, Length[v]}] \[Rule] v]] /. 
          \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar] \[Rule] 
            \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]])\)\)

noLeadingMinus[x_]:=
  If[Head[x]===Plus,
    If[Head[x\[LeftDoubleBracket]1\[RightDoubleBracket]]===Times,
      If[x\[LeftDoubleBracket]1,1\[RightDoubleBracket]<0,Expand[-x],x,x],
      If[x\[LeftDoubleBracket]1\[RightDoubleBracket]<0,Expand[-x],x,x]],x]



\!\(orthonormalize[b_List, innerProduct_] := \ 
    Module[{e, j, t}, 
      Do[e[j] = 
          b\[LeftDoubleBracket]j\[RightDoubleBracket] - 
            Plus@@\(\((
                  innerProduct[b\[LeftDoubleBracket]j\[RightDoubleBracket], 
                        e[#1]]\ e[#1]/t[#1]&)\)/@Range[j - 1]\); 
        t[j]\  = \ innerProduct[e[j], e[j]], {j, 1, Length[b]}]; 
      \(Together[e[#1]\/\@t[#1]]&\)/@Range[Length[b]]]\)





Inversion[\[Infinity]]:=0

\!\(Inversion[x_] := 
    \((If[Head[x] == Symbol, makeVector[x]]; 
      If[\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] === 0, 
        \[Infinity], 
        1/\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2\ x])\)\)

\!\(Inversion[x_, y_] := 
    \((If[Head[x] == Symbol, makeVector[x]]; 
      If[\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] === 0 && 
          y === 0, \[Infinity], \ 
        1/\((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2\  + \ 
                y\^2)\)\ \((x\  - \ y\ S)\)])\)\)



Kelvin[HoldPattern[+u__],x_]:=Plus@@Thread[Kelvin[{u},x]]

\!\(Kelvin[u_, x_List] := 
    \((unmakeVector[x]; 
      \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^\(2 - 
            Length[x]\)\ \((u /. Thread[x \[Rule] Inversion[x]])\))\)\)

\!\(Kelvin[u_, x_Symbol] := 
    \((makeVector[x]; 
      \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^\(2 - 
            Dimension[x]\)\ 
        \((u /. {x \[Rule] Inversion[x], 
              x\_j_ \[RuleDelayed] 
                x\_j\/\[LeftDoubleBracketingBar]x
                      \[RightDoubleBracketingBar]\^2})\))\)\)



Phi[{x_?vector,y_}]:=Phi[x,y]

Phi[z_List]:=(unmakeVector[z];
    If[Last[z]===-1&&
        \[LeftDoubleBracketingBar]Drop[z,-1]\[RightDoubleBracketingBar]===0,
      \[Infinity],
      2 Inversion[z+Append[Table[0,{Length[z]-1}],1]]-
        Append[Table[0,{Length[z]-1}],1]])

Phi[z_]:=(
    If[!(z===0||z===\[Infinity]||MatchQ[z,_. S]),makeVector[z];
      If[Dimension[S]=!= Dimension[z],SetDimension[S,Dimension[z]]]];
    2 Inversion[z-S]+S)

Phi[0,-1]=\[Infinity];

\!\(Phi[x_, y_] := \ 
    \((makeVector[x]; 
      unmakeVector[y]; {
        \(2\ x\)\/\(\[LeftDoubleBracketingBar]x
                \[RightDoubleBracketingBar]\^2 + \((y + 1)\)\^2\), 
        \(2\ \((y + 1)
                \)\)\/\(\[LeftDoubleBracketingBar]x
                  \[RightDoubleBracketingBar]\^2 + \((y + 1)\)\^2\) - 1})\)\)



\[LeftDoubleBracketingBar]S\[RightDoubleBracketingBar]^=1;

\!\(S . z_Symbol := \(-z\_\(Dimension[z]\)\) /; matrix[z] =!= True\)

\!\(a_Symbol . S := \(-a\_\(Dimension[a]\)\) /; matrix[z] =!= True\)

makeVector[S];

\!\(\(N[S\_j_Real]\  := \ S\_\(Round[j]\) /; \ Round[j]\  == \ j; \)\)

\!\(S\_j_ := \(-1\) /; j == Dimension[S]\)

\!\(S\_j_ := 0 /; j < Dimension[S]\)



KelvinM[HoldPattern[+u__],x_]:=Plus@@Thread[KelvinM[{u},x]]

\!\(KelvinM[u_, x_List] := 
    \((unmakeVector[x]; 
      2\^\(1\/2\ \((Length[x] - 2)\)\)\ 
        \((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2 + 
              2\ Last[x] + 1)\)\^\(2 - Length[x]\)\ 
        \((u /. Thread[x \[Rule] Phi[x]])\))\)\)

\!\(KelvinM[u_, z_Symbol] := 
    \((makeVector[z]; 
      2\^\(1\/2\ \((Dimension[z] - 2)\)\)\ 
        \[LeftDoubleBracketingBar]z - S\[RightDoubleBracketingBar]\^\(2 - 
            Dimension[z]\)\ 
        \((u /. {z \[Rule] Phi[z], z\_j_ \[RuleDelayed] Phi[z] . Delta[j]})\))
      \)\)

\!\(KelvinM[u_, x_Symbol, y_Symbol] := 
    \((makeVector[x]; unmakeVector[y]; 
      2\^\(1\/2\ \((Dimension[x] - 1)\)\)\ 
        \((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2 + 
              \((y + 1)\)\^2)\)\^\(1\/2\ \((1 - Dimension[x])\)\)\ 
        \((u /. {
              x \[Rule] 
                \(2\ x\)\/\(\[LeftDoubleBracketingBar]x
                        \[RightDoubleBracketingBar]\^2 + \((y + 1)\)\^2\), 
              y \[Rule] 
                \(2\ \((y + 1)
                        \)\)\/\(\[LeftDoubleBracketingBar]x
                          \[RightDoubleBracketingBar]\^2 + \((y + 1)\)\^2\) - 
                  1, x\_j_ \[RuleDelayed] 
                \(2\ x\_j
                      \)\/\(\[LeftDoubleBracketingBar]x
                        \[RightDoubleBracketingBar]\^2 + \((y + 1)\)\^2\)})\))
      \)\)





\!\(PoissonKernel[x_, z_] := 
    \((makeVector[x]; makeVector[z]; 
      \(1 - \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2\ 
            \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2\)\/\((
            1 - 2\ x . z + 
              \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2\ 
                \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2)
            \)\^\(1\/2\ Dimension[x + z]\))\)\)



\!\(PoissonKernelH[z_, w_] := 
    \((makeVector[z]; makeVector[w]; 
      \((2\ \((z . Delta[Dimension[z]] + w . Delta[Dimension[z]])\))\)/
        \((\((Dimension[z]\ Volume[Dimension[z]])\)\ 
            \((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2 + 
                  \[LeftDoubleBracketingBar]w\[RightDoubleBracketingBar]\^2 - 
                  2\ z . w + 
                  4\ z . Delta[Dimension[z]]\ w . Delta[Dimension[z]])\)^
              \((Dimension[z]\/2)\))\))\)\)

\!\(PoissonKernelH[x_, y_, t_, u_] := 
    \((makeVector[x]; makeVector[t]; unmakeVector[y]; unmakeVector[u]; 
      \((2\ \((y + u)\))\)/
        \((\((\((Dimension[x] + 1)\)\ Volume[Dimension[x] + 1])\)\ 
            \((\[LeftDoubleBracketingBar]x - t
                        \[RightDoubleBracketingBar]\^2 + \((y + u)\)\^2)
                \)\^\(1\/2\ \((Dimension[x] + 1)\)\))\))\)\)



\!\(BergmanKernel[x_, y_] := 
    \((makeVector[x]; makeVector[y]; 
      With[{n = Dimension[x + y]}, 
        \(n + \((8\ x . y - 2\ n - 4)\)\ 
              \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2\ 
              \[LeftDoubleBracketingBar]y\[RightDoubleBracketingBar]\^2 + 
            \((n - 4)\)\ 
              \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^4\ 
              \[LeftDoubleBracketingBar]y\[RightDoubleBracketingBar]\^4\)\/\(n
            \ Volume[n]\ 
            \((1 - 2\ x . y + 
                  \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\^2\ 
                    \[LeftDoubleBracketingBar]y\[RightDoubleBracketingBar]\^2)
                \)\^\(1 + n\/2\)\)])\)\)



\!\(BergmanKernelH[z_, w_] := 
    \((makeVector[z]; makeVector[w]; 
      With[{n = Dimension[z + w]}, 
        \((4\ \((
                n\ \((\((z . Delta[n])\)\^2 + \((w . Delta[n])\)\^2)\) - 
                  \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2 - 
                  \[LeftDoubleBracketingBar]w\[RightDoubleBracketingBar]\^2 + 
                  2\ z . w + \((2\ n - 4)\)\ z . Delta[n]\ w . Delta[n])\))\)/
          \((n\ Volume[n]\ 
              \((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]\^2 + 
                    \[LeftDoubleBracketingBar]w
                        \[RightDoubleBracketingBar]\^2 - 2\ z . w + 
                    4\ z . Delta[n]\ w . Delta[n])\)\^\(1 + n\/2\))\)])\)\)

\!\(BergmanKernelH[x_, y_, t_, u_] := 
    \((makeVector[x]; makeVector[t]; unmakeVector[t]; unmakeVector[u]; 
      With[{n = Dimension[x] + 1}, 
        \(4\ \((\((n - 1)\)\ \((y + u)\)\^2 - 
                \[LeftDoubleBracketingBar]x - t\[RightDoubleBracketingBar]\^2)
              \)\)\/\(n\ Volume[n]\ 
            \((\[LeftDoubleBracketingBar]x - t
                        \[RightDoubleBracketingBar]\^2 + \((y + u)\)\^2)
                \)\^\(1 + n\/2\)\)])\)\)





\!\(HarmonicConjugate[u_, x_Symbol, y_Symbol] := 
    \((unmakeVector[x]; unmakeVector[y]; 
      \[Integral]\_0\%y\[PartialD]\_x u \[DifferentialD]y - 
        \[Integral]\_0\%x\((\[PartialD]\_y u /. y \[Rule] 0)\) 
            \[DifferentialD]x)\)\)



BergmanProjection[0,x_]:=0

\!\(BergmanProjection[u_, x_Symbol] := 
    \((makeVector[x]; 
      Plus@@\(\((Plus@@#1&)\)/@
          \(\((\(Function[j, 
                      \(\((Dimension[x] + 
                                2\ #1\[LeftDoubleBracket]2
                                    \[RightDoubleBracket])\)\ #1
                              \[LeftDoubleBracket]1
                              \[RightDoubleBracket]\)\/\(#1
                              \[LeftDoubleBracket]2\[RightDoubleBracket] + 
                            j + Dimension[x]\)&]\)[#1\[LeftDoubleBracket]2
                      \[RightDoubleBracket]]/@#1\[LeftDoubleBracket]1
                    \[RightDoubleBracket]&)\)/@
            \(\(({polynomialDecomposition[#1\[LeftDoubleBracket]1
                        \[RightDoubleBracket], x], #1\[LeftDoubleBracket]2
                      \[RightDoubleBracket]}&)\)/@
              \(\(({Dirichlet[#1\[LeftDoubleBracket]1\[RightDoubleBracket], 
                        x], #1\[LeftDoubleBracket]2\[RightDoubleBracket]}&)
                  \)/@polynomialDecomposition[u, x]\)\)\)\))\)\)

BergmanProjection[u_,v_List]:=explicit[BergmanProjection,u,v]



Schwarz[0]=0;

\!\(Schwarz[x_] := 
    \((If[Head[x] == List, unmakeVector[x]]; 
        \(Module[{t, u}, 
              With[{s = 
                    \[Integral]\_0
                        \%\(\[Pi]\/2
                            \)\(Cos[u]\^\(Dimension[x] - 2\)\/\((
                                1 + 1\/t\^2 - \(2\ Sin[u]\)\/t)
                                \)\^\(Dimension[x]\/2\)\) \[DifferentialD]u}, 
                Together[
                  \((\((1 - 1\/t\^2)\)\ \((Dimension[x] - 1)\)\ 
                          Volume[Dimension[x] - 1]\ 
                          \((s - \((s /. t \[Rule] \(-t\))\))\))\)/
                      \((Dimension[x]\ Volume[Dimension[x]])\) /. 
                    t \[Rule] 
                      1\/\[LeftDoubleBracketingBar]x
                          \[RightDoubleBracketingBar]]]] /. 
            HoldPattern[\(+t_\)] \[RuleDelayed] 
              Factor[Select[\(+t\), FreeQ[#1, ArcTan]&]] + 
                Factor[Select[\(+t\), 
                    \(\[InvisibleSpace]\(! \((FreeQ[#1, ArcTan])\)\)\)&]]
                  \) /. \((\(-1\) + t_)\)\ \((1 + t_)\) \[Rule] 
            \(-\((1 - t\^2)\)\))\) /; EvenQ[Dimension[x]]\)

\!\(Schwarz[x_] := 
    \((If[Head[x] == List, unmakeVector[x]]; 
        Apart[\(Module[{t, u}, 
                With[{s = 
                      \[Integral]\_0
                          \%\(\[Pi]\/2
                              \)\(Cos[u]\^\(Dimension[x] - 2\)\/\((
                                  1 + t\^2 - 2\ Sin[u]\ t)
                                  \)\^\(Dimension[x]\/2\)\) 
                          \[DifferentialD]u}, 
                  Together[
                    PowerExpand[
                      \(\((\((1 - t\^2)\)\ \((Dimension[x] - 1)\)\ 
                                Volume[Dimension[x] - 1]\ 
                                \((s - \((s /. t \[Rule] \(-t\))\))\))\)/
                            \((Dimension[x]\ Volume[Dimension[x]])\) /. 
                          t \[Rule] 
                            \[LeftDoubleBracketingBar]x
                              \[RightDoubleBracketingBar]\) /. {
                          \((\(-1\) + 
                                  \[LeftDoubleBracketingBar]x
                                    \[RightDoubleBracketingBar])\)\^p_?EvenQ 
                            \[RuleDelayed] 
                            \((1 - 
                                  \[LeftDoubleBracketingBar]x
                                    \[RightDoubleBracketingBar])\)\^p, 
                          \((\(-1\) - 
                                  \[LeftDoubleBracketingBar]x
                                    \[RightDoubleBracketingBar])\)\^p_?EvenQ 
                            \[RuleDelayed] 
                            \((1 + 
                                  \[LeftDoubleBracketingBar]x
                                    \[RightDoubleBracketingBar])\)\^p}]]]] /. 
              HoldPattern[\(+t_\)] \[RuleDelayed] 
                Factor[Select[\(+t\), FreeQ[#1, \@_]&]] + 
                  Factor[Select[\(+t\), 
                      \(\[InvisibleSpace]\(! \((FreeQ[#1, \@_])\)\)\)&]]\) /. 
            \((\(-1\) + t_)\)\ \((1 + t_)\) \[Rule] \(-\((1 - t\^2)\)\)])
        \) /; OddQ[Dimension[x]]\)



IntegerQ[m_?IntegerQ+n_?IntegerQ]=True;

IntegerQ[m_?IntegerQ n_?IntegerQ]=True;









\!\(explicit[function_, f_, v_] := 
    \((unmakeVector[v]; 
      Module[{z}, setDimension[z, Length[v]]; 
        function[f /. Thread[v \[Rule] Table[\ z\_j, \ {j, Length[v]}]], 
            z] /. Prepend[Thread[Table[\ z\_j, \ {j, Length[v]}] \[Rule] v], 
            \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar] \[Rule] 
              \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]]])\)\)

\!\(explicit[function_, f_, g_, v_List] := 
    \((unmakeVector[v]; 
      Module[{z}, setDimension[z, Length[v]]; 
        function[f /. Thread[v \[Rule] Table[\ z\_j, \ {j, Length[v]}]], 
            g /. Thread[v \[Rule] Table[\ z\_j, \ {j, Length[v]}]], z] /. 
          Prepend[Thread[Table[\ z\_j, \ {j, Length[v]}] \[Rule] v], 
            \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar] \[Rule] 
              \[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]]])\)\)



Protect[Release[protectedWords]];

End[]

EndPackage[]

Print["* You can now use the functions in this package."]