BeginPackage["RationalDirichlet`"]

Print["Rational Dirichlet;  Version 1.0,  April 21, 1999."]

Unprotect[DefineReal, RatDirichlet, ExtRatDirichlet, RatHarmonicConjugate,
	Analytic, Laplacian, OnCircle]
Clear[DefineReal, RatDirichlet, ExtRatDirichlet, RatHarmonicConjugate, Analytic,
	Laplacian, OnCircle]

(* Usage definitions *)

RatDirichlet::usage =
	"RatDirichlet[f, x, y] solves the Dirichlet Problem for the polynomial"<>
	"or rational boundary function f given as an expression of x and y."

ExtRatDirichlet::usage =
	"ExtRatDirichlet[f, x, y] solves the Exterior Dirichlet Problem for" <>
	"the polynomial or rational boundary function f given as an expression" <>
	"of x and y."

RatHarmonicConjugate::usage = 
	"RatHarmonicConjugate[u, x, y] gives a harmonic conjugate for the" <>
	"harmonic function u of x and y."

Analytic::usage = 
	"Analytic[g, z] produces a analytic function of the complex" <> 
	"variable z on the unit disc that has the same real part as the" <>
	"expression g on the unit circle."

Laplacian::usage =
	"Laplacian[exp, {x1, x2, ..., xn}] computes the Laplacian of" <>
	"the function exp, with respect to the variables {x1, x2, ..., xn}." <>
	"Laplacian[exp] computes the Laplacian of exp with respect to {x, y}."

OnCircle::usage =
	"OnCircle[f, g, x, y] compares the values of the functions f and g on" <>
	"the disc, where x^2 + y^2 = 1, and returns True if they are equal or" <>
	"False if they are unequal."

(* end: usage definitions *)

Begin["`Private`"]

RatDirichlet::denominator = "Denominator of generated fraction is unfactorable." 
RatDirichlet::continuity = "Your function is not continuous on the circle."
RatHarmonicConjugate::harmonic = "Your function is not harmonic on the circle."

DefineReal[x_] :=
	Block[{p}, p = Unprotect[Conjugate,Re,Im];
	Conjugate[x] = x;
	Re[x] = x;
	Im[x] = 0;
	Protect[Release[p]]; ]

(* public functions *)

RatDirichlet[f_,xvar_,yvar_] :=
	Module[{g, G, F},
	DefineReal[xvar]; DefineReal[yvar];
	g = Simplify[ f/.{xvar -> (z + 1/z)/2, yvar -> (z - 1/z)/(2 I)} ];
	G = Analytic[g, z];
	F = Together[(Block[{z = (xvar + I yvar)}, G + Conjugate[G]]) / 2];
	ExpandDenominator[ Together[ Apart[ F ] ] ] ]

ExtRatDirichlet[f_,xvar_,yvar_] :=
	Module[{g, var, ExtG, F},
	DefineReal[xvar]; DefineReal[yvar];
	g = Simplify[ f/.{xvar -> (z + 1/z)/2, yvar -> (z - 1/z)/(2 I)} ];
	G[var_] = Analytic[g, z]/.{z -> var};
	ExtG = Together[ Apart[ Conjugate[ G[1 / Conjugate[z]] ] ] ];
	F = Together[(Block[{z = (xvar + I yvar)}, ExtG + Conjugate[ExtG]]) / 2];
	ExpandDenominator[ Together[ Apart[ F ] ] ] ]

RatHarmonicConjugate[u_,xvar_,yvar_] :=
	Module[{g, G, v, CL, vZero},
	g = Simplify[ u/.{xvar -> (z + 1/z)/2, yvar -> (z - 1/z)/(2 I)} ];
	G = Analytic[g, z];
	v = Together[ ( Block[{z = (xvar + I yvar)}, G - u] ) / I ];
	v = ExpandDenominator[ Together[ Apart[ v ] ] ];
	CL = Union[ Flatten[ CoefficientList[ Numerator[ v ], {xvar, yvar}] ] ];
	Do[ CL[[i]] = Head[ N[ CL[[i]] ] ], {i,1,Length[CL]} ];
	If[ MemberQ[CL, Complex]==True,
		Message[RatHarmonicConjugate::harmonic];
		Abort[],
		v ];
	vZero = v/.{xvar -> 0, yvar -> 0};
	If[ TrueQ[vZero != 0],
		ExpandDenominator[Together[Apart[ v - vZero ] ] ],
	  	v ] ]

