PHYSICAL OPTICS
Mathematica Documentation
Mathematica Source Code

(Hygncirc.ma)

Huygens' Principle allows one to figure out a given wavefront at any time in the future if the present position and shape of the wavefront are known. This is accomplished by assuming that every point in a given wavefront will serve as a source for secondary waves that radiate outward. At any given time, the shape and position of the new wavefront will be the tangents to all the secondary wavelets. This animation allows the instructor to show tothe students the propegation of a wavefront in space form a point source. This time interval between animation cells is a constant, thus resulting in a constant radius for the secondary waves generated from the secondary wavelets.

Exercises:
1. The teacher, after giving the codes to the students, can ask them to modify the program to generate the propegation of a plane wavefront using Huygens' postulates.
2. The present program does not take into effect the intensity of the secondary sources in the waves. How will the intensity of the wave front change for a point source and a plane wavefront?

Clear all variable names.
Clear["Global`*"]

Create a loop that produces wavelets that propegate outwards according to Huygens' Principle.
Do[{
arcs = Table[{ Circle[{Sqrt[r^2 - (r*Sin[Pi/(2*r)*y])^2], r*Sin[Pi/(2*r)*y]},1, {ArcTan[r*Sin[Pi/(2*r)*y]/ (Sqrt[(r+.000001)^2 - (r*Sin[Pi/(2*r)*y])^2])]-Pi/2, ArcTan[r*Sin[Pi/(2*r)*y]/ (Sqrt[(r+.000001)^2 - (r*Sin[Pi/(2*r)*y])^2])]+Pi/2}], Circle[-{Sqrt[r^2 - (r*Sin[Pi/(2*r)*y])^2], r*Sin[Pi/(2*r)*y]},1, {ArcTan[r*Sin[Pi/(2*r)*y]/ (Sqrt[(r+.000001)^2 - (r*Sin[Pi/(2*r)*y])^2])]-3*Pi/2, ArcTan[r*Sin[Pi/(2*r)*y]/ (Sqrt[(r+.000001)^2 - (r*Sin[Pi/(2*r)*y])^2])]-Pi/2}]}, {y,-r,r,2*r/(2^r)}]//N;
Show[Graphics[{ arcs,Circle[{0,0},r+1] }],Axes -> True, PlotRange->{{-10,10},{-10,10}}, AspectRatio -> Automatic, Ticks -> None]
},{r,1,6,1}]

(Refract.ma)

When light passes form one medium to another, which occurs in our everyday life, part of the light is bounced off (reflection) and the rest passes through the medium with a change in its direction (refraction). The amount of light bounced off is dependent on the angle of incidence and the refractive indicies of the two media. Snell' s Law allows us to determine the angle of refraction if the angle of incidence and the refractive index is known. Fresnel equations are used to determine the intensity of light reflected and refracted.

Exercises:
1. Create an animation using maple and mathematica.
2. Develop a spread sheet program to show the intensity of light reflected and refracted when
a) light passes from a denser medium to a rarer medium and
b) when light passes from a rarer medium to a denser medium.

Input Parameters
n1 & n2 are the two indicies of refraction
th eta1 is the angle of incidence
x0 is a starting point form which the ray is drawn (for best results, x0 should remain at -5)
Clear["Global`*"]
n1 = 1.2;
n2 = 2;
theta1 = 30;
x0 = -5;

