(* ::Package:: *)

BeginPackage["GridZeta`",{"ComplexUtilities`"}]


from::inverse = "Inverse functions are being used, and values will be missed.";

(* Conversions only involving hyperbolic space *)
tfrom\[Tau]::usage= "Convert \[Tau] \[Element] H to t=\!\(\*SqrtBox[\(q\)]\).";
qfrom\[Tau]::usage= "Convert \[Tau] \[Element] H to q.";
\[Tau]fromq::usage= "Convert q to (one value of) \[Tau] \[Element] H. (Warning: logarithm)";
\[Tau]fromt::usage= "Convert t=\!\(\*SqrtBox[\(q\)]\) to (one value of) \[Tau] \[Element] H. (Warning: logarithm)";
qfromt::usage = "Convert t to q=\!\(\*SuperscriptBox[\(t\), \(2\)]\).";
tfromq::usage = "Convert q to t=\!\(\*SqrtBox[\(q\)]\). (Warning: square root)";

(* Conversions involving elliptic parameters m and k *)
mfromk::usage = "Convert k to m=\!\(\*SuperscriptBox[\(k\), \(2\)]\).";
kfromm::usage = "Convert m=\!\(\*SuperscriptBox[\(k\), \(2\)]\) to k.  (Warning: square root)";

mfromq::usage = "Convert the nome q to m=\!\(\*SuperscriptBox[\(k\), \(2\)]\).  Uses InverseEllipticNomeQ.";
mfromt::usage = "Convert t=\!\(\*SqrtBox[\(q\)]\) to m=\!\(\*SuperscriptBox[\(k\), \(2\)]\).  Uses InverseEllipticNomeQ.";
mfrom\[Tau]::usage = "Convert m=\!\(\*SuperscriptBox[\(k\), \(2\)]\) to \[Tau] \[Element] H.  Uses InverseEllipticNomeQ.";

kfromq::usage = "Convert q to k. (Warning: square root, use kfromt if possible.)";
kfromt::usage = "Convert t=\!\(\*SqrtBox[\(q\)]\) to k. Uses QPochhammer.";
kfrom\[Tau]::usage = "Convert \[Tau] \[Element] H to k.  Uses QPochhammer.";

qfromm::usage = "Convert m=\!\(\*SuperscriptBox[\(k\), \(2\)]\) to q.  Uses EllipticNomeQ.";
\[Tau]fromm::usage = "Convert m=\!\(\*SuperscriptBox[\(k\), \(2\)]\) to \[Tau] \[Element] H.  Uses EllipticK.";
tfromm::usage = "Convert m=\!\(\*SuperscriptBox[\(k\), \(2\)]\) to t=\!\(\*SqrtBox[\(q\)]\).  Uses EllipticK.";
\[Tau]fromk::usage = "Convert k to \[Tau] \[Element] H. Built piecewise from \[Tau]fromm.";
tfromk::usage = "Convert k to t=\!\(\*SqrtBox[\(q\)]\). Uses \[Tau]fromk.";
qfromk::usage = "Convert k to q.  Uses qfromm[\!\(\*SuperscriptBox[\(k\), \(2\)]\)].";

(* Conversions Involving u *)
kfromu::usage = "Convert u to k";
mfromu::usage = "Convert u to m=\!\(\*SuperscriptBox[\(k\), \(2\)]\)";
ufromk::usage = "Convert k to a list of two u values.";
tfromu::usage = "Convert u to t=\!\(\*SqrtBox[\(q\)]\)";
ufromt::usage = "Convert t=\!\(\*SqrtBox[\(q\)]\)to a list of two u values.";


(* Modular group and fundamental domains *)
fundDomain\[Tau]m::usage = "Graphics directives describing the edge of the fundamental domain for m in the \[Tau] plane.";
fundDomain\[Tau]k::usage = "Graphics directives describing the edge of the fundamental domain for k in the \[Tau] plane.";
fundDomaintk::usage = "Graphics directives describing the edge of the fundamental domain for k in the t disk.";
fundDomainqm::usage = "Graphics directives describing the edge of the fundamental domain for m in the q disk.";
fundtk::usage = "fundtk[g] plots the fundamental domain for k in the t disk, after applying the mobius transformation given by 2x2 matrix g.";
fund\[Tau]::usage = "fund\[Tau][x,y] runs over the fundamental domain in the \[Tau] plane.\
 Using (x,y)\[Element](-1,1]\[Times][0,1) gives the fundamental domain for m,\
 and using (x,y)\[Element](-2,2]\[Times][0,1) gives the fundamental domain for k.";

