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

HarmonicTrap1D[];
sizes={4096};
timeStep=0.00025;
k=-0.5;
contactInteractionFactor=4 Pi k;
max=50ho;
min=-50ho;
initialize[];
gridReport

calcAllTable
p1=plotProjections;

AbsoluteTiming[evolve["ite",10,10]]

calcAllTable
p2=plotProjections;
Show[p1,p2,PlotRange->{{-3,3},All}]

Turn off harmonic potential

plot:=AppendTo[lplots,{time,plotProjections,interp[wf]}];
newPotential[(0&)];
reset;
AbsoluteTiming[evolve["ite",10,10,10]]

{plotM["meanX"],plotM["sigmaX"]}

wfsoliton=wf;
psi=interp[wfsoliton]
wfsoliton2=Map[psi[First[#]]&,mesh[],{dim}];
Norm[wfsoliton-wfsoliton2]

x0=-10;
wfshifted=Map[Piecewise[{{psi[First[#]-x0],Abs[First[#]-x0]<10}},0]&,mesh[],{dim}];

Scattering on a Gaussian potential bump

ss=0.1;
strength=4;
gaussian[r_List]:=strength(Pi^(-1/4)ss^(-1/2)Exp[-r[[1]]^2/(2ss^2)]);
Plot[gaussian[{x}],{x,-1,1},PlotRange->All]
newPotential[gaussian];

p=3; (* momentum *)
wf=wfshifted*Exp[I p meshX];
ListPlot[projectionX,PlotRange->All,Joined->True]

reset;
AbsoluteTiming[evolve["rte",15,50,10]]

plotM["meanX"]

GraphicsGrid[Table[{lplots[[i,2]]},{i,Length[lplots]}],ImageSize->3 72]

fitcosh[psi_,min_,max_]:=Module[{tab,fnc,sol,a,b,c},
tab=Table[{x,Abs[psi[x]]^2},{x,min,max,(max-min)/100}];
fnc=(a/Cosh[(x-c)/b])^2;
sol=FindFit[tab,fnc,{{a,1},{b,1},{c,(max+min)/2}},x];
fnc=fnc/.sol;
Print[fnc];
fnc
];

fnc=fitcosh[lplots[[-1,3]],-30,-24];
pfit=Plot[fnc,{x,-29,-25},PlotStyle->Red];
Show[lplots[[-1,2]],pfit,PlotRange->{{-30,-23},All}]

fnc=fitcosh[lplots[[-1,3]],28,35];
pfit=Plot[fnc,{x,30,33},PlotStyle->Red];
Show[lplots[[-1,2]],pfit,PlotRange->{{28,35},All}]