Calculate theta2 (Snell's Law)
theta2 = ArcSin[n1*Sin[(theta1*Pi)/180]/n2]//N

Animate two lines: one for the incoming ray and one for the refracted ray. Note that the refracted ray is drawn on top of the incoming ray.
Do[ Show[Graphics[{
Line[{{If[step <= -x0,x0,0],-x0*Tan[Pi/2 - theta1*Pi/180] },{If[step <= -x0,x0+step,0], -(x0+step)*Tan[Pi/2 - theta1*Pi/180]}}],
Line[{{If[step >= -x0,x0,0],-x0*Tan[Pi/2 - theta1*Pi/180] },{If[step >= -x0,0,0],0},
{If[step >= -x0,(x0+step),0], -(x0+step)*Tan[Pi/2 - theta2]}}] }],
Axes -> True, PlotRange -> {{-8,8},{-8,8}}, AspectRatio -> Automatic, Ticks -> None],
{step,0,-2*x0,(-x0/10)}]

(Refrref2.ma)

Similar to the problem above (Refract), this progam also draws a reflected ray in addition to a refracted ray.
There are several possible uses for this problem; and with a few changes this problem can be as hard or as easy as the instructor desires.

One posible modifiaction is to have the angle of incidence change as the animation runs. This would allow students to get a feel for what happens as the angle approaches the critical angle. However, should the angle be allowed to be greater than the critical angle, the current code would have to be changed, and/or new code would have to be written.  
Another possible use is involves Stokes Relations. By assigning values to the different amplitudes, studets could explore Stokes Relations and how they are related to Snell's law.
This problem could also be used in discussions about polarization, transverse electric and transverse magnetic modes, and Brewster's angle.

Clear["Global`*"]

Input the two indicies of refraction
n1 = 1; n2 = 1.5;

Calculate the critical angle (in degrees)
CriticalAngle = If[n2 > n1, ArcSin[n1/n2]*(180/Pi), ArcSin[n2/n1]*(180/Pi)]//N

Now input the desired angle of incidence (measured from the positive y-axis). NOTE: If the incidence angle is greater than the critical angle, total internal reflection will occur, and the program will output several error messages in between the graphs.
theta1 = 30;

Now theta2, the angle of refraction, can be calculated (in radians)
theta2 = ArcSin[n1*Sin[(theta1*Pi)/180]/n2]//N

Input a starting value for the x cordinate for the ray (a number between -8 & 0)
x0 = -3;

Animate the three lines. The complicated logic statements (multiple-layer IF statements) allow for the three lines to all be independent of each other. (None of them are drawn on top of each other like in Refract)
Do[ Show[Graphics[{
Thickness[.009],RGBColor[0,0.5,0], Line[{ {If[step <= -x0,x0,x0], If[step <= -x0,-x0*Tan[Pi/2 - theta1*Pi/180], -x0*Tan[Pi/2 - theta1*Pi/180]]}, {If[step <= -x0,x0+step,0], If[step <= -x0,-(x0+step)*Tan[Pi/2 - theta1*Pi/180],0]} }],
RGBColor[.2,.1,1], Line[{ {If[step >= -x0,0,0],0}, {If[step >= -x0,(x0+step),0], If[step >= -x0,-(x0+step)*Tan[Pi/2 - theta2],0]} }],
RGBColor[1,1,0], Line[{ {If[step >= -x0,0,0],0}, {If[step >= -x0,(x0+step),0], If[step >= -x0,(x0+step)*Tan[Pi/2 - theta1*Pi/180],0]} }], GrayLevel[0],PointSize[.0084], Point[{0,0}]
}],
Axes -> True, PlotRange -> {{-8,8},{-8,8}}, AspectRatio -> Automatic, Ticks -> None ],
{step,0,-2*x0,(-x0/10)}]

(Sphmir.ma)

Mirrors can be one of the simplest optical tools that we use every day but they can also be the msot complicted tools. A common example is a plane mirror, where as the side mirrors in an automobile are more complex. These are mirrors used in the headlight of a car which are carefully designed. The shape of the mirror and its size dictate the specific application of the mirror. The maple/mathematica exercises allow the students to experience in a simple way one of the problems in the imaging property of spherical and parabolic mirrors. The abberation the students will learn about is spherical abberation.

Exercises:
1. How far away do two rays, one close to the optical axis (paraxial ray) and the other, at quite a distance from the optic axis (marginal ray), that are parallel to the optical axis, intersect the optical axis for a spherical mirror.
2. What happens if the mirror is parabolic in shape?

Clear["Global`*"]

Radius of curvature of the mirror (best results obtained when greater than 50)
R = 60;

Height of the incoming ray from the optical axis
y = 10;

Compute DeltaX (distance along the optical axis from where the ray hits the mirror to the center of the mirror)
DeltaX = R - Sqrt[R^2 - y^2]

By performing several algebra steps, x' can be calculated.
R' = R - DeltaX
Theta = ArcSin[y/R]
Phi = ArcCos[y/R]
Theta' = Theta - Phi
x' = y*Tan[Theta']

Turn off one of Mathematica's warning messages which has no effect on the output.
Off[Graphics::gprim]

