getData[dir_]:=Dataset[With[{data=Import[#],metadata=ToExpression/@Rest[StringSplit[FileBaseName[#],"_"]]},<|"n"->metadata[[1]],"β0"->metadata[[2]],"β1"->metadata[[3]],"N"->metadata[[4]],"βs"->data[[1]],"f"->data[[2]],"ps"->data[[3]],"Es"->data[[4]],"Cs"->((data[[1]])^2)data[[5]],"ρs"->(#/Total[#]&)/@data[[6;;]]|>]&/@FileNames[FileNameJoin[{dir,"hmm_*.dat"}]]]

c[βi_,βj_][A_]:=2(((((βi)^2)+((βj)^2)))/(((βi-βj))^2))((InverseErfc[A])^2)

Ts[dat_]:=dat["βs"]

Cstmp[dat_]:=With[{ts=Ts[dat]},{{#,#+ts[[2]]-ts[[1]]}&/@ts,dat["ps"]}]

Cs[dat_]:=With[{ts=Ts[dat]},MapThread[{Mean[{#1[[1]],#1[[2]]}],c[#1[[1]],#1[[2]]][#2]}&,{MapThread[#1+{0,#2}&,{Most[ts],Differences[ts]}],dat["ps"]}]]

logZ0[n_]:=(n-1)(n/4+1)Log[2]+Sum[(k/2)Log[π]-Log[Gamma[k/2]],{k,2,n}]

logZ[n_][β_]:=logN[n]+n(n-1)/4Log[((2π E)/(n β))]

ΔSs[dat_]:=With[{ts=Ts[dat],cs=Cs[dat]},{Mean[#1],((#2(#1[[2]]-#1[[1]]))/(Mean[#1]))}&@@@Thread[{Rest@Rest[FoldList[{#1[[2]],#2}&,{0,0},ts]],Cs[dat][[All,2]]}]]

Cf[dat_][x_]:=Interpolation[Normal@Prepend[Thread[{Normal[dat["βs"]],Normal[dat["Cs"]]/Normal[dat["βs"]]/((dat["n"])^2)}],{0,0}],Method->"Spline"][x]

oldSs[dat_]:=With[{ss=ΔSs[dat]},Thread[{ss[[All,1]],Accumulate[ss[[All,2]]]}]]

Sf[dat_][x_]:=Module[{y},With[{c=Cf[dat][y]},Interpolation[{#,NIntegrate[c,{y,0,#}]}&/@Normal[dat["βs"]]][x]]]

ff[dat_,Δx_,δx_]:=Module[{x},Module[{n=dat["n"],β,xmax=dat["β1"]},
Solve[(Sf[dat][x]/.{x->xmax})==(logZ0[n]-logZ[n][xmax])/(n^2),logN[n]]
]]

ns=Import["https://oeis.org/A206711/b206711.txt","Table"][[2;;,2]]/2;

d=getData["~/doc/research/hadamard/data"];

d[All,{#n,#N,#β1}&]

sample8=d[SelectFirst[#n==8&&#N==128&]];
sample12=d[Select[#n==12&&#N==128&&#β1==50&]][[2]];
sample16=d[SelectFirst[#n==16&&#N==128&&#β1==50&]];
sample20=d[Select[#n==20&&#N==256&]][[3]];

c8[x_]=Cf[sample8][x];
c12[x_]=Cf[sample12][x];
c16[x_]=Cf[sample16][x];
c20[x_]=Cf[sample20][x];

Plot[{(1-3/π)(x^2),1/4,x c8[x],x c12[x],x c16[x],x c20[x],x c24[x]},{x,0.0001,8},AxesLabel->{β,"c(β)"},PlotRange->{All,{0,3.8}},PlotLegends->{"High-T expansion","Low-T expansion",n==8,n==12,n==16,n==20},PlotStyle->Join[{Black,{Black,Dashed}},ColorData[97]/@Range[5]]]

Show[LogLogPlot[{(1-3/π)x,1/(4x),c8[x],c12[x],c16[x],c20[x],c24[x]},{x,0.00001,50},AxesLabel->{β,"c(β) / β"},PlotLegends->{"High-T expansion","Low-T expansion",n==8,n==12,n==16,n==20},PlotRange->{All,{Automatic,1}},PlotStyle->Join[{Black,{Black,Dashed}},ColorData[97]/@Range[5]]]]

s8[x_]=Sf[sample8][x];
s12[x_]=Sf[sample12][x];
s16[x_]=Sf[sample16][x];
s20[x_]=Sf[sample20][x];

fit8=ff[sample8,1,0.01];
fit12=ff[sample12,1,0.01];
fit16=ff[sample16,1,0.01];
fit20=ff[sample20,1,0.01];

LogLinearPlot[Evaluate[{s8[x],(logZ0[8]-logZ[8][x])/(8^2)/.fit8[[1]],s12[x],(logZ0[12]-logZ[12][x])/(12^2)/.fit12[[1]],s16[x],(logZ0[16]-logZ[16][x])/(16^2)/.fit16[[1]],s20[x],(logZ0[20]-logZ[20][x])/(20^2)/.fit20[[1]]}],{x,1,10},AxesLabel->{β,(s[0]-s[β])},PlotStyle->Thread[{Flatten[ConstantArray[{None,Dashed},4]],ColorData[97]/@Riffle[Range[4],Range[4]]}],PlotLegends->{n==8,None,n==12,None,n==16,None,n==20}]

Exp[logN[8]]/ns[[2]]/.fit8[[1]]

Exp[logN[12]]/ns[[3]]/.fit12[[1]]

Exp[logN[16]]/ns[[4]]/.fit16[[1]]

ns[[5]]/Exp[logN[20]]/.fit20[[1]]

(Sf[sample20][50]-(logZ0[20]-logZ[20][50])/(20^2))/((logZ0[20]-logZ[20][50])/(20^2))/.logN[20]->Log[ns[[5]]]

ρn=0;
ρ8=sample8[Thread[{#βs,5000/(2ρn+1)Total/@#ρs[[All,5001-ρn;;5001+ρn]]/Sqrt@#n}]&];
ρ12=sample12[Thread[{#βs,5000/(2ρn+1)Total/@#ρs[[All,5001-ρn;;5001+ρn]]/Sqrt@#n}]&];
ρ16=sample16[Thread[{#βs,5000/(2ρn+1)Total/@#ρs[[All,5001-ρn;;5001+ρn]]/Sqrt@#n}]&];
ρ20=sample20[Thread[{#βs,5000/(2ρn+1)Total/@#ρs[[All,5001-ρn;;5001+ρn]]/Sqrt@#n}]&];

ρPlotOld=Show[ListLogPlot[{ρ8,ρ12,ρ16,ρ20,{{0,0}}},PlotRange->{{0,24},Automatic},Joined->True,PlotLegends->{N==8,N==12,N==16,N==20,"(Sqrt[β/π])(E^(-β))"},PlotStyle->Join[ColorData[97]/@Range[4],{{Black,Dashed}}],AxesLabel->{β,"(Subscript[ρ, 0])(β) /(Sqrt[N])"}],
LogPlot[{ Exp[-β](Sqrt[β/π])},{β,0,25},PlotStyle->{Black,Dashed}]]

E8=sample8[Thread[{#βs,#Es/((#n)^2)}]&];
E12=sample12[Thread[{#βs,#Es/((#n)^2)}]&];
E16=sample16[Thread[{#βs,#Es/((#n)^2)}]&];
E20=sample20[Thread[{#βs,#Es/((#n)^2)}]&];

ListPlot[{E8,E12,E16,E20},PlotRange->{{0,10},All},Joined->True,PlotLegends->{n==8,n==12,n==16},AxesLabel->{β,e[β]}]

colors={Blue,Magenta,Red};

cPlot=LogPlot[{x c12[x],x c16[x],x c20[x]},{x,0.0001,8},Frame->True,FrameStyle->Black,FrameLabel->{None,Style[c,16]},RotateLabel->False,PlotStyle->colors,PlotRange->{{-0.2,8.2},{.05,4}},FrameTicks->{{{.1,1},None},{{0,2,4,6,8},None}}]

ρPlot=ListLogPlot[{ρ12,ρ16,ρ20},Joined->True,Frame->True,FrameStyle->Black,FrameLabel->{Style[β,16],Style[(Subscript[ρ, 0]),16]},RotateLabel->False,PlotStyle->colors,FrameTicks->{{{.001,.01,.1},None},{{0,2,4,6,8},None}},PlotRange->{{-0.2,8.2},{0.0005,0.5}}]

cρPlot=GraphicsColumn[{cPlot,ρPlot}]

Export["~/doc/research/hadamard/paper/equivert.eps",cρPlot]