# I want to visualize this 3D like DensityPlot3D (About Diffraction of light)

Here is a code modeling the diffraction of light. I made this because I want simulate the diffraction of light and visualize the light's density. It's Fresnel diffraction, and I want to see the Density of light according to distance about slit and screen. 'z' meens distance about slit and screen. 'lamda' meens wavelength of light. This is a slit that has hole and disk. and 'a' meens diameter of hole and 'z' meens diameter of Disk. This is simulation about hole aed Disk, and it can make effect Disk and hole, because of Interference of light

z = 16;
a = 2/10^3;
m = 6/10^3;
λ = 532/10^9;
k = (2*Pi)/λ;
n = (k*q*m)/z;
v = (k*q*a)/z;
q = Sqrt[x^2 + y^2];
Subscript[U, 1] =
Sum[(a/q)^(4*i + 1)*BesselJ[4*i + 1, v], {i, 0, 3}] -
Sum[(a/q)^(4*i + 3)*BesselJ[4*i + 3, v], {i, 0, 3}];
Subscript[U, 2] =
Sum[(a/q)^(4*i + 2)*BesselJ[4*i + 2, v], {i, 0, 3}] -
Sum[(a/q)^(4*i + 4)*BesselJ[4*i + 4, v], {i, 0, 3}];
Subscript[V, 0] =
Sum[(q/m)^(4*i)*BesselJ[4*i, n], {i, 0, 3}] -
Sum[(q/m)^(4*i + 2)*BesselJ[4*i + 2, n], {i, 0, 3}];
Subscript[V, 1] =
Sum[(q/m)^(4*i + 1)*BesselJ[4*i + 1, n], {i, 0, 3}] -
Sum[(q/m)^(4*i + 3)*BesselJ[4*i + 3, n], {i, 0, 3}];


When I evaluate the following DensityPlot expression,

DensityPlot[
Subscript[U, 1]^2 + Subscript[U, 2]^2 +
Subscript[V, 0]^2 + Subscript[V, 1]^2 + 1 +
2*((Subscript[U, 1]*Subscript[V, 0] +
Subscript[U, 2]*Subscript[V, 1])*Sin[(k*(m^2 - a^2))/(2*z)] +
(Subscript[U, 2]*Subscript[V, 0] -
Subscript[U, 1]*Subscript[V, 1])*Cos[(k*(m^2 - a^2))/(2*z)] +
Subscript[U, 1]*Sin[(k*(q^2 + a^2))/(2*z)] -
Subscript[U, 2]*Cos[(k*(q^2 + a^2))/(2*z)] -
Subscript[V, 0]*Cos[(k*(q^2 + m^2))/(2*z)] -
Subscript[V, 1]*Sin[(k*(q^2 + m^2))/(2*z)]),
{x, -1./10^2, 1./10^2}, {y, -1./10^2, 1./10^2},
PlotRange -> {0, 6.}, PlotPoints -> 120, Mesh -> False,
Frame -> False, AspectRatio -> 1]


I get this image My problem/questions:

1. How can I visuailze this image using like that? 2. I want to see the 3D density plot about variable z. In this code, I fixed z = 16, but I want make this variable by using lever.
• Can you maybe mention where you got these formulae from? – J. M. is in limbo Nov 3 '17 at 17:35
• @Relativity I would like to be in context. It seems to me the Fresnel diffraction calculation in a screen at some distance, and you want to simulate how it evolves as you move the screen along the axis, isn’t it? Maybe by a zone plate? – José Antonio Díaz Navas Nov 3 '17 at 23:01
• As Jose and J.M. have said knowing context would help determine if there is a better way to visualize what you want to see. A possible alternative way is to do a DensityPlot with q and z on your axes instead of a DensityPlot3D with x,y and z. – qbit Nov 3 '17 at 23:17
• Yes, it's Fresnel diffraction, and I want to see the Density of light according to distance about slit and screen. 'z' meens distance about slit and screen. 'lamda' meens wavelength of light. This is a slit that has hole and disk. – Relativity Nov 3 '17 at 23:55
• and 'a' meens diameter of hole, 'm' meens diameter of Disk. – Relativity Nov 3 '17 at 23:57

This is my try. First, I must say that your plotting approach is a bit slow, so I have speed up the calculations a bit by approximating your $U1,U2,V0$ and $V1$ . My result is five times faster, using the same number of points.