Plot several incoming, parallel rays to show that they all don't cross at the focal point.
Show[Graphics[{
Table[{
Circle[{0,0},R,{-Pi/3,Pi/3}],
Hue[y/19],
Line[{ {-100,y}, {(R-(R-Sqrt[(R^2 - y^2)])),y}, {(R-(R-Sqrt[(R^2 - y^2)]))+(y*Tan[ArcSin[y/(R)] - ArcCos[y/(R)]]), 0}, {(R-(R-Sqrt[(R^2 - y^2)])) + 6*(y*Tan[ArcSin[y/(R)] - ArcCos[y/(R)]]), -5*y} }] }, {y,4,32,2}],
}],
Axes -> {True,False}, AspectRatio -> Automatic, Ticks -> {{{10^(-9),"R",.03},{R/2,"f",.03}}, None}, PlotRange -> {{-30,70},{-40,40}}]

(Parmir.ma)

Parmir shows that parallel rays can all meet at the focal point, so long as the mirror is parabolic. See Sphmir for more details.

Clear["Global`*"]

Radius of curvature of the mirror (best results obtained when greater than 50)
R = 60;

Height of the incoming ray from the optical axis
y = 20;

By performing several algebra steps, x' can be calculated.
DeltaX = y^2/(2*R)//N
Dis = Sqrt[(R - y^2/(2*R))^2 + y^2]//N
Theta = ArcSin[y/Sqrt[(R - y^2/(2*R))^2 + y^2]]//N
H = Sqrt[(R/2)^2 + Dis^2 - R*Dis*Cos[Theta]]//N
Phi = ArcCos[y/H]
x' = y*Tan[Phi]

Compare x' plus DeltaX to ROverTwo, the location of the focal point.
ROverTwo = y^2/(2*R) + y*Tan[ArcCos[y/(Sqrt[(R/2)^2 +
(Sqrt[(R - y^2/(2*R))^2 + y^2])^2 - R*(Sqrt[(R - y^2/(2*R))^2 + y^2])*Cos[Theta]])]]//N

Turn off three of Mathematica's warnings which do not affect the outcome of the graphs.
Off[Graphics::gprim] Off[Infinity::indet] Off[Graphics::gptn]

Show that several incoming, parallel rays will all meet at the focal point. (Note that the code here assumes that the lines will cross at the focal point; a good exercise for students would be to have them modify the code so that it proves that the parallel lines will all cross at the focal point.)
Plot[{Sqrt[2*R*x],-Sqrt[2*R*x]},{x,0,100},
Prolog -> Table[{
Hue[y/19],
Line[{
{100,y}, {y^2/(2*R),y}, {y^2/(2*R) + y*Tan[ArcCos[y/(Sqrt[(R/2)^2 +
(Sqrt[(R - y^2/(2*R))^2 + y^2])^2 - R*(Sqrt[(R - y^2/(2*R))^2 + y^2]) *Cos[ArcSin[y/Sqrt[(R - y^2/(2*R))^2 + y^2]]]])]]//N,0}
}] }, {y,-R,R,R/10}],
Axes -> {True,False}, AspectRatio -> Automatic,
Ticks -> {{{R,"R",.03},{R/2,"f",.03}}, None}, PlotRange -> {{-10,R+30},{-(R+10),(R+10)}}]

Refraction at any 2D lens surface.

Turn off some of Mathematica's warning messages and load the Implicit Plot package.
Clear["Global`*"]
Needs["Graphics`ImplicitPlot`"]
Off[General::spell1]
Off[General::spell]

Input parameters
1) The equation for the surface
2) The indicies of refraction for the lens and the surroundings
surface = 3*(x^2) + y^2 == 71*x;
n1 = 1.8;
n2 = 1.5;

Pick a point on the surface of the lens to be the point of intersection.
poi = {4,Sqrt[236]}

Test to see if the point is on the surface. A result of "True" means the point is on the lens surface. The point must be on the surface for the program to work.
surface/.{x -> poi[[1]], y-> poi[[2]]}//N

Find the normal slope at the point of intersection.
m = Dt[surface,x];
dydx = Solve[m, Dt[y,x]];
slope = dydx[[1,1,2]];
normal = -1/slope;
nslope = normal/.{x -> poi[[1]], y -> poi[[2]]}//N