generatorU::usage = "The 2x2 matrix (\[NoBreak]\*GridBox[{
{1, 2},
{0, 1}
}]\[NoBreak]), a generator for \[CapitalGamma](2).";
generatorV::usage = "The 2x2 matrix (\[NoBreak]\*GridBox[{
{1, 0},
{2, 1}
}]\[NoBreak]), a generator for \[CapitalGamma](2).";


(* Elliptic integrals *)
Kfromq::usage = "The elliptic integral K as a function of q.  Equal to EllipticK[mfromq[q]] for q \[Element] (-1,1), but analytic on the disk.";


(* Zeta Function Related *)

gridZetaf::usage = "This is f(t) = (1-\!\(\*SuperscriptBox[SubscriptBox[\(\[Theta]\), \(3\)], \(2\)]\)\!\(\*SuperscriptBox[SubscriptBox[\(\[Theta]\), \(4\)], \(4\)]\))/t from the zetagrid paper.";
gridZetaF::usage = "This is F(t), the primitive for zetaf[t] = f(t).";
gridZetat::usage = "This is the zeta function for the grid, as a function of t.  It is a list of the two values for each t.";
gridZeta::usage = "This is the zeta function for the grid, as a function of u.";
trLogDeltaCoefficient::usage = "trLogDeltaCoefficient[n] is the coefficient of \!\(\*SuperscriptBox[\(u\), \(2  n\)]\) in the power series for \!\(\*SubscriptBox[\(Tr\), \(\[Pi]\)]\)\!\(\*SubscriptBox[\(Log\[CapitalDelta]\), \(u\)]\), as a function of u.  Computed using the combinatorical formula.";


baduSet3D::usage = "baduSet3D[edgeform,height] gives the set on which u has issues, at given height (default is 1).";
baduSet::usage = "baduSet[edgeform] gives the set on which u has issues.";


Begin["`Private`"]


Kfromq[q_]:=(Pi/2)*EllipticTheta[3,0,q]^2


(* Conversions only involving hyperbolic space *)
tfrom\[Tau][\[Tau]_]:=Exp[I Pi \[Tau]/2]
qfrom\[Tau][\[Tau]_]:=Exp[I Pi \[Tau]]
\[Tau]fromq[q_]:=(Message[from::inverse]; Log[q]/(I Pi))
\[Tau]fromt[t_]:=(Message[from::inverse]; 2 Log[t]/(I Pi))
qfromt[t_]:=t^2
tfromq[q_]:=(Message[from::inverse]; Sqrt[q])

(* Conversions involving elliptic parameters m and k *)
mfromk[k_]:=k^2
kfromm[m_]:=(Message[from::inverse]; Sqrt[m])

mfromq[q_]:=InverseEllipticNomeQ[q]
mfromt[t_]:=InverseEllipticNomeQ[t^2]
mfrom\[Tau][\[Tau]_]:=InverseEllipticNomeQ[qfrom\[Tau][\[Tau]]]

kfromq[q_]:=kfromm[mfromq[q]]
kfromt[t_]:=(4t)/(QPochhammer[-t^2,t^4]QPochhammer[t^4,t^8])^4
kfrom\[Tau][\[Tau]_]:=kfromt[tfrom\[Tau][\[Tau]]]

qfromm[m_]:=EllipticNomeQ[m]
\[Tau]fromm[m_]:=(I EllipticK[1-m])/EllipticK[m]
tfromm[m_]:=tfrom\[Tau][\[Tau]fromm[m]]

qfromk[k_]:=qfromm[k^2]
\[Tau]fromk[k_]:= If[Re[k]>=0,\[Tau]fromm[k^2],
				If[Im[k]>=0,\[Tau]fromm[k^2]+2,
							\[Tau]fromm[k^2]-2]]
tfromk[k_]:=tfrom\[Tau][\[Tau]fromk[k]]

(* Conversions Involving u *)
kfromu[u_]:=(4u)/(1+3u^2)
mfromu[u_]:=kfromu[u]^2
ufromk[k_]:={(2-Sqrt[4-3 k^2])/(3 k),(2+Sqrt[4-3 k^2])/(3 k)}
tfromu[u_]:=tfromk[kfromu[u]]
ufromt[t_]:=ufromk[kfromt[t]]