Second, the use of Manipulate is sloppy as the time to generate the plot is about 6 sec for each value of $z$. Therefore, I have been in favor to see an animation of how the diffraction pattern evolves as $z$ increases. Thus, I have generated the images, and animated them. 3D plots are generated a bot more slowly, but if you want to see them just change DensityPlot by Plot3Dand modify the code and the options accordingly.

 Clear[z, a, m, \[Lambda], k, n, v, q, U1, U2, V0, V1, interp1, interp2, interp4, interp3, f, g, img, imgs];

a = 2/10^3;
m = 6/10^3;
\[Lambda] = 532/10^9;
k = (2*Pi)/\[Lambda];

imgs = {};

Do[
n = (k*q*m)/z;
v = (k*q*a)/z;

U1 = Sum[(a/q)^(4*i + 1)*BesselJ[4*i + 1, v], {i, 0, 3}] -Sum[(a/q)^(4*i + 3)*BesselJ[4*i + 3, v], {i, 0, 3}];
U2 = Sum[(a/q)^(4*i + 2)*BesselJ[4*i + 2, v], {i, 0, 3}] - Sum[(a/q)^(4*i + 4)*BesselJ[4*i + 4, v], {i, 0, 3}];
V0 = Sum[(q/m)^(4*i)*BesselJ[4*i, n], {i, 0, 3}] - Sum[(q/m)^(4*i + 2)*BesselJ[4*i + 2, n], {i, 0, 3}];
V1 = Sum[(q/m)^(4*i + 1)*BesselJ[4*i + 1, n], {i, 0, 3}] - Sum[(q/m)^(4*i + 3)*BesselJ[4*i + 3, n], {i, 0, 3}];

interp1 = Quiet@FunctionInterpolation[U1, {q, -0.03, 0.03}, AccuracyGoal -> 5];
interp2 = Quiet@FunctionInterpolation[U2, {q, -0.03, 0.03}, AccuracyGoal -> 5];
interp3 = Quiet@FunctionInterpolation[V0, {q, -0.03, 0.03}, AccuracyGoal -> 5];
interp4 = Quiet@FunctionInterpolation[V1, {q, -0.03, 0.03}, AccuracyGoal -> 5];

f = Subscript[U, 1]^2 + Subscript[U, 2]^2 + Subscript[V, 0]^2 + Subscript[V, 1]^2 + 1 + 2*((Subscript[U, 1]*Subscript[V, 0] + Subscript[U, 2]*Subscript[V, 1])*Sin[(k*(m^2 - a^2))/(2*z)] + (Subscript[U, 2]*Subscript[V, 0] - Subscript[U, 1]*Subscript[V, 1])*Cos[(k*(m^2 - a^2))/(2*z)] + Subscript[U, 1]*Sin[(k*(q^2 + a^2))/(2*z)] - Subscript[U, 2]*Cos[(k*(q^2 + a^2))/(2*z)]-Subscript[V, 0]*Cos[(k*(q^2 + m^2))/(2*z)] - Subscript[V, 1]*Sin[(k*(q^2 + m^2))/(2*z)])
//. {Subscript[U,1] -> interp1[q], Subscript[U, 2] -> interp2[q], Subscript[V, 0] -> interp3[q], Subscript[V, 1] -> interp4[q], q -> Sqrt[x^2 + y^2] } // N // FullSimplify ;

img = Image@DensityPlot[Evaluate@f, {x, -1/10^2, 1/10^2}, {y, -1/10^2, 1/10^2},PlotRange -> {{-.008, 0.008}, {-.008, 0.008}, {0, 6.}}, PlotPoints -> 60, Mesh -> False, Frame -> False, RegionFunction -> (#1^2 + #2^2 <= 0.006^2 &),PlotLabel -> Style["z = " <> ToString[z] , 24]];

(*img=Image@Plot3D[Evaluate@(f),{x,-1./10^2,1./10^2},{y,-1./10^2,1./10^2},PlotRange -> {{-.008, 0.008}, {-.008, 0.008}, All},PlotPoints -> 60, Mesh -> False, Boxed - >False,Axes -> False, BoxRatios -> {1, 1,2}, ColorFunction -> "Rainbow", RegionFunction -> (#1^2+#2^2<=0.006^2&), PlotLabel -> Style["z = " <> ToString[z], 24]];*)

imgs = AppendTo[imgs, img],
{z, 8, 35}
];
ListAnimate[imgs, 2] Animated gifs:  As a minor point $V0$ and $V1$ functions oscillate fast and with very high values near 0.006 and beyond, so I have restricted the size a bit less that you have specified. If not, in 3D plots there will be higher values hindering the visibility of the center of the pattern. However, you can modify the code as you wish.

Finally, my experience in the field says to me that it is better to work numerically with diffraction patterns via Fourier transforms. As you probably know, Fresnel diffraction can be done via Fourier transform of the field in the aperture multiplied by a quadratic phase factor in the aperture.

• Thanks to you^^ – Relativity Nov 5 '17 at 10:12
• Sorry, I copy your code and push shift+enter, but mathematica says error. – Relativity Nov 5 '17 at 10:59
• Clear[z, a, m, [Lambda], k, n, v, q, U1, U2, V0, V1, interp1, interp2, interp4, interp3, f, g, img, imgs]; a = 2/10^3; m = 6/10^3; [Lambda] = 532/10^9; k = (2*Pi)/[Lambda]; imgs = {}; Do[n = (kqm)/z; v = (kqa)/z; U1 = Sum[(a/q)^(4*i + 1)*BesselJ[4*i + 1, v], {i, 0, 3}] - Sum[(a/q)^(4*i + 3)*BesselJ[4*i + 3, v], {i, 0, 3}]; U2 = Sum[(a/q)^(4*i + 2)*BesselJ[4*i + 2, v], {i, 0, 3}] - Sum[(a/q)^(4*i + 4)*BesselJ[4*i + 4, v], {i, 0, 3}]; – Relativity Nov 5 '17 at 11:00
• V0 = Sum[(q/m)^(4*i)*BesselJ[4*i, n], {i, 0, 3}] - Sum[(q/m)^(4*i + 2)*BesselJ[4*i + 2, n], {i, 0, 3}]; V1 = Sum[(q/m)^(4*i + 1)*BesselJ[4*i + 1, n], {i, 0, 3}] - Sum[(q/m)^(4*i + 3)*BesselJ[4*i + 3, n], {i, 0, 3}]; interp1 = Quiet@FunctionInterpolation[U1, {q, -0.03, 0.03}, AccuracyGoal -> 5]; interp2 = Quiet@FunctionInterpolation[U2, {q, -0.03, 0.03}, AccuracyGoal -> 5]; interp3 = Quiet@FunctionInterpolation[V0, {q, -0.03, 0.03}, AccuracyGoal -> 5]; interp4 = Quiet@FunctionInterpolation[V1, {q, -0.03, 0.03}, AccuracyGoal -> 5]; – Relativity Nov 5 '17 at 11:01
• f = Subscript[U, 1]^2 + Subscript[U, 2]^2 + Subscript[V, 0]^2 + Subscript[V, 1]^2 + 1 + 2*((Subscript[U, 1]*Subscript[V, 0] + Subscript[U, 2]*Subscript[V, 1])* Sin[(k*(m^2 - a^2))/(2*z)] + (Subscript[U, 2]* Subscript[V, 0] - Subscript[U, 1]*Subscript[V, 1])* Cos[(k*(m^2 - a^2))/(2*z)] + Subscript[U, 1]*Sin[(k*(q^2 + a^2))/(2*z)] - Subscript[U, 2]*Cos[(k*(q^2 + a^2))/(2*z)] - – Relativity Nov 5 '17 at 11:03

Not sure exactly what plot you want, but can't you just use DensityPlot3D with what you have? First comment out or delete your z = 16 in your first code set, then

    DensityPlot3D[
Subscript[U, 1]^2 + Subscript[U, 2]^2 + Subscript[V, 0]^2 +
Subscript[V, 1]^2 + 1 +
2*((Subscript[U, 1]*Subscript[V, 0] +
Subscript[U, 2]*Subscript[V, 1])*
Sin[(k*(m^2 - a^2))/(2*z)] + (Subscript[U, 2]*Subscript[V, 0] -
Subscript[U, 1]*Subscript[V, 1])*Cos[(k*(m^2 - a^2))/(2*z)] +
Subscript[U, 1]*Sin[(k*(q^2 + a^2))/(2*z)] -
Subscript[U, 2]*Cos[(k*(q^2 + a^2))/(2*z)] -
Subscript[V, 0]*Cos[(k*(q^2 + m^2))/(2*z)] -
Subscript[V, 1]*Sin[(k*(q^2 + m^2))/(2*z)]), {x, -1./10^2,
1./10^2}, {y, -1./10^2, 1./10^2}, {z, 0, 10},
PlotRange -> {0,
6.}(*,PlotPoints\[Rule]120,Mesh\[Rule]False,Frame\[Rule]False,\
AspectRatio\[Rule]1*)] • Thank you! but I don't wanna that graph... Sorry. – Relativity Nov 4 '17 at 7:18
• I'll edit my question... – Relativity Nov 4 '17 at 7:18