Line equation for the normal.
normalline = (y - poi[[2]]) == nslope*(x - poi[[1]])//N

Now an equation for the ray must be created. The one existing condition is that the ray line must pass through the point of intersection. Thus, only a slope is needed to create the equation.
rslope = .2;
rayline = (y - poi[[2]]) == rslope*(x - poi[[1]])//N

The angle between the two lines (i1) must now be found. The logic statemsnts insure that the correct angle is found.
i1 = If[Abs[ArcTan[(nslope-rslope)/(1+rslope*nslope)]] < Pi/2, Abs[ArcTan[(nslope-rslope)/(1+rslope*nslope)]], Pi - Abs[ArcTan[(nslope-rslope)/(1+rslope*nslope)]]]//N

Calculate the angle of refraction (i2) by use of Snell's Law.
i2 = ArcSin[(n1*Sin[i1])/n2]

Use this angle to derive the slope (and line equation) of the ray beyond the lens surface.
newslope = (nslope - Tan[-i2])/(1 + Tan[-i2]*nslope)
refractline = (y - poi[[2]]) == newslope*(x - poi[[1]])//N

Use the line equations (rayline & refractline) to come up with starting and stoping points so that Mathematica's Line function can be used.
xstart = -100;
ystart = Solve[rayline/.x->xstart][[1,1,2]]
xstop = 100;
ystop = Solve[refractline/.x->xstop][[1,1,2]]

Generate the different graph types, then display all of them at once.
graph1 = Graphics[Line[{{xstart,ystart},{poi[[1]],poi[[2]]},{xstop,ystop}}]];
graph2 = ImplicitPlot[{surface,normalline},{x,poi[[1]]-10,poi[[1]]+10}, {y,poi[[2]]-10,poi[[2]]+10},
DisplayFunction -> Identity];
Show[{graph1,graph2}, PlotRange -> {{poi[[1]]-20,poi[[1]]+20}, {poi[[2]]-20,poi[[2]]+20}}, Axes -> True, AspectRatio -> Automatic]

(Lens7.ma)

Clear["Global`*"]

Opens the Mathematica package ImplicitPlot
Needs["Graphics`ImplicitPlot`"]

Turn off several Mathematica warning functions that do not change the output of the program.
Off[Solve::ifun]
Off[LinearSolve::sing]
Off[General::stop]
Off[LinearSolve::luc]
Off[General::spell1]

Initial Conditions
r1 & r2 are the two radii of the lens (should be positive); n is the index of refraction for the lens;
NumberOfRays is the number of rays to be drawn through the lens.
r1 = 60;
r2 = 80;
n = 1.8;
NumberOfRays = 15;

Calculations for drawing the lens

Parametric equations for lens surface 1
xr1 = -r1*Cos[t1]+Min[r1,r2]/1.5;
yr1 = r1*Sin[t1];

Parametric equations for lens surface 2
xr2 = r2*Cos[t2]-Min[r1,r2]/1.5;
yr2 = r2*Sin[t2];

Find the two points of intersection for the two surfaces.
int = FindRoot[{xr1==xr2, yr1 == yr2}, {t1,{.1,Pi}},{t2,{.1,Pi}}, MaxIterations -> 3*Max[r1,r2]]//N;

Calculate the lens height
LensHeight = yr1/.t1 -> int[[1,2]];

Exact Ray Trace

rh = Table[{
i1 = Solve[RayHeight == yr1, t1]//N;
i1 = i1[[1,1,2]];
i2 = ArcSin[Sin[i1/n]]//N;
t' = r1 - Sqrt[r1^2 - LensHeight^2];
deltat' = r1 - Sqrt[r1^2 - RayHeight^2]//N;
t1' = t' - deltat';
b = t1'*Tan[i1-i2];
m = -b/t1';
x1 = xr1/.t1->i1;
line = y == m*(x-x1)+RayHeight//Simplify;
circle = (x + Min[r1,r2]/1.5)^2 + y^2 == (r2)^2;
at = Solve[{line,circle},{x,y}];
xspot1 = xr1/.t1 -> i1,
yspot1 = RayHeight,
xspot2 = at[[1,2,2]],
yspot2 = at[[1,1,2]],
T = Sqrt[(xspot1-xspot2)^2 + (yspot1-yspot2)^2];
d = T*Sin[i1-i2];
i3 = ArcSin[(RayHeight-d)/r2] + (i1-i2);
i4 = ArcSin[n*Sin[i3]];
i5 = Abs[i4 - ArcSin[(RayHeight-d)/r2]];
xdis = yspot2/Tan[i5]
},{RayHeight,.1,LensHeight,LensHeight/NumberOfRays}];

