"Computation of the heavy quark potential from classical statistical simulations of the Wilson loop.";

SetDirectory[NotebookDirectory[]]

FitFunc0[x_,G_,M_,A_,B_]:=A*((G*Cos[B]-(M-x)*Sin[B])/((x-M)^2+G^2));

FitFunc[x_,G_,M_,A_,B_,K_,L_,F_]:=A*((G*Cos[B]-(M-x)*Sin[B])/((x-M)^2+G^2))+K+L*(M-x)+F*(M-x)^2;

"The simulations according to which parameters shall used in the analysis:";

NumBeta=3;

Betas={4,6,8,10,12,16,20};
BetasChar = {"4","6", "8", "10","12", "16","20"};

Betas={12,16,20};
BetasChar = {"12","16","20"};

Nconfs={600,200,1200};

Nxs=Table[32,{k,1,NumBeta}]

Parameterdts=Table[0.1,{k,1,NumBeta}]

MeasureEach=Table[2,{k,1,NumBeta}]

dts=MeasureEach*Parameterdts

dtsChar=Table["01",{k,1,NumBeta}]

as=Table[1,{k,1,NumBeta}]

asChar=Table["1",{k,1,NumBeta}]

WlExtents=Table[7,{k,1,NumBeta}]

NJack =10;

"Check whether there are any missing runs";

For[BetaCnt=1,BetaCnt≤NumBeta, BetaCnt++,

(* Current set of paramters *)

CurBeta=Betas[[BetaCnt]];
CurBetasChar=BetasChar[[BetaCnt]];

Cura=as[[BetaCnt]];
CuraChar=asChar[[BetaCnt]];

Curdt=dts[[BetaCnt]];
CurdtChar=dtsChar[[BetaCnt]];

CurNx=Nxs[[BetaCnt]];
CurWlExtent=WlExtents[[BetaCnt]];
CurNconf=Nconfs[[BetaCnt]];

Print["Checking for Beta: " , CurBeta," whether runs are missing"];

For[r=1,r≤CurWlExtent,r++, Print["R=",r,"\n"];MissingRuns[r]={};For[RunCnt=0,RunCnt≤CurNconf-1,RunCnt++,If[Not[FileExistsQ[ToString@StringForm["../ColSource/V``B``/R``/results_THRM_S0_S``_N``_a``_dt``_B``_run``/WilsonLoop.dat",CurNx,CurBetasChar,r,r,CurNx,CuraChar,CurdtChar,CurBetasChar,RunCnt]]],AppendTo[MissingRuns[r],RunCnt],If[FileByteCount[ToString@StringForm["../ColSource/V``B``/R``/results_THRM_S0_S``_N``_a``_dt``_B``_run``/WilsonLoop.dat",CurNx,CurBetasChar,r,r,CurNx,CuraChar,CurdtChar,CurBetasChar,RunCnt]]<293000,AppendTo[MissingRuns[r],RunCnt]] ]
];
Print[MissingRuns[r]]
];
];

"Check the whether all runs have finished, check out fitting ranges";

WindowingUpToRTmp=Table[2,{k,1,NumBeta}]

