NotebookEvaluate[NotebookDirectory[]<>"/split_step.nb"];

defaultscale=8;
HarmonicTrap1D[];
sizes={1024};
k=0;
contactInteractionFactor=4 Pi k;
initialize[];
gridReport

res={};
do[ts_]:=Module[{},
timeStep=ts;
initialize[];reset;
evolve["ite",100 ,10];
AppendTo[res,Join[{timeStep},energies,{psi}]];
];
doall:={
do[0.0001];
do[0.000316];
do[0.001];
do[0.00316];
do[0.01];
do[0.0316];
do[0.1];
};

testGSenergy:=Module[{refE,l,f},
refE=res[[1,2]];
l=Sort@Map[#-{0,refE}&,res[[All,{1,2}]]];
Print[ListLogLogPlot[l]];
l=Log[Drop[l,1]];
f=Fit[l,{1,x},x];
Print[f];
];
testVirial:=Module[{l,f},
l=Sort@res[[All,{1,6}]];
l=Abs[l];
Print[ListLogLogPlot[l]];
l=Log[Drop[l,1]];
f=Fit[l,{1,x},x];
Print[f];
];
testWFn:=Module[{psi0,l,f},
psi0=res[[1,-1]];
l=Sort@Table[{res[[i,1]],NIntegrate[(Abs[res[[i,-1]][x]]^2-Abs[psi0[x]]^2)^2,{x,-8,8}]},{i,Length[res]}];
Print[ListLogLogPlot[l]];
l=Log[Drop[l,1]];
f=Fit[l,{1,x},x];
Print[f];
];

res={};
algorithm="firstOrder";
k=0;
contactInteractionFactor=4 Pi k;
initialize[];
reset;
doall;

testGSenergy

testVirial

testWFn

res={};
algorithm="secondOrder";
k=0;
contactInteractionFactor=4 Pi k;
initialize[];
reset;
doall;

testGSenergy

testVirial

testWFn

res={};
algorithm="firstOrder";
k=-0.2;
contactInteractionFactor=4 Pi k;
initialize[];
reset;
doall;

testGSenergy

testVirial

testWFn

res={};
algorithm="secondOrder";
k=-0.2;
contactInteractionFactor=4 Pi k;
initialize[];
reset;
doall;

testGSenergy

testVirial

testWFn