NIntegrate[1,{x,-1,1},{y,-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)])}]

r[rs_,x_,y_]:=(Sqrt[Norm[rs]+(x^2)+(y^2)-2Part[rs,1]x-2Part[rs,2]y]);

int[rp_,rs_,x_,y_,λ_]:=((Exp[I ((2π)/λ)(r[rs,x,y]+r[rp,x,y])])/(r[rp,x,y]r[rs,x,y]))(((Part[rs,3])/(r[rs,x,y]))-((Part[rp,3])/(r[rp,x,y])));

λ=500.0×(10^(-9));
rs={0.0,0.0,-0.50};
rp={0.0,0.0,0.50};
L=(1/(1/0.50+1/0.50))
R0=0.0008;
up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-R0,R0},{y,-(Sqrt[(R0^2)-(x^2)]),(Sqrt[(R0^2)-(x^2)])}]
up=(I/(2λ))NIntegrate[int[rp,rs,ρ Cos[θ],ρ Sin[θ],λ]ρ,{ρ,0,R0},{θ,0,2π}]
R1=(Sqrt[L λ]);
R2=(Sqrt[2])(Sqrt[L λ]);
R3=(Sqrt[3])(Sqrt[L λ]);
R4=(Sqrt[4])(Sqrt[L λ]);
R5=(Sqrt[5])(Sqrt[L λ]);
(((R0/R1))^2)
up=(I/(2λ))NIntegrate[int[rp,rs,ρ Cos[θ],ρ Sin[θ],λ]ρ,{ρ,0,R1},{θ,0,2π}]
up=(I/(2λ))NIntegrate[int[rp,rs,ρ Cos[θ],ρ Sin[θ],λ]ρ,{ρ,0,R2},{θ,0,2π}]
up=(I/(2λ))NIntegrate[int[rp,rs,ρ Cos[θ],ρ Sin[θ],λ]ρ,{ρ,0,R3},{θ,0,2π}]
up=(I/(2λ))NIntegrate[int[rp,rs,ρ Cos[θ],ρ Sin[θ],λ]ρ,{ρ,0,R4},{θ,0,2π}]
up=(I/(2λ))NIntegrate[int[rp,rs,ρ Cos[θ],ρ Sin[θ],λ]ρ,{ρ,0,R5},{θ,0,2π}]

up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,,},{y,-0.0005,0.0005}]

rs={0.0,0.0,-50.0};
n=400;
data={};
ypmin=-0.010;
ypmax=0.010;
λ=500.0×(10^(-9));
Do[
yp=ypmin+(ypmax-ypmin)i/n;
rp={0.0,yp,50.0};
up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-0.0005,0.0005},{y,-0.0005,0.0005}];
data=Append[data,{yp,((Abs[up])^2)}];
If[Mod[i,100]==0,Print[i]];
,{i,0,n}]

ListPlot[data,Joined->True,PlotRange->All]

rs={0.0,0.0,-50.0};
n=400;
data={};
ypmin=-0.010;
ypmax=0.010;
s=0.002;
λ=500.0×(10^(-9));
Do[
yp=ypmin+(ypmax-ypmin)i/n;
rp={0.0,yp,50.0};
up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-0.0005,0.0005},{y,-s-0.0005,-s+0.0005}]+(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-0.0005,0.0005},{y,s-0.0005,s+0.0005}];
data=Append[data,{yp,((Abs[up])^2)}];
If[Mod[i,100]==0,Print[i]];
,{i,0,n}]

ListPlot[data,Joined->True,PlotRange->All]

rs={0.0,0.0,-0.50};
n=400;
data={};
ypmin=-0.0050;
ypmax=0.0050;
λ=500.0×(10^(-9));
Do[
yp=ypmin+(ypmax-ypmin)i/n;
rp={0.0,yp,0.50};
up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-0.0005,0.0005},{y,-0.0005,0.0005}];
data=Append[data,{yp,((Abs[up])^2)}];
If[Mod[i,100]==0,Print[i]];
,{i,0,n}]

ListPlot[data,Joined->True,PlotRange->All]

rs={0.0,0.0,-0.50};
n=400;
data={};
data2={};
ypmin=-0.0050;
ypmax=0.0050;
λ=500.0×(10^(-9));
λ2=λ/2;
Do[
yp=ypmin+(ypmax-ypmin)i/n;
rp={0.0,yp,0.50};
up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-0.0005,0.0005},{y,-(Sqrt[(0.0005^2)-(x^2)]),(Sqrt[(0.0005^2)-(x^2)])},AccuracyGoal->8];
up2=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-λ2,λ2},{y,-(Sqrt[((λ2)^2)-(x^2)]),(Sqrt[((λ2)^2)-(x^2)])},AccuracyGoal->8];
data=Append[data,{yp,((Abs[up])^2)}];
data2=Append[data2,{yp,((Abs[up2])^2)}];
If[Mod[i,100]==0,Print[i]];
,{i,0,n}]

ListPlot[data,Joined->True,PlotRange->All]

ListPlot[data2,Joined->True,PlotRange->All]

rs={0.0,0.0,-50.0};
n=30;
data={};
xpmin=-0.008;
xpmax=0.008;
ypmin=-0.008;
ypmax=0.008;
λ=500.0×(10^(-9));
Do[
xp=xpmin+(xpmax-xpmin)j/n;
Do[
yp=ypmin+(ypmax-ypmin)i/n;
rp={xp,yp,50.0};
up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-0.0005,0.0005},{y,-0.0005,0.0005},AccuracyGoal->8];
data=Append[data,{xp,yp,((Abs[up])^2)}];
,{i,0,n}];
If[Mod[j,5]==0,Print[j]];
,{j,0,n}]

ListPlot3D[data,PlotRange->All]

ListDensityPlot[data,InterpolationOrder->3]

rs={0.0,0.0,-50.0};
n=30;
data={};
xpmin=-0.008;
xpmax=0.008;
ypmin=-0.008;
ypmax=0.008;
λ=500.0×(10^(-9));
Do[
xp=xpmin+(xpmax-xpmin)j/n;
Do[
yp=ypmin+(ypmax-ypmin)i/n;
rp={xp,yp,50.0};
up=(I/(2λ))NIntegrate[int[rp,rs,x,y,λ],{x,-0.0005,0.0005},{y,-(Sqrt[(0.0005^2)-(x^2)]),(Sqrt[(0.0005^2)-(x^2)])},AccuracyGoal->8];
data=Append[data,{xp,yp,((Abs[up])^2)}];
,{i,0,n}];
If[Mod[j,5]==0,Print[j]];
,{j,0,n}]

ListPlot3D[data,PlotRange->All]
ListDensityPlot[data,InterpolationOrder->3]