(*This calculates the bk, ck, and dk for the antisymmetric OPs, plus some graphs of these values *)

(* If you don't want to recalculate the values of the OP (takes over an hour), skip down to the section that begins "Start Here" *)

(*number of iterations for the boundary values of the Pj, take a little more than twice the number of OPs you want to calculate *)
n=210; 
(*initial data for the Pjs*)
m=5;
Array[alpha,n];
Array[beta,n];
Array[nu,n];
Array[gamma,n];
alpha[0]=1;
beta[0]=-1/2;
nu[0]=0;
alpha[1]=1/6;
gamma[0]=1/2;
(* This calculates the boundary values of the Pj1*)
Do[alpha[j] = (4/(5^j-5))*Sum[alpha[j-l]*alpha[l], {l,1,j-1}],{j,2,n,1}];

(*boundary values of the Pj2 *)
Do[beta[j] = (2/(15*(5^j-1)))*Sum[ (3*5^(j-l) - 5^(l+1)+6)*alpha[j-l]*beta[l],{l,0,j-1}],{j,1,n,1}];

(*boundary values of the normal derivatives of the Pj1*)Do[nu[j]=((5^j+1)/2)*alpha[j]+2*Sum[nu[l]*beta[j-l],{l,1,j-1}],{j,1,n,1}];

(*boundary values of the Pj3*)Do[gamma[j]=3*alpha[j+1],{j, 1,n,1}];

(*This cell creates the inner product matrix for the recurrence relation*)
(*Let Num be the number of OPs you want to make*)
Num=100;
Symmetrize[mat_]:=mat+Transpose[mat]-DiagonalMatrix[Flatten[MapIndexed[#1[[#2]]&,mat]]];
Mn=Table[0,{j,0,Num},{k,0,Num}];
Do[Mn[[j+1,k+1]]=18*Sum[ alpha[j-l+1]*nu[k+l+2]-alpha[k+l+2]*nu[j-l+1],{l,0,j}], {j,0,Num},{k,j,Num}];
Mn=Symmetrize[Mn];

(*Below are the initial values for the recurrence for antisymmertric OPs*)
omegas=Table[0,{j,1,Num},{i,1,Num}];
ds = Table[0,{j,1,Num+1}];
bs = Table[0,{j,1,Num}];
cs =Table[0,{j,1,Num}];
fs = Table[0,{j,1,Num}];
zeta = Table[0,{j,1,Num}];
zeta[[1]] = -1/6;
ds[[1]]=30;
ds[[2]]=(271188000/449);
bs[[1]] =-(79/540);
(*The omegas matrix is the matrix of coefficents of OPk in terms of the Pj.  omegas[[k,j]] is the coefficient of Pj for OPk *)
omegas[[1,1]] = 1;
omegas[[2,1]] = -(11/540);
omegas[[2,2]] = 1;

(*Recurrence for the antisymmetric omegas, bk, ck, dk, and gammas*)
Do[
If[j==2, h = Table[0,{j,1,Num}] , h =omegas[[j-2]]];
g = omegas[[j-1]];
f = Table[omegas[[j-1,k-1]],{k,Num}];
f[[1]]=zeta[[j-1]];
(*omegas[[j]] = N[(f - bs[[j-1]]*g - cs[[j-1]]*h),100000];*)
omegas[[j]] = (f - bs[[j-1]]*g - cs[[j-1]]*h);

zeta[[j]] = -Sum[omegas[[j,g]]*gamma[g],{g,1,j}]/gamma[0];
v=Table[If[k==1,zeta[[j]],omegas[[j,k-1]]],{k,Num+1}];
u = Table[If[k==Num+1,0,omegas[[j,k]]],{k,Num+1}];
w=Mn.u;
ds[[j]] = 1/(u.w );
bs[[j]] = ds[[j]]*v.w;
cs[[j]] = ds[[j-1]]/ds[[j]];
Print[j];
,{j,2,Num}];
MatrixForm[N[omegas],100000];
MatrixForm[N[ds],100000];

(* Save values from Recurrence so we can use them later for plotting *)DumpSave["pk-recursion.mx", "Global`"]

(*Saving values to file for use in Matlab code (i.e. plotting).  The function N[] rounds off the the [blank] digit *)
Export["pk-recursion.mx",N[omegas],"Table"];
Export["pk-bs",N[bs],"Table"];
Export["pk-cs",N[cs],"Table"];
Export["pk-ds",N[ds],"Table"];

(* Start Here *)
(*if you only want to make graphs, and not recalculate the bk, ck, dk and omegas*)
<<"pk-recursion.mx"

(* Determinant of the Jacobi Matrix *)
DJ = Table[0,{j,1,Num}];
DJ[[1]]=1;
DJ[[2]]=bs[[1]];
Do[DJ[[j]] = bs[[j]]*DJ[[j-1]] - cs[[j]]*DJ[[j-2]],{j,3,Num,1}]; 
LDJ = Log[Abs[DJ]];
min=Min[Abs[DJ]];

(* Plot of the Det (Jk) in log scale *)
pLDJ1=ListPlot[LDJ,PlotLabel->"sequence of Log( Det(Jk) ) for antisymmetric OP",(*PlotRange->{{0,Num},{0,-2500}},*)FrameLabel->{"k","Det (Subscript[J, k])"},RotateLabel -> True,PlotMarkers-> {Automatic,7},FrameStyle->Directive[Black,FontSize->12],ImageSize -> {400},Frame -> True]
pDJ1=ListLogPlot[Abs[DJ],PlotLabel->"sequence of Det (Subscript[J, k]) for Antisymmetric OP",(*PlotRange->{{1,Num},{1,min}},*)FrameLabel->{"k","Det (Subscript[J, k])"},RotateLabel -> True,PlotMarkers-> {Automatic,7},FrameStyle->Directive[Black,FontSize->12],ImageSize -> {400},Frame -> True]
Export["pk_logDetJk.eps",pDJ1];
Export["pk_logDetJk.png",pDJ1];

(* Plots ratio of d(k+1)^2 to dk^2 *)
ratio = Table[ds[[j+1]]/ds[[j]],{j,1,Num-1}];
ratioinv = Table[ds[[j]]/ds[[j+1]],{j,1,Num-1}];
(*ratiol = Log[ratio];*)
p0=ListLogPlot[ratio,PlotLabel->"Sequence of ((Subscript[d, k + 1])^2)/((Subscript[d, k])^2)  (Antisymmetric)",Frame -> True,FrameLabel->{Text[k],Text[((Subscript[d, k+1])^2)/((Subscript[d, k])^2)]},RotateLabel -> True,ImageSize-> 400,FrameStyle->Directive[Black,FontSize->12]]
p0inv=ListLogPlot[ratioinv,PlotLabel->"Sequence of ((Subscript[d, k])^2)/((Subscript[d, k + 1])^2)  (Antisymmetric)",Frame -> True,FrameLabel->{Text[k],Text[((Subscript[d, k])^2)/((Subscript[d, k+1])^2)]},RotateLabel -> True,ImageSize-> 400,FrameStyle->Directive[Black,FontSize->12]]
Export["Mathematica Output\\pkdkratio1.eps",p0];
Export["Mathematica Output\\pkdkratio2.eps",p0inv];
Export["Mathematica Output\\pkdkratio2.png",p0inv];

(* Plots of the bk, ck, and dk in log scale *)
p1=ListLogPlot[ds,PlotLabel->"Sequence of Log[((Subscript[d, k])^2)] (Antisymmetric)",Frame -> True,FrameLabel->{Text[k],Text[Log[((Subscript[d, k])^2)]]},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],ImageSize-> 400]
p2=ListLogPlot[-bs,PlotLabel->"Sequence of Log[-(Subscript[b, k])] (Antisymmetric)",Frame -> True,FrameLabel->{Text[k],Text[Log[-(Subscript[b, k])]]},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],ImageSize-> 400]
p3=ListLogPlot[cs,PlotLabel->"Sequence of Log[(Subscript[c, k])] (Antisymmetric)",Frame -> True,FrameLabel->{Text[k],Text[Log[(Subscript[c, k])]]},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],ImageSize-> 400]
(*scs = Table[Sqrt[cs[[j]]],{j,1,Num}];
p4=ListPlot[scs,PlotLabel->"sequence of sqrt[ck]",PlotRange->{{0,Num},{0,.0001}},AxesOrigin->{0,0},AxesLabel->{"k","Sqrt[ck]"}]*)
Export["Mathematica Output\\pkds.eps",p1];
Export["Mathematica Output\\pkbs.eps",p2];
Export["Mathematica Output\\pkcs.eps",p3];
Export["Mathematica Output\\pkds.png",p1];
Export["Mathematica Output\\pkbs.png",p2];
Export["Mathematica Output\\pkcs.png",p3];
(*Export["pk_scs.png",p4];*)

(* Fits of Log(dk^2)/Log(k) *)dsk = Table[Log[ds[[j]]]/Log[j],{j,2,Num}];
N[line=LinearModelFit[dsk,{1,x},x]]
p5=Show[ListPlot[dsk,PlotLabel->"Linear Fit of Log[dk^2]/Log[k] (Antisymmetric)",Frame -> True,FrameLabel->{Text[k],Text[Log[((Subscript[d, k])^2)]/Log[k]]},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],ImageSize-> 400],Plot[line,{x,1,Num},PlotStyle -> Black]]
N[line[{"RSquared","ANOVATable"}]]
(*dsll = Table[Log[Log[ds[[j]]]],{j,1,Num}];
g = Fit[dsll,{1,Log[x]},x]
p6=Show[ListPlot[dsll,PlotLabel->"Logarithm Fit of Log[Log[dk^2]]",AxesLabel->{"k","Log[Log[dk^2]]"}],Plot[g,{x,1,Num}]]
(*line[{"Rsquared","FitResiduals","ANOVATable"}]*)*)
Export["pkdklinearfit.png",p5];
Export["pkdklinearfit.eps",p5];
(*Export["dklogfit.png",p6];*)