Generate (but do not display) the plot for the first lens surface.
SurfaceOne = ParametricPlot[{xr1, yr1}, {t1,-int[[1,2]], int[[1,2]]}, DisplayFunction -> Identity];

Generate (but do not display) the plot for the second lens surface.
SurfaceTwo = ParametricPlot[{xr2, yr2}, {t2,-int[[2,2]], int[[2,2]]}, DisplayFunction -> Identity];

Generate (but do not display) the plot for the rays.
Rays = Table[Graphics[{Hue[r/NumberOfRays],Line[{{-100,rh[[r,2]]},{rh[[r,1]],rh[[r,2]]}, {rh[[r,3]],rh[[r,4]]},{rh[[r,5]]+rh[[r,3]],0},{rh[[r,3]]+(3*rh[[r,5]]),-2*rh[[r,4]]}}]},
DisplayFunction -> Identity], {r,1,NumberOfRays,1}];

Display all three plots together.
Show[ SurfaceOne,SurfaceTwo,Rays, AspectRatio -> Automatic,
PlotRange -> {{-1.4*Max[r1,r2],2*Max[r1,r2]}, {-1*Max[r1,r2],1*Max[r1,r2]}},
DisplayFunction -> $DisplayFunction, Ticks -> Automatic, Axes -> {True, False}]

(Basinfr.ma)

Interference, the superposition of waves, allows us to observe shimmering colors from oil spills on wet surfaces, complicated ripples in a pool of water, to a low frequency beat when tuning two musical instruments. Using a point source that emit waves radially, the total intensity from the two sources on a screen placed a fixed distance from the source can be determined. The following problems can be generated from this idea.
a) Obtain the interference pattern of two point sources by creating two independent sources of light. Calculate the irradiances due to each independent source, the irradiance due to the two sources combined, the maximum and minimum irradiances, and the fringe constant. Graph the interference pattern in terms of the phase difference phi.
b) Use the interference pattern of part a) in a Young's Double slit experiment. Graph the interference pattern in terms of y, the distance along the screen. Perform a check (m*lambda = d*sin(theta)) on your work.

Part A
Clear["Global`*"]
Off[General::spell1]
Off[General::spell]

Generate two point sources of the form E0*Cos[wt + e0]
E1 = E01*Cos[w1*t + e1];
E2 = E02*Cos[w2*t + e2];
E01 = 4;
E02 = 2;
w1 = 3*Pi;
w2 = 2*Pi;
e1 = Pi/3;
e2 = Pi/4;

Comput the irradiances due to each independent source. (Note: the constants e0 and c have been left off)
I1 = (1/2)*E01^2 I2 = (1/2)*E02^2

Compute the interference term (I12) in terms of the phase difference phi.
I12 = 2*Sqrt[I1*I2]*Cos[phi]

The total irradiance is thus:
Itotal = I1+I2+I12

Compute the maximum and minimum irradiances & the fringe constant
Imax = Itotal/.{phi->0}
Imin = Itotal/.{phi->Pi}
fc = (Imax - Imin)/(Imax + Imin)

Graph the interference pattern
Plot[Itotal,{phi,-5*Pi,5*Pi}, PlotRange -> {Imin-1, Imax+1}, AxesOrigin -> {0,Imin-1}, AxesLabel -> {"Phi","Irradiance"}]

Part B
Assume now that this is the interference pattern for a Young's double-slit experment with two light sources, seperated by .1 mm, at 550 nm and a screen distance of 1 m. Graph the intereference pattern in terms of y (the distance along the screen).
d = .1*10^(-3);
L = 1;
lambda = 550*10^(-9);
phi = 2*Pi*(y*d)/(L*lambda)

Plot[Itotal,{y,-.01,.01}, PlotRange -> {Imin-1, Imax+1}, AxesOrigin -> {0,Imin-1}, AxesLabel -> {"y","Irradiance"}]

