(* Random walk of a collection of jlast walkers in 1D *)

(*Initialisation*)
nlast=100;jlast=1000;(*total number of walkers*)
L=50;L1=-L;L2=L;
xc=IntegerPart[(L1+L2)/2];

For[j=1,j≤jlast,j++,
x[j]=xc;
];(*end of j*)
avex[0]=Sum[x[j],{j,1,jlast}]/jlast;
avexsq[0]=Sum[x[j]^2,{j,1,jlast}]/jlast;
tnx[0]=Table[{x[j],0},{j,1,jlast}];
tx[0]=Table[x[j],{j,1,jlast}];
(*End of Initialisation*)

For[n=1,n≤nlast,n++,
For[j=1,j≤jlast,j++,
ran=Random[Real];
If[ran≤0.5,x[j]=x[j]-1
,(*else*)
x[j]=x[j]+1
](*end if *)
];(*end for j*)

tnx[n]=Table[{x[j],0},{j,1,jlast}];
tx[n]=Table[x[j],{j,1,jlast}];

(*evaluate the average of x[j] and xsq[j] at fixed n*)
avex[n]=Sum[x[j],{j,1,jlast}]/jlast;
avexsq[n]=Sum[x[j]^2,{j,1,jlast}]/jlast;
(*evaluate the average of x[j] and xsq[j] at fixed n*)
];(*end for n *)
tnavex=Table[{n,avex[n]},{n,0,nlast}];
tnavexsq=Table[{n,avexsq[n]},{n,0,nlast}];
lptnavexsq=ListPlot[tnavexsq,PlotLabel->"avexsq, vs n"];

yth[n_]:=2*Dx*n;Dx=0.5;
tny=Table[{n,yth[n]},{n,0,nlast}];
lptny=ListPlot[tny,PlotLabel->"yth vs n",PlotStyle->Red];
Show[{lptny,lptnavexsq}];
(*Manipulate[
ListPlot[tnx[n],PlotRange->{{L1,L2},{-1,1}}],
{n,0,nlast,1}]*)
Manipulate[Show[Histogram[tx[n],50]],{n,0,nlast,1}];

Manipulate[
ListPlot[N[BinCounts[tx[n],{L1,L2,2}]/jlast],PlotRange->{0,1},Filling->Axis]
,{n,0,nlast,1}]

BinCounts[tx[4],{L1,L2,1.0}]
BinCounts[tx[4],{L1,L2,1.0}][[55]]

(*This is an example of sorting the walkers into bins of width bw. The horizontal axis is expressed in values of xbin. The vertical axis is the normalised number of walker in bin xbin *)
totalbin=50;(**Total number of bins you decide to use*)
bw=(L2-L1)/totalbin;
xbin[j_]:=N[L1+j*bw];(*1 <= j ≤ totalbin)*)
NormalisedBC[j_,n_]:=BinCounts[tx[n],{L1,L2,bw}][[j]]/jlast;(* 1<=j≤totalbin *)
(*aveNormalisedBC[j_,n_]:=(NormalisedBC[j-1,n]+NormalisedBC[j+1,n]+NormalisedBC[j,n])/3.0;*)
(* 2<=j≤totalbin-1 *)
aveNormalisedBC[j_,n_]:=N[NormalisedBC[j,n]];
txbinnwinbinj[n_]:=Table[{xbin[j],aveNormalisedBC[j,n]},{j,2,totalbin-1}];

Manipulate[ListPlot[txbinnwinbinj[n],Joined->True,Filling->Axis,Mesh->All,PlotRange->{{L1,L2},{0,1}}],{n,0,nlast,1}]

(*Use FindFit to find the parameters A (normalisation constant), and width (=sigma) of the curve*)
lptxbinnwinbinj=ListPlot[txbinnwinbinj[100],Mesh->All,PlotStyle->Blue]
sol=FindFit[txbinnwinbinj[100], A*Exp[-(x -(xc+bw/2))^2/(2wid^2)], 
    {A, wid}, x];
Print[sol];

(* Plot function *)
A1=A/.sol;wid1->wid/.sol;
y[x1_]:=A1*Exp[-(x1 -(xc+bw/2))^2/(2wid1^2)]/.{A1->A/.sol,wid1->wid/.sol};
Print["y[x]= ",y[x]];
txy=Table[{xbin[ind],y[xbin[ind]]},{ind,1,totalbin}];
lptxy=ListPlot[txy,Joined->True,Mesh->All,PlotStyle->Red]
Show[{lptxbinnwinbinj,lptxy}]
(*ListPlot[{Table[{k,y[x]},{x,L1,L2,bw}],txbinnwinbinj[100]},Joined->True,Mesh->Automatic,PlotStyle->{Red,Black},PlotRange->All]*)

txbinnwinbinj[100]

ListPlot[txbinnwinbinj[100]]

sol=FindFit[txbinnwinbinj[100], A*Exp[-(x -(xc+bw/2))^2/(2wid^2)], 
    {A, wid}, x];
Print[sol];

sol

A1=A/.sol

wid1=wid/.sol