(* even and odd halves of the bk, ck, and dk *)
Num=200;
(*bsodd = Table[{2*j-1,Log[-bs[[2*j-1]]]},{j,1,Num/2}];
bseven = Table[{2*j,Log[-bs[[2*j]]]},{j,1,Num/2}];
csodd = Table[{2*j-1,Log[cs[[2*j-1]]]},{j,1,Num/2}];
cseven = Table[{2*j,Log[cs[[2*j]]]},{j,1,Num/2}];
dsodd = Table[{2*j-1,Log[ds[[2*j-1]]]},{j,1,Num/2}];
dseven = Table[{2*j,Log[ds[[2*j]]]},{j,1,Num/2}];*)

bsodd = Table[{2*j-1,-bs[[2*j-1]]},{j,1,Num/2}];
bseven = Table[{2*j,-bs[[2*j]]},{j,1,Num/2}];
csodd = Table[{2*j-1,cs[[2*j-1]]},{j,1,Num/2}];
cseven = Table[{2*j,cs[[2*j]]},{j,1,Num/2}];
dsodd = Table[{2*j-1,ds[[2*j-1]]},{j,1,Num/2}];
dseven = Table[{2*j,ds[[2*j]]},{j,1,Num/2}];
mb=Min[-bs];
Mb=Max[-bs];
mc=Min[cs];
Mc=Max[cs];
md=Min[ds];
Md=Max[ds];

