(* ::Package:: *)

(* Choose p for type of potential *)
(*p = 1 for square potential barrier *)
(*p = 2 for square LJ *)
(*p = 3 for infinite well *)
(*p = 4 for 1-D harmonic oscillator *)
(*p = 5 for 1-D harmonic oscillator *)
p=5;


(*square potential barrier*)
If[p==1,
(************ parameters *******************)
DeltaEorig=0.1;
tol=10^(-4); iota=10^(-4);b=1000;
NN=5000;(*default 10000*)
Eini=0.5;k=1.0;nlast=5;
(************ parameters *******************)


(** (length/2) = xlastorig; *)
(** (w/2)+ L + skd) + (xlastorig-xlast) = (length/2)*)
(** skd = xlast - L = gap between xlastorig and xlast **)
(*length=xlastorig-xiniorig;*)
h=100;H=10^5;
(* Begin of sum of the coefficients in skd, L, w must be less but close to 1.0 *)
skd=0.1;(*by default set coefficient to (length/2) to 0.01 *)
L=1-0.1;
w=0.2;
xlastorig=L + (w/2) + skd;
xiniorig=-xlastorig;
(* End of sum of the coefficients in skd, L, w must be less but close to 1.0 *)

Potential[x_]:=Piecewise[{
{h,0<=Abs[x]<=w/2},{0,w/2<Abs[x]<(w/2)+L},{H,(w/2)+L<=Abs[x]<xlastorig}
}];

xlast = xlastorig;
xini = -xlast;

(* Find xmin, x where V is minimum *)
(* Either define xMin manually or use Minimize[] *)
xmin=Minimize[Potential[x],x>0,x]; xMin=x/.xmin[[2]][[1]];
(*xMin=0.3;*)
iMin=IntegerPart[xMin/Deltax];


Deltax=xlastorig/NN;
iiniorig=IntegerPart[xiniorig/Deltax];
ilastorig=IntegerPart[xlastorig/Deltax];

For[i=iiniorig,i<=ilastorig,i++,
xdum=i*Deltax;
V[i]=Potential[xdum]
];(*end for i*)

Print[Plot[Potential[x],{x,xiniorig,xlastorig}]];
];(*end if p=1*)


(* LJ potential *)
If[p==2,
(************ parameters *******************)
iota=10^(-4);
xiniorig=0.1;
xlastorig=5.0;
DeltaEorig=0.1;
tol=10^(-5); b=50;
NN=5000;(*default 10000*)
Eini=-2.0;k=0.9;nlast=5;
(************ parameters *******************)

sigma=1;epsilon=10;
Potential[xdum_]:=4epsilon*((sigma/xdum)^12-(sigma/xdum)^6);
V[0]=Potential[iota];

(* Find xmin, x where V is minimum *)
(* Either define xMin manually or use Minimize[] *)
xmin=Minimize[Potential[x],x>0,x]; xMin=x/.xmin[[2]][[1]];
(*xMin=0;*)
iMin=IntegerPart[xMin/Deltax];


(* Tune the values of xini and xlast appropririately *)
xini=xiniorig+0.5;xlast=xlastorig-0.5;
(* End of tuning the values of xini and xlast appropririately *)


Print[Plot[Potential[x],{x,xiniorig,xlastorig}]]
];(*end if p=2*)



(*  infinite square well *)
If[p==3,
(************ parameters *******************)
DeltaEorig=0.1;
tol=10^(-4); iota=10^(-4);b=1000;
NN=7000;(*default 10000*)
Eini=1.2;k=1.5;nlast=5;
(************ parameters *******************)


(** (length/2) = xlastorig; *)
(** (w/2)+ L + skd) + (xlastorig-xlast) = (length/2)*)
(** skd = xlast - L = gap between xlastorig and xlast **)
(*length=xlastorig-xiniorig;*)
h=0;H=10^3;
(* Begin of sum of the coefficients in skd, L, w must be less but close to 1.0 *)
skd=0.1;(*by default set coefficient to (length/2) to 0.01 *)
L=1;
w=0.0;
xlastorig=L + (w/2) + skd;
xiniorig=-xlastorig;
(* End of sum of the coefficients in skd, L, w must be less but close to 1.0 *)

Potential[x_]:=Piecewise[{
{h,0<=Abs[x]<=w/2},{0,w/2<Abs[x]<(w/2)+L},{H,(w/2)+L<=Abs[x]<xlastorig}
}];

xlast = L+(w/2);
xini = -xlast;

(* Find xmin, x where V is minimum *)
(* Either define xMin manually or use Minimize[] *)
(*xmin=Minimize[Potential[x],x>0,x]; xMin=x/.xmin[[2]][[1]];*)
xMin=0.0;
iMin=IntegerPart[xMin/Deltax];


Deltax=xlastorig/NN;
iiniorig=IntegerPart[xiniorig/Deltax];
ilastorig=IntegerPart[xlastorig/Deltax];

For[i=iiniorig,i<=ilastorig,i++,
xdum=i*Deltax;
V[i]=Potential[xdum]
];(*end for i*)

Print[Plot[Potential[x],{x,xiniorig,xlastorig}]];
];(*end if p=1*)


(*  1-D harmonic oscillator *)
If[p==4,
(************ parameters *******************)
iota=10^(-4);
xiniorig=-10;
xlastorig=-xiniorig;
DeltaEorig=0.1;
tol=10^(-5); b=1000;
NN=10000;(*default 10000*)
Eini=-0.1; k=0.1;nlast=5;
(************ parameters *******************)

K=1;
Potential[xdum_]:=K*(xdum^2)/2;
(*V[0]=Potential[iota];*)

(* Find xmin, x where V is minimum *)
(* Either define xMin manually or use Minimize[] *)
xmin=Minimize[Potential[x],x]; xMin=x/.xmin[[2]][[1]];
(*xMin=0;*)
iMin=IntegerPart[xMin/Deltax];

(* Tune the values of xini and xlast appropririately *)
xini=xiniorig;xlast=xlastorig;
(* End of tuning the values of xini and xlast appropririately *)


Print[Plot[Potential[x],{x,xiniorig,xlastorig}]]
];(*end if p=4*)


(*  Coulomb potentail  *)
If[p==5,
(************ parameters *******************)
iota=10^(-4);
xiniorig=-20;
xlastorig=20;
DeltaEorig=0.1;
tol=10^(-5); b=1000;
NN=10000;(*default 10000*)
Eini=-15; k=0.1;nlast=5;
(************ parameters *******************)

K=1;
Potential[xdum_]:=-1/xdum^2;
V[0]=Potential[10^2*iota];

(* Find xmin, x where V is minimum *)
(* Either define xMin manually or use Minimize[] *)
(*xmin=Minimize[Potential[x],x]; xMin=x/.xmin[[2]][[1]];*)
(*xMin=(xiniorig+xlastorig)/2;*)
xMin=1;
iMin=IntegerPart[xMin/Deltax];

(* Tune the values of xini and xlast appropririately *)
xini=xiniorig;xlast=xlastorig;
(* End of tuning the values of xini and xlast appropririately *)


Print[Plot[Potential[x],{x,xiniorig,xlastorig}]]
];(*end if p=5*)
