Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 52

Mathematica implementation of cubic theta function

BeginPackage["CubicTheta`"];

Unprotect[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];

ClearAll[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];

Begin["`Private`"];

CubicTheta[u_, v_, q_] /;
(VectorQ[{u, v, q}, NumericQ] && Abs[q] < 1 && Precision[{u, v, q}] < Infinity) :=
Module[{p = q^3, t, a = u + v, b = u - v},
	t = EllipticTheta[3, a, p]*EllipticTheta[3, b, q];
	If[!PossibleZeroQ[p], t += EllipticTheta[2, a, p]*EllipticTheta[2, b, q]*q/(p^(1/4)*q^(1/4))];
	Return[t];
];

CubicThetaA[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
(QPochhammer[q]^3 + 9*q*QPochhammer[q^9]^3) / QPochhammer[q^3];

CubicThetaA /: HoldPattern[Derivative[n_][CubicThetaA]] /; n > 0 := Derivative[n - 1][CubicThetaAPrime];

CubicThetaB[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
QPochhammer[q]^3 / QPochhammer[q^3];

CubicThetaB /: HoldPattern[Derivative[n_][CubicThetaB]] /; n > 0 := Derivative[n - 1][CubicThetaBPrime];

CubicThetaC[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
3 * q^(1/3) * QPochhammer[q^3]^3 / QPochhammer[q];

CubicThetaC /: HoldPattern[Derivative[n_][CubicThetaC]] /; n > 0 := Derivative[n - 1][CubicThetaCPrime];

CubicThetaAPrime[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
With[{a = CubicThetaA[q], c = CubicThetaC[q]},
	(c^3/3 - a^3/12 + a*WeierstrassZeta[1, WeierstrassInvariants[{1, -I*Log[q]/(2*Pi)}]]/Pi^2)/q
];

CubicThetaAPrime /: HoldPattern[Derivative[1][CubicThetaAPrime]] :=
(With[{a = CubicThetaA[#], b = CubicThetaB[#], c = CubicThetaC[#], ap = CubicThetaAPrime[#]},
	2*b^3*c^3/(a*(3*#)^2) - ap/# + 2*ap^2/a
] &);

CubicThetaBPrime[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
With[{a = CubicThetaA[q], b = CubicThetaB[q]},
	b*(-a^2/12 + WeierstrassZeta[1, WeierstrassInvariants[{1, -I*Log[q]/(2*Pi)}]]/Pi^2)/q
];

CubicThetaBPrime /: HoldPattern[Derivative[1][CubicThetaBPrime]] :=
(With[{a = CubicThetaA[#], b = CubicThetaB[#], c = CubicThetaC[#], bp = CubicThetaBPrime[#]},
	-a*b*c^3/(3*#)^2 - bp/# + 2*bp^2/b
] &);

CubicThetaCPrime[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
With[{a = CubicThetaA[q], c = CubicThetaC[q]},
	c*(a^2/4 + WeierstrassZeta[1, WeierstrassInvariants[{1, -I*Log[q]/(2*Pi)}]]/Pi^2)/q
];

CubicThetaCPrime /: HoldPattern[Derivative[1][CubicThetaCPrime]] :=
(With[{a = CubicThetaA[#], b = CubicThetaB[#], c = CubicThetaC[#], cp = CubicThetaCPrime[#]},
	-a*b^3*c/(3*#)^2 - cp/# + 2*cp^2/c
] &);

End[];

SetAttributes[{CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime}, {Listable, NumericFunction, ReadProtected}];

Protect[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];

EndPackage[];

Last edited by lanxiyu (2024-09-28 22:03:52)