For[BetaCnt=1,BetaCnt≤NumBeta, BetaCnt++,

(* Current set of paramters *)

CurBeta=Betas[[BetaCnt]];
CurBetasChar=BetasChar[[BetaCnt]];

Cura=as[[BetaCnt]];
CuraChar=asChar[[BetaCnt]];

Curdt=dts[[BetaCnt]];
CurdtChar=dtsChar[[BetaCnt]];

CurNx=Nxs[[BetaCnt]];

CurNconf=Nconfs[[BetaCnt]];

CurWindowUpToR=WindowingUpToRTmp[[BetaCnt]];
CurReadEveryTimeStep=2;
CurWlExtent=WlExtents[[BetaCnt]];

ClearAll[CurRawWilsonData];

Print["Checking for Beta: " , CurBeta," whether all runs have finished"];

CurRawWilsonData=ParallelTable[Table[Import[ToString@StringForm["../ColSource/V``B``/R``/results_THRM_S0_S``_N``_a``_dt``_B``_run``/WilsonLoop.dat",CurNx,CurBetasChar,r,r,CurNx,CuraChar,CurdtChar,CurBetasChar,RunCnt],"Table"][[2;;;;CurReadEveryTimeStep]],{RunCnt,0,CurNconf-1}],{r,1,CurWlExtent}];

Print["Beta: ", CurBeta, " checking for missing runs"];

For[r=1,r≤CurWlExtent,r++, Print["R=",r,"\n"];MissingRuns[r]={};For[RunCnt=1,RunCnt≤CurNconf,RunCnt++,If[Not[Dimensions[CurRawWilsonData[[r]][[RunCnt]]][[1]]==6000/CurReadEveryTimeStep+1],AppendTo[MissingRuns[r],RunCnt]] ;
];
Print[MissingRuns[r]]
];

CurWilsonLoops=ParallelTable[Map[{#[[2]],#[[3]]+I*#[[4]],#[[4+3]]+I*#[[4+4]]}&,Transpose[Join[Transpose[ Mean[CurRawWilsonData[[r]]] ],Transpose[StandardDeviation[CurRawWilsonData[[r]]]/Sqrt[CurNconf] ]]] ] ,{r,1,CurWlExtent}];

Print[GraphicsColumn[Table[  ListPlot[ Map[{Re[#[[1]]]/Curdt,Around[Re[#[[2]]],Re[#[[3]]]]}&,CurWilsonLoops[[r]][[;;All;;4,{1,2,3}]]]
,AxesLabel->{"t","Re[W]"},AxesStyle->{Gray, AbsoluteThickness[3],FontSize->20},ImageSize->800,PlotLabel->Text[Style["Real-part of the DFT of the Wilson Loop",FontSize->20]],AxesOrigin->{0,0},PlotRange->All]  ,{r,1,CurWindowUpToR+1}  ]]];

(*Print[Show[ListPlot[Table[Map[{#[[1]]/Curdt,Around[Re[#[[2]]],Re[#[[3]]]]}&,CurWilsonLoops[[r]][[;;All;;4,{1,2,3}]]],{r,1,1}]
,Joined->True,AxesLabel->{"lines in file","Re[W]"},AxesStyle->{{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20}},ImageSize->800,PlotLabel->Text[Style["Real-part of the Wilson Loop",FontSize->20]]],PlotRange->All]];*)
 
CurDFTLen=(Length[CurRawWilsonData[[1]][[1]]]-1)*2;
Print[CurDFTLen];
 Frequencies=Table[If[k≤CurDFTLen/2,(Pi/(Curdt*CurReadEveryTimeStep))*k/(CurDFTLen/2),(Pi/(Curdt*CurReadEveryTimeStep))*(-CurDFTLen+k)/(CurDFTLen/2)],{k,0,CurDFTLen-1}];

CurRawWilsonDataDFT=ParallelTable[Table[  Table[Exp[I*l*Pi],{l,0,CurDFTLen-1}]*Fourier[Join[ Reverse[Map[#[[3]]+I*#[[4]]&,CurRawWilsonData[[r]][[RunCnt]]]][[;;-2]],Map[#[[3]]+I*#[[4]]&,CurRawWilsonData[[r]][[RunCnt]]][[;;-2]] ]*If[r<=CurWindowUpToR,Table[Sin[Pi*x/CurDFTLen]^2,{x,0,CurDFTLen-1} ],1]   ],{r,1,CurWlExtent}],{RunCnt,1,CurNconf}];

TmpMeanDFTRaw=Mean[CurRawWilsonDataDFT];
    TmpReStdDevDFTRaw=StandardDeviation[Re[CurRawWilsonDataDFT]];
    TmpImStdDevDFTRaw=StandardDeviation[Im[CurRawWilsonDataDFT]];

CurWilsonLoopsDFT=ParallelTable[MapThread[{#1+I*#1,#2,#3+I*#4}&,{Frequencies,TmpMeanDFTRaw[[r]],TmpReStdDevDFTRaw[[r]]/Sqrt[CurNconf],TmpImStdDevDFTRaw[[r]]/Sqrt[CurNconf]}] ,{r,1,CurWlExtent}];

Print[GraphicsColumn[Table[  ListPlot[ Map[{Re[#[[1]]]/Frequencies[[2]],Around[Re[#[[2]]],Re[#[[3]]]]}&,CurWilsonLoopsDFT[[r]][[;;100,{1,2,3}]]]
,AxesLabel->{"t","Re[W]"},AxesStyle->{Gray, AbsoluteThickness[3],FontSize->20},ImageSize->800,PlotLabel->Text[Style["Real-part of the DFT of the Wilson Loop",FontSize->20]],AxesOrigin->{0,0},PlotRange->{{0,100},All}]  ,{r,1,CurWlExtent}  ]]];

]

"Select the Read-in and the fitting ranges";

MaxReadIn={All,All,All};

ReadEveryTimeStep=Table[2,{k,1,NumBeta}]

WindowingUpToR={2,2,2}

FitRanges={{{10,12, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15},
                      {30, 35, 40, 43,45, 45, 45, 45, 45, 45, 45, 45}},{{15,15, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18},{28, 40, 40, 37,42, 42, 42, 42, 42, 42, 42, 42}},{{15,13, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15},
{27, 38, 38, 40,42, 42, 45, 45, 45, 45, 45, 45}}}

"Compute expectation values of Basic observables";

For[BetaCnt=1,BetaCnt≤NumBeta, BetaCnt++,

(* Current set of paramters *)

CurBeta=Betas[[BetaCnt]];
CurBetasChar=BetasChar[[BetaCnt]];

Cura=as[[BetaCnt]];
CuraChar=asChar[[BetaCnt]];

Curdt=dts[[BetaCnt]];
CurdtChar=dtsChar[[BetaCnt]];

CurNx=Nxs[[BetaCnt]];

CurNconf=Nconfs[[BetaCnt]];

CurMaxReadIn=MaxReadIn[[BetaCnt]];
CurWlExtent=WlExtents[[BetaCnt]];

Print["Basic Observables for Beta: ", CurBeta];

(* Read in plaquette, Energy and Gauss Law at intermediate r=Nx/4 *)

CurRawBasicData=Table[Import[ToString@StringForm["../ColSource/V``B``/R``/results_THRM_S0_S``_N``_a``_dt``_B``_run``/BasicOvservables.dat",CurNx,CurBetasChar,Min[Max[WlExtents],Floor[CurNx/4]],Min[Max[WlExtents],Floor[CurNx/4]],CurNx,CuraChar,CurdtChar,CurBetasChar,RunCnt],"Table"],{RunCnt,0,CurNconf-1}];

CurBasicData=Map[{#[[2]],#[[3]],#[[7+3]],#[[4]],#[[7+4]],#[[3]]+#[[4]],Sqrt[#[[7+3]]^2+#[[7+4]]^2],#[[7]],#[[7+7]]}&,Transpose[Join[Transpose[Mean[CurRawBasicData]],Transpose[StandardDeviation[CurRawBasicData]/Sqrt[CurNconf]]]] ] ;

(* Output Energy *)

(*Print[Show[ListPlot[{CurBasicData[[;;;;4,{1,2}]],CurBasicData[[;;;;10,{1,4}]]},Frame->{True,True,False,False}, Axes->{False,False,False,False},FrameLabel->{"t","E"},FrameStyle->{{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20}},ImageSize->800,PlotStyle->AbsolutePointSize[10],PlotMarkers->{Automatic,12},PlotLabel->Text[Style["Electric and Magnetic Energy",FontSize->20]]],PlotRange->All]];

Print[Show[ListPlot[CurBasicData[[;;;;4,{1,6}]],Frame->{True,True,False,False}, Axes->{False,False,False,False},FrameLabel->{"t","E"},FrameStyle->{{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20}},ImageSize->800,PlotStyle->AbsolutePointSize[10],PlotMarkers->{Automatic,12},PlotLabel->Text[Style["Total Energy",FontSize->20]]],PlotRange->All]];*)

Quiet[DeleteFile[ToString@StringForm["V``B``ColSourceNrg.dat",CurNx,CurBetasChar]]];
EnergyOutput=OpenAppend[ToString@StringForm["V``B``ColSourceNrg.dat",CurNx,CurBetasChar],FormatType->StandardForm,PageWidth->Infinity];

For[i=1,i<Length[CurBasicData],i++,

Write[EnergyOutput,CForm[ CurBasicData[[i,1 ]] ], " " , CForm[ CurBasicData[[i,2 ]] ], " ", CForm[ CurBasicData[[i,3]] ], " " ,CForm[CurBasicData[[i,4 ]] ], " ", CForm[ CurBasicData[[i,5 ]] ]    ];
];

Close[EnergyOutput];

(* Output Gauss Law *)

(*Print[Show[ListPlot[CurBasicData[[;;;;10,{1,8}]],Frame->{True,True,False,False}, Axes->{False,False,False,False},FrameLabel->{"t","ΔG"},FrameStyle->{{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20}},ImageSize->800,PlotStyle->AbsolutePointSize[10],PlotMarkers->{Automatic,12},PlotLabel->Text[Style["Deviation from Gauss Law",FontSize->20]]],PlotRange->All]];*)

Quiet[DeleteFile[ToString@StringForm["V``B``ColSourceGaussLaw.dat",CurNx,CurBetasChar]]];
GaussLawOutput=OpenAppend[ToString@StringForm["V``B``ColSourceGaussLaw.dat",CurNx,CurBetasChar],FormatType->StandardForm,PageWidth->Infinity];

For[i=1,i<Length[CurBasicData],i++,

Write[GaussLawOutput,CForm[ CurBasicData[[i,1 ]] ], " " , CForm[ CurBasicData[[i,8 ]] ], " ", CForm[ CurBasicData[[i,9]] ]    ];
];

Close[GaussLawOutput];

(* Output Thermalization steps *)

CurRawThrmControlData=Table[Import[ToString@StringForm["../ColSource/V``B``/R``/results_THRM_S0_S``_N``_a``_dt``_B``_run``/MooreThermalizationCntrol.dat",CurNx,CurBetasChar,Min[Max[WlExtents],Floor[CurNx/4]],Min[Max[WlExtents],Floor[CurNx/4]],CurNx,CuraChar,CurdtChar,CurBetasChar,RunCnt],"Table"][[2;;]],{RunCnt,0,10}];

(*Print[ListPlot[Table[CurRawThrmControlData[[RunCnt]][[All,{1,4}]],{RunCnt,1,CurNconf}],Frame->{True,True,False,False}, Axes->{False,False,False,False},FrameLabel->{"t","E"},FrameStyle->{{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20}},ImageSize->800,PlotStyle->Table[{AbsolutePointSize[10]},{RunCnt,1,CurNconf}],PlotMarkers->{Automatic,12},PlotLabel->Text[Style["Thermalization Check",FontSize->20]]]];*)

Quiet[DeleteFile[ToString@StringForm["V``B``ColSourceTherm.dat",CurNx,CurBetasChar]]];
ThermOutput=OpenAppend[ToString@StringForm["V``B``ColSourceTherm.dat",CurNx,CurBetasChar],FormatType->StandardForm,PageWidth->Infinity];
For[r=1,r≤10,r++,
For[i=1,i≤19,i++,

Write[ThermOutput, r, " ",CForm[ CurRawThrmControlData[[r]][[i,1]] ], " " , CForm[ CurRawThrmControlData[[r]][[i,4]] ]      ];
]
];
Close[ThermOutput];

]

"Carry out the analysis";

For[ReV={};ImV={};ReVBW={};ImVBW={};ImVExp={};BetaCnt=3,BetaCnt≤NumBeta, BetaCnt++,

(* Current set of paramters *)

CurBeta=Betas[[BetaCnt]];
CurBetasChar=BetasChar[[BetaCnt]];

Cura=as[[BetaCnt]];
CuraChar=asChar[[BetaCnt]];

Curdt=dts[[BetaCnt]];
CurdtChar=dtsChar[[BetaCnt]];

CurNx=Nxs[[BetaCnt]];

CurNconf=Nconfs[[BetaCnt]];

CurWindowUpToR=WindowingUpToR[[BetaCnt]];
CurReadEveryTimeStep=ReadEveryTimeStep[[BetaCnt]];
CurWlExtent=WlExtents[[BetaCnt]];

CurMaxReadIn=MaxReadIn[[BetaCnt]];

Print["Beta: ", CurBeta];

(* Read in the raw Wilson loop data *)

CurJackNumSamp=CurNconf-CurNconf/NJack;

ClearAll[CurRawWilsonData];

CurRawWilsonData=ParallelTable[Table[Import[ToString@StringForm["../ColSource/V``B``/R``/results_THRM_S0_S``_N``_a``_dt``_B``_run``/WilsonLoop.dat",CurNx,CurBetasChar,r,r,CurNx,CuraChar,CurdtChar,CurBetasChar,RunCnt],"Table"][[2;;CurMaxReadIn;;CurReadEveryTimeStep]],{RunCnt,0,CurNconf-1}],{r,1,CurWlExtent}];

Print["Wilson Loop Raw"];

(*Print[CurRawWilsonData];*)

CurDFTLen=(Length[CurRawWilsonData[[1]][[1]]]-1)*2;
Print[CurDFTLen];

ClearAll[CurRawWilsonDataDFT];

CurRawWilsonDataDFT=ParallelTable[Table[  Table[Exp[I*l*Pi],{l,0,CurDFTLen-1}]*Fourier[Join[ Reverse[Map[#[[3]]+I*#[[4]]&,CurRawWilsonData[[r]][[RunCnt]]]][[;;-2]],Map[#[[3]]+I*#[[4]]&,CurRawWilsonData[[r]][[RunCnt]]][[;;-2]] ]*If[r<=CurWindowUpToR,Table[Sin[Pi*x/CurDFTLen]^2,{x,0,CurDFTLen-1} ],1]      ],{r,1,CurWlExtent}],{RunCnt,1,CurNconf}];

Frequencies=Table[If[k≤CurDFTLen/2,(Pi/(Curdt*CurReadEveryTimeStep))*k/(CurDFTLen/2),(Pi/(Curdt*CurReadEveryTimeStep))*(-CurDFTLen+k)/(CurDFTLen/2)],{k,0,CurDFTLen-1}];

(*Print[Frequencies];*)

(* Compute the Jackknife subaverages mean and variance  *)

ClearAll[TmpMeanDFTRaw,TmpReStdDevDFTRaw,TmpImStdDevDFTRaw];
ClearAll[CurWilsonLoopsJack,CurWilsonLoopsDFTJack];

CurWilsonLoopsJack={};
CurWilsonLoopsDFTJack={};
SetSharedVariable[CurWilsonLoopsJack,CurWilsonLoopsDFTJack];

ParallelDo[

    CurJackWilsonLoopRaw=If[JckBin==1,Transpose[Transpose[CurRawWilsonData][[1+(CurNconf/NJack);;]]],
       If[JckBin==NJack,Transpose[Transpose[CurRawWilsonData][[1;;(NJack-1)*(CurNconf/NJack)]]],
    Join[  Transpose[Transpose[CurRawWilsonData][[1;;(JckBin-1)*(CurNconf/NJack)]]],
  Transpose[Transpose[CurRawWilsonData][[(JckBin)*(CurNconf/NJack)+1;;]]]  ] ]];

	Print["Jackknife Bin :", JckBin];

    AppendTo[CurWilsonLoopsJack,
      Table[Map[{#[[2]],#[[3]]+I*#[[4]],#[[4+3]]+I*#[[4+4]]}&,
 Transpose[Join[Transpose[ Mean[CurJackWilsonLoopRaw[[r]]] ],Transpose[StandardDeviation[CurJackWilsonLoopRaw[[r]]]/Sqrt[CurJackNumSamp] ]]] ] ,{r,1,CurWlExtent}] ];

(*Print[CurWilsonLoopsJack];*)

     CurJackWilsonLoopDFTRaw=If[JckBin==1,CurRawWilsonDataDFT[[1+(CurNconf/NJack);;]],
    If[JckBin==NJack,CurRawWilsonDataDFT[[1;;(NJack-1)*(CurNconf/NJack)]],
  Join[  CurRawWilsonDataDFT[[1;;(JckBin-1)*(CurNconf/NJack)]],CurRawWilsonDataDFT[[(JckBin)*(CurNconf/NJack)+1;;]]  ] ]];

    TmpMeanDFTRaw=Mean[CurJackWilsonLoopDFTRaw];
    TmpReStdDevDFTRaw=StandardDeviation[Re[CurJackWilsonLoopDFTRaw]];
    TmpImStdDevDFTRaw=StandardDeviation[Im[CurJackWilsonLoopDFTRaw]];

     AppendTo[CurWilsonLoopsDFTJack,
    Table[MapThread[{#1+I*#1,#2,#3+I*#4}&,{Frequencies,TmpMeanDFTRaw[[r]],TmpReStdDevDFTRaw[[r]]/Sqrt[CurJackNumSamp], TmpImStdDevDFTRaw[[r]]/Sqrt[CurJackNumSamp]}] ,{r,1,CurWlExtent}]  ];

  Clear[TmpMeanDFTRaw,TmpReStdDevDFTRaw,TmpImStdDevDFTRaw];

(*Print[CurJackWilsonLoopDFTRaw];*)

,{JckBin,1,NJack}];

(*Compute the mean and variance/errors of the primary quantity Wilson loop & output to HD *)

CurWilsonLoops=ParallelTable[Map[{#[[2]],#[[3]]+I*#[[4]],#[[4+3]]+I*#[[4+4]]}&,Transpose[Join[Transpose[ Mean[CurRawWilsonData[[r]]] ],Transpose[StandardDeviation[CurRawWilsonData[[r]]]/Sqrt[CurNconf] ]]] ] ,{r,1,CurWlExtent}];

Quiet[DeleteFile[ToString@StringForm["V``B``ColSourceWilsonCoord.dat",CurNx,CurBetasChar]]];
WLCoordOutput=OpenAppend[ToString@StringForm["V``B``ColSourceWilsonCoord.dat",CurNx,CurBetasChar],FormatType->StandardForm,PageWidth->Infinity];

For[r=1,r<=Length[CurWilsonLoops],r++,
For[i=1,i<Length[CurWilsonLoops[[r]]],i++,

Write[WLCoordOutput, r, " ",CForm[ CurWilsonLoops[[r]][[i,1]] ], " " , CForm[ Re[CurWilsonLoops[[r]][[i,2]]] ], " ", CForm[ Re[CurWilsonLoops[[r]][[i,3]]] ], " " , CForm[ Im[CurWilsonLoops[[r]][[i,2]]] ], " ", CForm[ Im[CurWilsonLoops[[r]][[i,3]]] ]    ];
]
];

Close[WLCoordOutput];

(*Compute the mean and variance/errors of the primary quantity Wilson loop in Fourier space & output to HD *)

    TmpMeanDFTRaw=Mean[CurRawWilsonDataDFT];
    TmpReStdDevDFTRaw=StandardDeviation[Re[CurRawWilsonDataDFT]];
    TmpImStdDevDFTRaw=StandardDeviation[Im[CurRawWilsonDataDFT]];

CurWilsonLoopsDFT=ParallelTable[MapThread[{#1+I*#1,#2,#3+I*#4}&,{Frequencies,TmpMeanDFTRaw[[r]],TmpReStdDevDFTRaw[[r]]/Sqrt[CurNconf],TmpImStdDevDFTRaw[[r]]/Sqrt[CurNconf]}] ,{r,1,CurWlExtent}];

Clear[TmpMeanDFTRaw,TmpReStdDevDFTRaw,TmpImStdDevDFTRaw];

Quiet[DeleteFile[ToString@StringForm["V``B``ColSourceWilsonDFT.dat",CurNx,CurBetasChar]]];
WLDFTOutput=OpenAppend[ToString@StringForm["V``B``ColSourceWilsonDFT.dat",CurNx,CurBetasChar],FormatType->StandardForm,PageWidth->Infinity];

For[r=1,r<=Length[CurWilsonLoopsDFT],r++,
For[i=1,i<=Length[CurWilsonLoopsDFT[[r]]],i++,

Write[WLDFTOutput, r, " ",CForm[ Re[CurWilsonLoopsDFT[[r]][[i,1]]] ], " " , CForm[ Re[CurWilsonLoopsDFT[[r]][[i,2]]] ], " ", CForm[ Re[CurWilsonLoopsDFT[[r]][[i,3]]] ], " " , CForm[ Im[CurWilsonLoopsDFT[[r]][[i,2]]] ], " ", CForm[ Im[CurWilsonLoopsDFT[[r]][[i,3]]] ]    ];

]
];

Close[WLDFTOutput];

Print[Show[ListPlot[Table[Map[{#[[1]]/Curdt,Around[Re[#[[2]]],Re[#[[3]]]]}&,CurWilsonLoops[[r]][[;;All;;4,{1,2,3}]]],{r,1,4}]
,AxesLabel->{"t/dt","Re[W]"},AxesStyle->{{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20},{Gray, AbsoluteThickness[3],FontSize->20}},ImageSize->800,PlotLabel->Text[Style["Real-part of the Wilson Loop",FontSize->20]]],PlotRange->All]];

(* Carry out the Breit-Wigner fit to the DFT data as well as simple exponential fit to the original Wilson loop   *)


ClearAll[CurPeakFitJack,CurBWPeakFitJack];
CurBWPeakFitJack={};
CurPeakFitJack={};
SetSharedVariable[CurPeakFitJack,CurBWPeakFitJack];

ParallelDo[

Print["Starting Fits for Jackknife bin:",JckBin];

CurFitRangeLowByEye=FitRanges[[BetaCnt]][[1]];
CurFitRangeUpByEye=FitRanges[[BetaCnt]][[2]];

CurMaxPositionInt=Table[Position[Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[All,2]]  ],Max[ Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[CurFitRangeLowByEye[[r]];;CurFitRangeUpByEye[[r]],2]]  ] ]][[1]][[1]],{r,1,CurWlExtent}];
(*Print[CurMaxPositionInt];*)

CurFitRangeLow=CurMaxPositionInt-Table[7,{r,1,CurWlExtent}];
CurFitRangeUp=CurMaxPositionInt+Table[7,{r,1,CurWlExtent}];

CurMaxPos=Table[Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[  CurMaxPositionInt[[r]],1]]],{r,1,CurWlExtent}];
(*Print[CurMaxPos];*)

CurMaxVal=Table[Max[ Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[CurFitRangeLow[[r]];;CurFitRangeUp[[r]],2]]  ] ],{r,1,CurWlExtent}];;

CurBWOnlyPeakFitBin=Table[ Quiet[NonlinearModelFit[ Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[CurFitRangeLow[[r]];;CurFitRangeUp[[r]],{1,2}]]  ],{FitFunc0[x,G,M,A,B]},{{G,0.001},{M,CurMaxPos[[r]]},{A,1},{B,0.0000001}},x,Method->LevenbergMarquardt,MaxIterations->5000,Weights->1/Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[CurFitRangeLow[[r]];;CurFitRangeUp[[r]],3]]  ]^2]],{r,1,CurWlExtent}];

CurPeakFitBin=
Table[ Quiet[NonlinearModelFit[ Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[CurFitRangeLow[[r]];;CurFitRangeUp[[r]],{1,2}]]  ],{FitFunc[x,G,M,A,B,K,L,F]},{{G,0.01},{M,CurMaxPos[[r]]},{A,1},{B,0.0000001},{K,0.0000001},{L,0.0000001},{F,0.000001}},x,Method->LevenbergMarquardt,MaxIterations->5000,Weights->1/Re[ CurWilsonLoopsDFTJack[[JckBin]][[r]][[CurFitRangeLow[[r]];;CurFitRangeUp[[r]],3]]  ]^2]],{r,1,CurWlExtent}];


AppendTo[CurPeakFitJack,CurPeakFitBin];
AppendTo[CurBWPeakFitJack,CurBWOnlyPeakFitBin];

Print[Show[ListPlot[Table[Map[{Re[#[[1]]],Around[Re[#[[2]]],Re[#[[3]]]]}&,CurWilsonLoopsDFTJack[[JckBin]][[r]][[All,{1,2,3}]]],{r,1,4}]
,AxesLabel->{"t","Re[W]"},AxesStyle->{Gray, AbsoluteThickness[3],FontSize->20},PlotLabel->Text[Style["Real-part of the DFT of the Wilson Loop",FontSize->10]]  ],Plot[Evaluate@Table[FitFunc[x,G,M,A,B,K,L,F]/.CurPeakFitBin[[r]]["BestFitParameters"],{r,1,4}],{x,0,0.3},PlotStyle->{AbsoluteThickness[3]},PlotRange->{{ 0,0.3  },{0,CurMaxVal[[1]]}}],PlotRange->{{ 0,0.3   },{0,4}},ImageSize->400]];

,{JckBin,1,NJack}];

(*Achtung 1->NJack*)


(*Output the results from the fits for later plotting to HD *)
 
Quiet[DeleteFile[ToString@StringForm["V``B``ColSourceWilsonDFTFitParms.dat",CurNx,CurBetasChar]]];
WLDFTOutput=OpenAppend[ToString@StringForm["V``B``ColSourceWilsonDFTFitParms.dat",CurNx,CurBetasChar],FormatType->StandardForm,PageWidth->Infinity];

For[JckBin=1,JckBin<=NJack,JckBin++,
For[r=1,r<=Length[CurWilsonLoopsDFT],r++,
Write[WLDFTOutput,JckBin, " ", r, " " ,CForm[ (M/.CurPeakFitJack[[JckBin]][[r]]["BestFitParameters"]) ], " ",CForm[ (G/.CurPeakFitJack[[JckBin]][[r]]["BestFitParameters"]) ], " " ,CForm[ (A/.CurPeakFitJack[[JckBin]][[r]]["BestFitParameters"]) ], " ",CForm[ (B/.CurPeakFitJack[[JckBin]][[r]]["BestFitParameters"]) ], " ",CForm[ (K/.CurPeakFitJack[[JckBin]][[r]]["BestFitParameters"]) ], " ", CForm[ (L/.CurPeakFitJack[[JckBin]][[r]]["BestFitParameters"]) ], " ",CForm[ (F/.CurPeakFitJack[[JckBin]][[r]]["BestFitParameters"]) ], " "    ];
];
];
Close[WLDFTOutput];

(* If we find differences between the fit values of the direct exponential fit and the BW it may be related to ringing 
   artifacts introduced by the not perfect decay of the Wilson loop over time.  *)

(*Print[CurNaiveExpFit[[1]]["BestFitParameters"],CurPeakFitBin[[1]]["BestFitParameters"],CurBWOnlyPeakFitBin[[1]]["BestFitParameters"]];
Print[CurNaiveExpFit[[2]]["BestFitParameters"],CurPeakFitBin[[2]]["BestFitParameters"],CurBWOnlyPeakFitBin[[2]]["BestFitParameters"]];*)

(*Now collect the extracted values and compose the result for the potential*) 

AppendTo[ReV,{ CurBeta, Table[  {Cura*r, Mean[Table[(M/.CurPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]] , Sqrt[((NJack-1)/NJack)*Variance[Table[(M/.CurPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]]],Mean[Table[(CurPeakFitJack[[k]][[r]]["ParameterErrors"][[2]]),{k,1,NJack}]]   }       ,{r,1,CurWlExtent}]}];

AppendTo[ImV,{ CurBeta, Table[  {Cura*r, Mean[Table[(Abs[G]/.CurPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]] , Sqrt[((NJack-1)/NJack)*Variance[Table[(Abs[G]/.CurPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]]],Mean[Table[(CurPeakFitJack[[k]][[r]]["ParameterErrors"][[1]]),{k,1,NJack}]]   }       ,{r,1,CurWlExtent}]}];

AppendTo[ReVBW,{ CurBeta, Table[  {Cura*r, Mean[Table[(M/.CurBWPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]] , Sqrt[((NJack-1)/NJack)*Variance[Table[(M/.CurBWPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]]], Mean[Table[(CurBWPeakFitJack[[k]][[r]]["ParameterErrors"][[2]]),{k,1,NJack}]]    }       ,{r,1,CurWlExtent}]}];

AppendTo[ImVBW,{ CurBeta, Table[  {Cura*r, Mean[Table[(Abs[G]/.CurBWPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]] , Sqrt[((NJack-1)/NJack)*Variance[Table[(Abs[G]/.CurBWPeakFitJack[[k]][[r]]["BestFitParameters"]),{k,1,NJack}]]],Mean[Table[(CurBWPeakFitJack[[k]][[r]]["ParameterErrors"][[1]]),{k,1,NJack}]]   }       ,{r,1,CurWlExtent}]}];

]

"Plot the extracted potential values and outpout to HD.";

ListPlot[Table[Map[{#[[1]],Around[#[[2]],Sqrt[#[[3]]^2+#[[4]]^2] ]*20}&,ImV[[b]][[2]]],{b,1,1}]]

ListPlot[Table[Map[{#[[1]],Around[#[[2]],  Sqrt[#[[3]]^2+#[[4]]^2] ] }&,ReV[[b]][[2]]],{b,1,1}]]

Export["ImVColSourceMathematica.dat",Compress[ImV],"String"]

Export["ImVBWColSourceMathematica.dat",Compress[ImVBW],"String"]

Export["ReVColSourceMathematica.dat",Compress[ReV],"String"]

Export["ReVBWColSourceMathematica.dat",Compress[ReVBW],"String"]

OverallNx=32;

Quiet[DeleteFile[ToString@StringForm["V``ColSourcePot.dat",OverallNx]]];
PotOutput=OpenAppend[ToString@StringForm["V``ColSourcePot.dat",OverallNx],FormatType->StandardForm,PageWidth->Infinity];

For[b=1,b≤1,b++,
For[r=1,r<=Length[ImV[[b]][[2]]],r++,

Write[PotOutput, 3," ", Betas[[3]], " ", r, " ", CForm[ ReV[[b]][[2]][[r,2]] ], " " , CForm[  ImV[[b]][[2]][[r,3]] ], " " , CForm[  ImV[[b]][[2]][[r,4]] ], " ", CForm[ ImV[[b]][[2]][[r,2]] ], " ", CForm[ ImV[[b]][[2]][[r,3]] ], " ", CForm[ ImV[[b]][[2]][[r,4]] ]    ];

]
];
Close[PotOutput];

"Debye mass";

SigmaLat=Gamma[1/24]^2*Gamma[11/24]^2*(Sqrt[3]-1)/(48*Pi^2);

Nc=3;

Sqrt[N[(SigmaLat/(4*Pi))*2*Nc/12]]

Sqrt[N[(SigmaLat/(4*Pi))*2*Nc/16]]

Sqrt[N[(SigmaLat/(4*Pi))*2*Nc/20]]