Prove that this graph is correct. (i.e. show that the known conditions at the maximum and/or minimum are true)
FindMinimum[Itotal,{y,{.002,.006}}]
Equal[d*.00275/L == (1-(1/2))*lambda]

(Cpxinfr.ma)

A more advance look at interference; using the ideas formulated in basinfr:
a) Compare the interference patterns for two and three point sources (all sources have equal amplitude)
b) Change the amplitude in one of the three sources; what happens? (what happens when the amplitudes of the sources are the binomial coefficients?)
c) Change the number of point sources from three to ten and note the change. Keep the distance between the sources constant. (Multiple Interference)
d) Make the spacing between the point sources in the case of three point sources unequal and observe the changes. If the distance between the sources is different by a factor of fifteen, what do you observe? (Principle of carrier waves)

Solution

Load one of Mathematica's packages
Clear["Global`*"]
Needs["Algebra`Trigonometry`"]

Definitions:
d is the spacing between the sources (micrometers)
L is the wavelength of the light used (micrometers)
Z is the distance of observation
t is the distance along the screen
d = 6.33;
L = .633;
Z = 10;
t = y;

The point source at the center is a[1]. The even numbered sources lie above a[1] and the odd numbered sources lie below a[1]. The distances of the even numbered sources are positive; the odd sources are at a negative distance.
The two Table functions produce 5 point sources each. The top function produces the odd numbered sources.
Table[a[2*n-1] = Exp[-I*(n-1)*(Pi*d*t)/(L*Z)], {n,1,5,1}];
Table[a[2*n] = Exp[I*(n)*(Pi*d*t)/(L*Z)], {n,1,5,1}];

Part A
The interfernce pattern can be obtained by squaring the sum of the point sources.
Two point sources
TwoSources = Abs[Sum[a[n1],{n1,1,2}]]
TwoSourcesSquared = Simplify[TwoSources^2/2^4]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll Plot[TwoSourcesSquared,{y,0,8}]

Three point sources (Same procedure as for two sources)
ThreeSources = Abs[Sum[a[n1],{n1,1,3}]];
ThreeSourcesSquared = Simplify[ThreeSources^2/2^4]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll;
Plot[ThreeSourcesSquared, {y,0,8}]

Part B
Change the amplitude in one of the point sources (try different sources)
DifferentAmplitudes = Abs[a[1]+4*a[2]+a[3]]
DifferentAmplitudesSquared = Simplify[DifferentAmplitudes^2/2^4]/.{Rule -> Equal}// ComplexToTrig//ExpandAll Plot[DifferentAmplitudesSquared,{y,0,8}]

Part C
Ten point sources (Again, the same procedure as two or three sources)
TenSources = Abs[Sum[a[n1],{n1,1,10}]];
TenSourcesSquared = Simplify[TenSources^2/2^7]/. {Rule -> Equal}//ComplexToTrig//ExpandAll; Plot[TenSourcesSquared,{y,0,8}, PlotPoints -> 100]

Part D
Make the point spacing unequal. This is accomplished by multiplying the 'd' term by a constant
a[1] = Exp[-I*(1-1)*(Pi*d*t)/(L*Z)];
a[2] = Exp[I*(2-1)*(Pi*15*d*t)/(L*Z)];
a[3] = Exp[-I*(2-1)*(Pi*d*t)/(L*Z)];
SpacingUnequal = Abs[Sum[a[n1],{n1,1,3}]]; SpacingUnequalSquared = Simplify[SpacingUnequal^2/2^7]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll; Plot[SpacingUnequalSquared,{y,0,2}, PlotPoints -> 100]

(Difintr.ma)

Displaying diffraction as interference of equispaced point sources.

By increasing the number of point sources, the interfernce pattern produced can be used to show diffraction.

Solution

Load one of Mathematica's packages
Clear["Global`*"]
Needs["Algebra`Trigonometry`"]
Definitions:
d is the spacing between the sources (micrometers)
L is the wavelength of the light used (micrometers) Z is the distance of observation
t is the distance along the screen
d = 6.33;
L = .633;
Z = 10;
t = y;

