n=1.59/1.33;  (*Refraction Index Ratio*)
r=0.5;                (* Particle Radius*)
σ=0.1;                 (* StandardDeviation*)


Deltab[x_,r_]:=(x/r)√(1-((((n x)/r))^2))-((n x)/r)√(1-(((x/r))^2)) (*DeltaBubble*)
Int[x_,σ_]:=((E^(-((((x-r)/σ))^2)))+(E^(-((((x+r)/σ))^2)))) (*Laser Intensity Profile*)


Plot[{Deltab[x,r]},{x,-r,r},PlotRange->All]   (*Plot DeltaBubble*)
Plot[{Int[x,σ]},{x,-4,4},PlotRange->All] (*Plot Laser Intensity Profile*)
Plot[{Deltab[x,r]*Int[x,σ]},{x,-r,r},PlotRange->All] (*Plot ProductFunction [when integrated it gives the force]]*)
 
Fb[x0_,σ_,a_]:=-Re[(∑_(i=1)^10000)Deltab[(-a+(i 2 a)/10000),a]Int[(-a+(i 2 a)/10000)-x0,σ]((2a)/10000)]
(*Numerical Integral -Computation of the Force-*)

(*Data export for Force and Potential*)
sigma=0.1;radio =0.50; delta = 0.01;long=16; 
(*delta=step,radio=radius,long=range to obtain data, *)
Tabfuerza= Table[{pos1, Fb[pos1, sigma, radio]}, {pos1, -(long/(2)), (long/(2)), delta}];
(*Fueza=Force*)
gfuerza = ListPlot[Tabfuerza,PlotRange->All];
countmon = Count[Flatten[Tabfuerza], _Real]/2;
Potencial= Table[1, {i, 1, countmon}];
(*Potencial=Potencia*)
For[i = 2, i < countmon +1, i++, Potencial[[i]] = -Tabfuerza[[i, 2]]*delta + Potencial[[i - 1]]]
Tabpotencial = Table[{Tabfuerza[[i, 1]], Potencial[[i]]}, {i, 1, countmon}];
gpotencial = ListPlot[Tabpotencial,PlotRange->All];
Export["Force.dat",Tabfuerza];
Export["Potential.dat",Tabpotencial];
GraphicsGrid[{{gfuerza, gpotencial}}]