Analytic[exp_,var_] :=
	Module[{g, p, q, Gp, gr, qz, W, Win = {}, Wout, var2, TempOutside,
		TempInside, qz2, qOut, qIn, Win2 = {}, Wout2 = {}, NewInside},

	g = Together[Apart[exp]];
	p = Numerator[g]; q = Denominator[g];
	Gp = Simplify[ PolynomialQuotient[p, q, var] ];
	gr = Simplify[ PolynomialRemainder[p, q, var] / q ];

	qz = Solve[Denominator[gr]==0, var];
	If[ MemberQ[qz/.{{x_Rule}->True,ToRules[y_]->False}, False], 
		Message[RatDirichlet::denominator];
		Abort[] ];

	W = Union[ Table[var/.qz[[i]],{i,1,Length[qz]}] ];
	If[ MemberQ[ Expand[Abs[W]], 1], Message[RatDirichlet::continuity] ];
	Do[Win = Union[Win, If[ N[Abs[(W[[i]])]]<1, {W[[i]]}, {}] ],
		{i,1,Length[W]}];
	Wout = Complement[W, Win];

	If[ Length[Win] <= Length[Wout],
		Inside[var2_] = Together[
			Sum[ Normal[Series[gr, {var, Win[[i]], -1}] ],
			{i,1,Length[Win]}] ]/.{var -> var2};
		TempOutside = Together[ Apart[ Together [ gr - Inside[var] ] ] ];
		qz2 = Solve[Denominator[TempOutside]==0, var];
		qOut = Union[ Table[var/.qz2[[i]],{i,1,Length[qz2]}] ];
		Do[Win2 = Union[Win2, If[ N[Abs[(qOut[[i]])]]<1, {qOut[[i]]}, {}] ],
			{i,1,Length[qOut]}];
		Extra = Product[(var - Win2[[i]]), {i,1,Length[Win2]}];
		Outside[var2_] = (Apart[ Numerator[TempOutside] / Extra] /
			Apart[ Denominator[TempOutside] / Extra])/.{var -> var2};
		NewInside = Conjugate[ Inside[1 / Conjugate[var]] ];
		Simplify[ Together[ Gp + Outside[var] + NewInside ] ],

		Outside[var2_] = Together[
			Sum[ Normal[Series[gr, {var, Wout[[i]], -1}] ],
			{i,1,Length[Wout]}] ]/.{var -> var2};
		TempInside = Together[ Apart[ Together [ gr - Outside[var] ] ] ];
		qz2 = Solve[Denominator[TempInside]==0, var];
		qIn = Union[ Table[var/.qz2[[i]],{i,1,Length[qz2]}] ];
		Do[Wout2 = Union[Wout2, If[ N[Abs[(qIn[[i]])]]<1, {qIn[[i]]}, {}] ],
			{i,1,Length[qIn]}];
		Extra = Product[(var - Wout2[[i]]), {i,1,Length[Wout2]}];
		Inside[var2_] = (Apart[ Numerator[TempInside] / Extra] /
			Apart[ Denominator[TempInside] / Extra])/.{var -> var2};
		NewInside = Conjugate[ Inside[1 / Conjugate[var]] ];
		Simplify[ Together[ Gp + Outside[var] + NewInside ] ] ] ]

Laplacian[exp_,vars_:{Global`x, Global`y}] :=
	Module[{Ans},
	Ans = Sum[D[exp, {vars[[i]], 2}], {i,1,Length[vars]}];
	ExpandNumerator[Together[Ans]] ]

OnCircle[f_,g_,xvar_,yvar_] :=
	Module[{h, n, Convert = {}},
	h = Together[f - g];
	n = Max[ Exponent[Numerator[h], yvar], Exponent[Denominator[h], yvar] ];
	Do[Convert = Join[Convert, { yvar (yvar^2)^i -> yvar (1 - xvar^2)^i},
		{(yvar^2)^i -> (1 - xvar^2)^i}], {i,Floor[n/2],1,-1}];
	TrueQ[ ExpandNumerator[ h/.Convert ] == 0 ] ]

(* End public functions *)

(* Mathematica rule improvements *)

DefineReal[Global`x]
DefineReal[Global`y]

Unprotect[Conjugate, Power]

Conjugate[r_ + s_]:=Conjugate[r] + Conjugate[s]
Conjugate[r_ - s_]:=Conjugate[r] - Conjugate[s]
Conjugate[r_ * s_]:=Conjugate[r] * Conjugate[s]
Conjugate[r_ / s_]:=Conjugate[r] / Conjugate[s]
Conjugate[r_^s_]:=Conjugate[r]^Conjugate[s]

(-1)^(s_Rational):=Cos[s Pi] + I Sin[s Pi]

Protect[ Release[Conjugate, Power] ]

(* End Mathematical rule improvements *)

End[]

Protect[DefineReal, RatDirichlet, RatHarmonicConjugate, Analytic, Laplacian,
	OnCircle]

EndPackage[]