pb=ListLogPlot[{{bsodd},{bseven}},PlotStyle->{Blue,Red}, Joined -> True,Frame -> True,FrameLabel->{"k",Text[Log[-(Subscript[b, k])]]},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],PlotLabel->"the odd and even sequence of (Subscript[b, k]) for the (Subscript[p, k]) (Antisymmetric)",PlotRange ->{{1,200},{mb,Mb}},ImageSize-> 400]
pc=ListLogPlot[{{csodd},{cseven}},PlotStyle->{Blue,Red}, Joined -> True,Frame -> True,FrameLabel->{"k",Text[Log[(Subscript[c, k])]]},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],PlotLabel->"the odd and even sequence of (Subscript[c, k]) for the (Subscript[p, k]) (Antisymmetric)",PlotRange ->{{1,200},{mc,Mc}},ImageSize-> 400]
pd=ListLogPlot[{{dsodd},{dseven}},PlotStyle->{Blue,Red},Joined -> True, Frame -> True,FrameLabel->{"k",Text[Log[((Subscript[d, k])^2)]]},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],PlotLabel->"the odd and even sequence of ((Subscript[d, k])^2) for the (Subscript[p, k]) (Antisymmetric)",PlotRange ->{{1,200},{md,Md}},ImageSize-> 400]
Export["pkbkoddeven.png",pb];
Export["pkckoddeven.png",pc];
Export["pkdkoddeven.png",pd];

(* Plots of the dynamics? of the boundary points of the OP *)
gammas = Table[gamma[j-1],{j,1,Num}];
bv = N[omegas.gammas,1000000];
dyn =Table[{bv[[j]],bv[[j+1]]},{j,1,Num-1}];
N[bvmin=Min[bv]]
N[bvmax=Max[bv]]

(* Logs of the boundary points of the OP*)
bvl = Log[bv];
dynl = Log[dyn];

(* Plots of the dynamics at the boundary points*)
p7=ListLogPlot[bv,PlotLabel-> "Dynamics: (Subscript[p, k])((Subscript[q, 1])) vs k",Frame -> True,FrameLabel->{"k","(Subscript[p, k])((Subscript[q, 1]))"},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],ImageSize-> 400]
p8=ListLogLogPlot[{dyn},PlotLabel-> "Dynamics: (Subscript[p, k + 1])((Subscript[q, 1])) vs (Subscript[p, k])((Subscript[q, 1]))",Frame -> True,FrameLabel->{"(Subscript[p, k])((Subscript[q, 1]))","(Subscript[p, k + 1])((Subscript[q, 1]))"},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],ImageSize-> 400,PlotRange ->{{bvmin,bvmax},{bvmin,bvmax}}]
Export["pkvsk.png",p7];
Export["pkdyn.png",p8];
Export["pkvsk.eps",p7];
Export["pkdyn.eps",p8];

plined=LinearModelFit[dynl,{1,x},x];
N[plined]
p9=Show[ListPlot[dynl,PlotLabel->"Linear Fit of Log[(Subscript[p, k])((Subscript[q, 1]))]/Log[(Subscript[p, k - 1])((Subscript[q, 1]))]",Frame -> True,FrameLabel->{"Log[(Subscript[p, k])((Subscript[q, 1]))]","Log[(Subscript[p, k - 1])((Subscript[q, 1]))]"},RotateLabel-> True, FrameStyle->Directive[Black,FontSize->12],ImageSize-> 400],Plot[plined[x],{x,-3000,Num},PlotStyle -> Black]]
N[plined[{"RSquared","ANOVATable"}]]
Export["pkdynlinearfit.eps",p9];
Export["pkdynlinearfit.png",p9];