fundDomain\[Tau]m = {Arrow[{{1,0},{1,4}}],Arrow[{{-1,0},{-1,4}}],
				Circle[{1/2,0},1/2,{0,Pi}],Circle[{-1/2,0},1/2,{0,Pi}]};
fundDomain\[Tau]k = {Arrow[{{2,0},{2,4}}],Arrow[{{-2,0},{-2,4}}],
				Circle[{1/2,0},1/2,{0,Pi}],Circle[{-1/2,0},1/2,{0,Pi}],
				Circle[{3/2,0},1/2,{0,Pi}],Circle[{-3/2,0},1/2,{0,Pi}]};
fundDomaintk = (ParametricPlot[
					Evaluate[Table[ComplexSplit[tfrom\[Tau][c/2+Exp[2I x]/2]],{c,{1,3,-3,-1}}]],
					{x,0,\[Pi]/2}]
				//Graphics)[[1]]//Flatten//Cases[#,_Line]&;
fundDomainqm = (ParametricPlot[
					Evaluate[Table[ComplexSplit[qfrom\[Tau][c/2+Exp[2I x]/2]],{c,{1,-1}}]],
					{x,0,\[Pi]/2}]
				//Graphics)[[1]]//Flatten//Cases[#,_Line]&;
fundtk[g_,options___]:=ParametricPlot[Evaluate[
	ComplexSplit/@
		(tfrom\[Tau]/@
			(MobiusTransformation[g,#]&/@
				{-2+I t^3,-3/2+(1/2)Exp[I t],-1/2+(1/2)Exp[I t],
				1/2+(1/2)Exp[I t],3/2+(1/2)Exp[I t],2+I t^3})
		 )
	],{t,0,Pi},Axes->None,options]

(* smoothedge numerically solves the heat equation to get a smooth approximation to the edge of the \[Tau] fundamental domain,
   which consists of half circles and has cusps at each integer. *)
smoothedge=f/.NDSolve[{D[f[x,t],t] == D[f[x,t],x,x],
						f[x,0]==Sqrt[1/4-(Abs[x]-1/2)^2]-1/2,
						f[-1/2,t]==0,
						f[1/2,t]==0},
						f,{x,-1/2,1/2},{t,0,1}][[1]];
(* To produce the fundamental domain for k in the \[Tau] half plane, replicate the smoothed edge periodically and
   add some ad-hoc terms to stretch it out nicely over the domain *)
fund\[Tau][x_,y_]:=x+I * (smoothedge[Mod[x-1/2,1]-1/2,y^4]+1/(2-2y))


(* Some routines for working with the modular group \[CapitalGamma](2) *)
generatorU = {{1,2},{0,1}};
generatorV = {{1,0},{2,1}};


Kfromq[q_]:=(Pi/2)*EllipticTheta[3,0,q]^2


(* This gives the grid zeta function using elliptic theta functions *)
gridZetaf[t_]:=(1-(EllipticTheta[3,0,t^2]^2)(EllipticTheta[4,0,t^2]^4))/t;
gridZetaF[t_]:=Evaluate[Normal[Integrate[Series[gridZetaf[t],{t,0,256}],t]]];
gridZetat[t_]:=t Exp[-gridZetaF[t]]/((#-#^3)&@ufromt[t]);
gridZeta[u_]:=((#*Exp[-gridZetaF[#]])&@tfromu[u])/(u-u^3);

trLogDeltaCoefficient[n_] := -Sum[(-3)^(n-k)/(n+k) Binomial[n+k,2k]Binomial[2k,k]^2,{k,0,n}]


baduSet3D[edge_List,height_:1] := Graphics3D[Join[
		{{FaceForm[None],EdgeForm[edge],
		Polygon[Table[{(1/Sqrt[3])Cos[x],(1/Sqrt[3])Sin[x],height},{x,0,2Pi,Pi/32}]]}},
		edge,{Line[{{-1,0,height},{-1/3,0,height}}],Line[{{1,0,height},{1/3,0,height}}]}
		]]
baduSet3D[edge_,height_:1]:=baduSet3D[{edge},height]

baduSet[edge_List] := Graphics[Join[
		{{FaceForm[None],EdgeForm[edge],
		Polygon[Table[(1/Sqrt[3]){Cos[x],Sin[x]},{x,0,2Pi,Pi/32}]]}},
		edge,{Line[{{-1,0},{-1/3,0}}],Line[{{1,0},{1/3,0}}]}
		]]
baduSet[edge_]:=baduSet[{edge}]


End[] (*`Private`*)


EndPackage[]