The point source at the center is a[1]. The even numbered sources lie above a[1] and the odd numbered sources lie below a[1]. The distances of the even numbered sources are positive; the odd sources are at a negative distance.
The two Table functions produce 25 point sources each. The top function produces the odd numbered sources.
Table[a[2*n-1] = 2*Exp[-I*(n-1)*(Pi*d*t)/(L*Z)], {n,1,25,1}];
Table[a[2*n] = 2*Exp[I*(n)*(Pi*d*t)/(L*Z)], {n,1,25,1}];

The interfernce pattern can be obtained by squaring the sum of the point sources.

Two point sources
TwoSources = Abs[Sum[a[n1],{n1,1,2}]]
TwoSourcesSquared = Simplify[TwoSources^2/(2*2)^2]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll Plot[TwoSourcesSquared,{y,0,8}]

Three point sources (Same procedure as for two sources)
ThreeSources = Abs[Sum[a[n1],{n1,1,3}]];
ThreeSourcesSquared = Simplify[ThreeSources^2/(2*3)^2]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll; Plot[ThreeSourcesSquared, {y,0,8}]

Four point sources
FourSources = Abs[Sum[a[n1],{n1,1,4}]];
FourSourcesSquared = Simplify[FourSources^2/(2*4)^2]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll; Plot[FourSourcesSquared, {y,0,8}]

Seven point sources
SevenSources = Abs[Sum[a[n1],{n1,1,7}]];
SevenSourcesSquared = Simplify[SevenSources^2/(2*7)^2]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll; Plot[SevenSourcesSquared, {y,0,8}, PlotRange -> {0,1}]

Fifteen point sources
FifteenSources = Abs[Sum[a[n1],{n1,1,15}]];
FifteenSourcesSquared = Simplify[FifteenSources^2/(2*15)^2]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll; Plot[FifteenSourcesSquared, {y,0,4}, PlotPoints -> 50, PlotRange -> {0,1}]

Twenty-five point sources
TwentyFiveSources = Abs[Sum[a[n1],{n1,1,25}]]
TwentyFiveSourcesSquared = Simplify[TwentyFiveSources^2/(2*25)^2]/.{Rule -> Equal}//ComplexToTrig//ExpandAll Plot[TwentyFiveSourcesSquared, {y,0,2}, PlotPoints -> 50, PlotRange -> {0,1}]

Fifty point sources
FiftySources = Abs[Sum[a[n1],{n1,1,50}]];
FiftySourcesSquared = Simplify[FiftySources^2/(2*50)^2]/.{ Rule -> Equal}//ComplexToTrig//ExpandAll; Plot[FiftySourcesSquared, {y,0,2}, PlotPoints -> 50, PlotRange -> {0,1}]

Now compare the 25 and 50 point-source interfernce patterns to the normal diffraction pattern for a rectangular aperture.
Note: the 25 and 50 point sources must be recalculated so that they are in phase with each other.

25 Sources
Table[s25[2*n-1] = 2*Exp[-I*(n-1)/25*(Pi*d*t)/(L*Z)], {n,1,25,1}];
Table[s25[2*n] = 2*Exp[I*(n)/25*(Pi*d*t)/(L*Z)], {n,1,25,1}];
source25 = Abs[Sum[s25[n1],{n1,1,25}]];
source25sq = Simplify[ source25^2/(2*25)^2]/.{Rule -> Equal}//ComplexToTrig//ExpandAll;

50 Sources
Table[s50[2*n-1] = 2*Exp[-I*(n-1)/50*(Pi*d*t)/(L*Z)], {n,1,50,1}];
Table[s50[2*n] = 2*Exp[I*(n)/50*(Pi*d*t)/(L*Z)], {n,1,50,1}];
source50 = Abs[Sum[s50[n1],{n1,1,50}]];
source50sq = Simplify[ source50^2/(2*50)^2]/.{Rule -> Equal}//ComplexToTrig//ExpandAll;

Normal diffraction pattern
norm = (Sin[(.5*Pi*d*t)/(L*Z)]/((.5*Pi*d*t)/(L*Z)))^2

Plot all 3 functions together.
Plot[{norm,source25sq, source50sq}, {y,10^(-8),8}, PlotPoints -> 50, PlotRange -> {0,1}]

(Partcoh1.ma)

This program uses double slits with finite seperation and slit widths and an extended source.
We generate partial coherence through a double slit by adding incoherently interference patterns from parts of an extended source.
a = width of either slit
d = center-to-center spacing of slits lam = wavelength of light
k = 2*Pi/lam
w = half-width of extended source
L1 is the distance of the screen
Ls is the distance of the source from the double slit
Ys is the distance of the source above the z-axis
y is the coordinate at the screen

Clear["Global`*"]

lam = 1;
k = 2*Pi/lam;
a = 1;
d = 4;
ga = k*a/2*(y/L1)
gd = k*d/2*(y/L1 + ys/Ls)

L1 = 10;
Ls = 5;
ys = 0;
f = ((Sin[ga]/ga)*Cos[gd])^2

Plot[f,{y,-20,20}, PlotRange -> {0,1}, PlotPoints -> 40, PlotLabel -> "Double Slit, Full Coherence", AxesLabel -> {"y","Irradiance"}]

ga = k*a/2*(y/L1);
gd = k*d/2*(y/L1 + ys1/Ls);
f1 = ((Sin[ga]/ga)*Cos[gd])^2

Off[NIntegrate::ploss] ListPlot[Table[{y,(1/.8)*NIntegrate[f1,{ys1,-.4,.4}]},{y,-20,20,.08}], PlotJoined -> True, PlotRange -> {0,1}, PlotLabel -> "Partial Coherence, w = +/- 0.8"]

(Arago2.ma)

Plotting the focii of a zone plate.

Clear["Global`*"]
L = wavelength of light
r0 = radius of the first zone
zon = number of zones
L = 0.633;
r0 = 1.8;
zon = 15;

Create a table that generates the radii of the zones (Don't worry when the answers come back "Null", that's what their supposed to do.)
Table[{r1[n_] := r1[n] = Sqrt[r0*L*n+n^2*L^2/4]}, {n,1,2*zon}]
Table[{a[m_] := a[m] = (40/Sqrt[z^2+(y-r1[2*m-1])^2]) * Exp[I*2*Pi/L* Sqrt[z^2+(y-r1[2*m-1])^2]]},{m,1,zon}]

f3 = Abs[Sum[a[n1],{n1,1,zon}]] Sqrt[(-1.11336 + y) + z ]

Needs["Algebra`Trigonometry`"]
f4 = Simplify[f3^2]/.{Rule -> Equal}//ComplexToTrig//ExpandAll

Plot3D[f4, {z,9,30}, {y,-2,2}, PlotPoints -> 30, PlotRange -> {0,200}, AxesLabel -> {"Z","Y","Irradiance"}]

(Zoneamp.ma)

STILL UNDER CONSTRUCTION

Clear["Global`*"]

L = 0.633; r = 5; fzon = 1; zon = 8;

Table[{r1[n_] := r1[n] = Sqrt[r*L*n+n^2*L^2/4]}, {n,1,2*zon}]
Table[{a[m_] := a[m] = (40/Sqrt[z^2+(y-r1[2*m-1])^2]) * Exp[I*2*Pi/L*Sqrt[z^2+(y-r1[2*m-1])^2]]}, {m,fzon,zon}]

f3 = Abs[Sum[a[n1],{n1,fzon,zon}]]

Needs["Algebra`Trigonometry`"]
f4 = Simplify[f3^2]/.{Rule -> Equal}//ComplexToTrig//ExpandAll

y = 0;
Plot[f4, {z,0,30}, PlotRange -> {0, 7500}]

(Zoneplat2.ma)

Create a program that draws both positive and negative zone plates of different sizes.

Load one of Mathematica's packages
Clear["Global`*"]
Off[Graphics::gprim]

Initial radius & wavelength
r0 = 5*10^(-2);
lambda = 550*10^(-9);

Zone Plate Type & Size
Note: Type = 1 (Positive Plate); Type = -1 (Negative Plate)
The larger the size, the more rings the plate will have.
Size must be an integer greater than 1.
type = -1;
size = 35;

Plot the zone plate
Show[Graphics[{
Table[{ GrayLevel[If[type==1,If[EvenQ[N],0,1],If[type == -1,If[EvenQ[N],1,0],.5]]], Disk[{0,0},Sqrt[N*r0*lambda]] },{N,size,1,-1}]
}],
AspectRatio -> Automatic]