Birefringence in dispersive media

by Evi Postelnicu

I. Interference condition

Light is comprised of electromagnetic waves with electric and magnetic fields oscillating orthogonal to each other. Light propogates in the direction perpendicular to the electric and magnetic fields. When electromagnetic waves encounter obstacles or slits, their path is diverted, or the waves are diffracted. Diffraction is most pronounced when obstacles, or slits, are close in size to the wavelength. When light encounters a diffraction grating of N slits of width d, the condition that light interferes constructively on a screen a large distance away is that the difference in path length that the light travels is equal to an integer wavelength. The path difference between two rays diffracted from different slits at angle θ, d*sin(θ), is solved for using trigonometry, the distance between the slits, and the diffracted angle. The optical path length difference is just the physical path difference multiplied by ((2*π)/λ) , where λ is the wavelength of incident light. If light interferes constructively, a bright fringe is formed on the screen. Otherwise, if the path length does not abide by the condition for constructive interference, light interferes destructively and a dark fringe will appear on the screen. The intensity profile of the light on the screen can be calculated for the diffracted angle, theta, of light of wavelength lambda. The equation for the intensity of interference is derived from superposition of the electric fields of the interfering waves, and is proportional to the electric field squared. The intensity of N-slit interference is thus proportional to (((((sin(((N δ)/2)))/(sin((δ/2))))))^2), where δ = ((2×π×d×sin(θ))/λ). N-slit interference shows the greatest intensity for locations on the screen directly across from a screen because light travels the least to arrive at that location. However there are N-1 peaks surrounding the main peak which correspond constructive interference as well where the light has to travel further and so the intensity is lower.

delta[d_,theta_,lambda_]:=2*Pi*d*Sin[theta]/lambda

interferenceI[N_,d_,theta_,lambda_]:=(Sin[N*delta[d,theta,lambda]/2]/(delta[d,theta,lambda]/2))^2

Ex. Interference packet

Manipulate[Plot[interferenceI[N,180,angle,300]/N^2,{angle,-2*Pi,2*Pi},PlotRange-> {0,1},ImageSize->Large,Frame->True,FrameLabel->{Style["Angle[radians]",14],Style["Intensity[power/m^2]",14]}],{N,2,5}]

GraphicsRow[{Graphics[{Line[{{0,0},{0,2}}],Line[{{0,2.5},{0,4.5}}],Line[{{0,5},{0,7}}],Line[{{0,7.5},{0,9.5}}],Line[{{0,10},{0,12}}],Line[{{0,12.5},{0,14.5}}],Arrowheads[Medium],Arrow[{{0,2.25},{2,2.25}}],Arrow[{{0,4.75},{2,2.25}}],Arrow[{{0,7.25},{2,2.25}}],Arrow[{{0,9.75},{2,2.25}}],Arrow[{{0,12.25},{2,2.25}}]}],ImageRotate[Plot[interferenceI[5,150,angle,300]/25,{angle,-2*Pi,2*Pi},PlotRange-> {0,1}],Pi/2]}]

II. Single slit diffraction

The effects on the intensity profile of light passing through a single slit are analyzed in single slit diffraction. The behavior of light passing through one slit of width D is the same as that of light passing through two very narrow slits of width D/2. However, in this case the optical path length between the rays is now (D/2)*sin(θ) because the slits are half as wide and so the waves are diffracted less. The intensity profile is obtained in the same way as it was for interference, by summing up the individual components of the electric field and taking the square in order to have something proportional to intensity. For single-slit diffraction, the intensity is proportional to (((((sin(β))/β)))^2) where β = ((π×D×sin(θ))/λ). The optical path length difference, β, is 1/2 that of n-slit interference due to the two-slit construction out of one slit, as discussed previously and illustrated in the figure.

betaSingleSlit[widthOfSlit_,theta_,lambda_]:=Pi*widthOfSlit*Sin[theta]/lambda

diffractionI[widthOfSlit_,theta_,lambda_]:=(Sin[betaSingleSlit[widthOfSlit,theta,lambda]]/betaSingleSlit[widthOfSlit,theta,lambda])^2

GraphicsRow[{Graphics[{Line[{{0,0},{0,1}}],Line[{{0,2},{0,3}}],{Dashed,Line[{{-0.5,1.5},{2,1.5}}]},Arrowheads[Medium],Arrow[{{-0.5,1.55},{-0.5,1.9}}],Arrow[{{-0.5,1.45},{-0.5,1}}],Text[Style["D/2",14],{-0.3,1.7}],Text[Style["D/2",14],{-0.3,1.3}],Text[Style["θ",14],{0.1,1.7}],Text[Style["β",14],{0.3,1.1}],Dashed,Arrow[{{0,2},{2,2.3}}],Arrow[{{0,1},{2,2.}}],Line[{{0,2},{0.5,1.36}}],Line[{{0,2},{0,1}}]}],Rotate[Plot[diffractionI[600,angle,300],{angle,-(π/2),(π/2)},ImageSize->Large,Frame->True,FrameLabel->{Style["Angle[radians]",14],Style["Intensity[power/m^2]",14]}],(π/2)]}]

III. N-slit diffraction

The combination of both N-slit interference and single-slit diffraction provides the intensity profile for N-slit diffraction, or the intensity profile from a diffraction grating of N slits of width, D, and distance between slits, d. The intensity profile allows for constructive interference due to single-slit diffraction and N-slit interference. The N-slit diffraction profile is actually like a lot of little N-slit interference packets underneath the larger single-slit diffraction packet, as shown in the plot of nSlitDiffraction. The intensity of N-slit diffraction is thus simply I = (Subscript[I, 0])× (((((sin(β))/β)))^2)×(((((sin(((N δ)/2)))/(sin((δ/2))))))^2), where (Subscript[I, 0]) is the intensity of the central maximum.

nSlitDiffractionI[widthOfSlit_,d_,theta_,lambda_,N_]:=interferenceI[N,d,theta,lambda]*diffractionI[widthOfSlit,theta,lambda]

Plot[{interferenceI[3,1000,angle,300],nSlitDiffractionI[600,1000,angle,300,3],diffractionI[600,angle,300]},{angle,-Pi/2,Pi/2},PlotRange-> {0,1},PlotLegends-> {"N-slit interference","N-slit diffraction","single-slit diffraction"},ImageSize->Large,Frame->True,FrameLabel->{Style["Angle[radians]",14],Style["Intensity[power/m^2]",14]}]

IV. Diffraction grating equation

N-slit diffraction provides the intensity profile for various allowable diffracted angles from perpendicular incident light. The Bragg condition for constructive interference then provides the angles at which constructive interference actually occurs. The grating equation profits from this and relates the path length difference between the incident wave and the reflected wave to an integer number times the wavelength such that: m×λ = d (sin(α) + sin(β)). (see figure 2) The integer number is referred to as the order of diffraction. The diffraction grating equation tells us the direction in which light will travel after passing through a diffraction grating and being diffracted by the angle β given the incident angle, α, the distance between slits, d, and the wavelength of light. For reflective diffraction gratings, the convention is that the angle is positive if you precess counterclockwise to obtain the angle from the normal to the wavefront. If you rotate clockwise from the normal to the wavefront, the angle is positive. Angles are always measured with respect to the normal of the material (or in this case, the diffraction grating). The grating equation was used to solve for the diffracted angle, β, of different incident angles and wavelengths.

gratingEquation=m*lambda==d (Sin[alpha]+Sin[beta]);

V. Dispersion due to grating

Diffraction gratings act as prisms to disperse white light because of the different angles at which wavelengths are refracted. In order to determine where the light is dispersed, at what intensity, and the dependency of location on wavelength, I first solved the grating equation β in terms of the incident angle, α, the distance between slits, d, and the wavelength, λ. I then calculated the different positions different orders of the same wavelength would exist at by taking the tangent of the diffracted angle and multiplying by the distance to the screen we are dispersing light onto, which I let be one meter.

betaSol=beta/.Solve[gratingEquation,beta][[2]]/.C[1]->0

betaSols[order_,distanceBetweenSlits_,incidentAngle_,wavelength_]:=betaSol/.{m->order,d-> distanceBetweenSlits,alpha->incidentAngle,lambda->wavelength}

yDistance[order_,distanceBetweenSlits_,incidentAngle_,wavelength_]:=Tan[betaSols[order,distanceBetweenSlits,incidentAngle,wavelength]]

spectralIntensity[800,600,5,5000,1,-Pi/3]

Re[yDistance[2,5000,-Pi/3,800]]

Example: Dispersion of white light, where each part of light is one specific wavelength. Distance between slits, d = 5000 nm, width of slits, a = 800 nm, distance to screen = 1 m, number of slits = 5.

In order to show the dispersion accurately with the correct intensity for each spectral order, the spectral intensity needs to be calculated from the summation of the plane waves intefering from each pair of slits. A single slit in a far field regime creates the distribution of the electric field given by E = (Subscript[E, 0]) ×((sin(β/2))/(β/2))×cos(kr-ωt), β = 2 ((πa)/λ)sin(θ). The total electric field from N slits is the sum of the indiviudal contributions. As the electric field travels through the medium, it picks up a phase due to difference light path length. We can take the center slit to have 0 phase shift. Each subsequent slip will have a phase shift (Subscript[Δϕ, i]), such that (Subscript[E, i])=(Subscript[E, 0]) ×((sin(β/2))/(β/2))×cos(kr-ωt+(Subscript[Δϕ, i])). Consider an even number of slits, with the first pairs of slits are at distances ± d/2 from the center, the second pair at distances ± 3d/2, etc. The phase difference between paths from the centers of the slits and the center line are, (Subscript[Δϕ, 1])=±ϕ/2, (Subscript[Δϕ, 2])=±3ϕ/2,...,(Subscript[Δϕ, p])=±(p-(1/2))ϕ, where p is the pair number. In order to sum over each electric field from the slits, we sum over the contribution of two slits in each pair, and then sum over all the pairs. 
E= (Subscript[E, 0])×((sin(β/2))/(β/2))×(∑_(p=1)^(N/2))[cos(kr-wt+(Subscript[Δϕ, p]))-cos(kr-wt-(Subscript[Δϕ, p]))]
E= 2(Subscript[E, 0])×((sin(β/2))/(β/2))×(∑_(p=1)^(N/2))[cos((Subscript[Δϕ, p]))]*cos(kr-wt)
The spectral intensity is then I= 4(Subscript[I, 0])×(((((sin(β/2))/(β/2))))^2)×(([(∑_(p=1)^(N/2))[(2p-1)π×(d/λ)sin(θ)]])^2).

spectralIntensity[a_,wavelength_,Number_,distanceBetweenSlits_,order_,incidentAngle_]:=N[4*(Sin[Pi*Sin[betaSols[order,distanceBetweenSlits,incidentAngle,wavelength]]*(a/wavelength)]/(Pi*Sin[betaSols[order,distanceBetweenSlits,incidentAngle,wavelength]]*(a/wavelength)))^2*(Sum[Cos[(2*p-1)*Pi*distanceBetweenSlits*Sin[betaSols[order,distanceBetweenSlits,incidentAngle,wavelength]]/wavelength],{p,1,number/2}])^2/.{number->Number}]

color[400]=Purple;color[475]=Blue;color[510]=Green;color[570]=Yellow;color[590]=Orange;color[650]=Red;

Show[Graphics[{Line[{{0,.2},{0,.7}}],Line[{{0,-.3},{0,-0.8}}],Line[{{0,-1.3},{0,-1.8}}],Line[{{0,1.2},{0,1.7}}],Line[{{0,2.2},{0,2.7}}],Line[{{0,-2.3},{0,-2.8}}],Thin,Table[{Dashed,Opacity[1,Black],Line[{{3,10*yDistance[0,5000,0,wavelength]-0.25},{3,10*yDistance[0,5000,0,wavelength]+0.25}}],Thick,Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,0]/25,color[wavelength]],Line[{{3,10*yDistance[order,5000,0,wavelength]-0.25},{3,10*yDistance[order,5000,0,wavelength]+0.25}}]},{order,1,2,1}],1],Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,0]/25,color[wavelength]],Line[{{3,10*yDistance[order,5000,0,wavelength]-0.25},{3,10*yDistance[order,5000,0,wavelength]+0.25}}]},{order,-2,-1,1}],1]},{wavelength,{400,475,510,570,590,650}}]}],ImageSize->Large]

Show[Graphics[{Line[{{0,.25},{0,.75}}],Line[{{0,-.25},{0,-0.75}}],Line[{{0,-1.25},{0,-1.75}}],Line[{{0,1.25},{0,1.75}}],Line[{{0,2.25},{0,2.75}}],Line[{{0,-2.25},{0,-2.75}}],Thin,Table[{{Dashed,Circle[{5,5*yDistance[0,5000,0,wavelength]},.1]},Thin,Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,0]/25,color[wavelength]],Arrowheads[Small],Arrow[{{0,0},{5,0}}],Arrow[{{0,1},{5,0}}],Arrow[{{0,2},{5,0}}],Arrow[{{0,-1},{5,0}}],Arrow[{{0,-2},{5,0}}],Arrow[{{0,0},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,1},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,2},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,-1},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,-2},{5,5*yDistance[order,5000,0,wavelength]}}],Disk[{5,5*yDistance[order,5000,0,wavelength]},.1]},{order,1,3,1}],1],Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,0]/25,color[wavelength]],Arrowheads[Small],Arrow[{{0,0},{5,0}}],Arrow[{{0,1},{5,0}}],Arrow[{{0,2},{5,0}}],Arrow[{{0,-1},{5,0}}],Arrow[{{0,-2},{5,0}}],Arrow[{{0,0},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,1},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,2},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,-1},{5,5*yDistance[order,5000,0,wavelength]}}],Arrow[{{0,-2},{5,5*yDistance[order,5000,0,wavelength]}}],Disk[{5,5*yDistance[order,5000,0,wavelength]},.1]},{order,-3,-1,1}],1]},{wavelength,{400,475,510,570,590,650}}]}],ImageSize->Large]

options1=Table[Show[Graphics[{Table[{{Dashed,Circle[{5,5*yDistance[0,5000,incidentAngle,wavelength]},.1]},Thin,Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,incidentAngle]/25,color[wavelength]],Disk[{5,5*yDistance[order,5000,incidentAngle,wavelength]},.1]},{order,1,3,1}],1],Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,incidentAngle]/25,color[wavelength]],Disk[{5,5*yDistance[order,5000,incidentAngle,wavelength]},.1]},{order,-3,-1,1}],1]},{wavelength,{400,475,510,570,590,650}}]}],ImageSize->Large],{incidentAngle,-Pi/3,Pi/3,Pi/96}]

Animate[options1[[i]],{i,20,50,1}]

allOptions2=Table[Show[Graphics[{Line[{{0,.25},{0,.75}}],Line[{{0,-.25},{0,-0.75}}],Line[{{0,-1.25},{0,-1.75}}],Line[{{0,1.25},{0,1.75}}],Line[{{0,2.25},{0,2.75}}],Line[{{0,-2.25},{0,-2.75}}],Thin,Table[{{Dashed,Circle[{5,5*Re[yDistance[0,5000,incidentAngle,wavelength]]},.1]},Thin,Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,incidentAngle]/25,color[wavelength]],Arrowheads[Small],Arrow[{{0,0},{5,0}}],Arrow[{{0,1},{5,0}}],Arrow[{{0,2},{5,0}}],Arrow[{{0,-1},{5,0}}],Arrow[{{0,-2},{5,0}}],Arrow[{{0,0},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,1},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,2},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,-1},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,-2},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Disk[{5,Re[5*yDistance[order,5000,incidentAngle,wavelength]]},.1]},{order,1,3,1}],1],Flatten[Table[{Opacity[spectralIntensity[800,wavelength,5,5000,order,incidentAngle]/25,color[wavelength]],Arrowheads[Small],Arrow[{{0,0},{5,0}}],Arrow[{{0,1},{5,0}}],Arrow[{{0,2},{5,0}}],Arrow[{{0,-1},{5,0}}],Arrow[{{0,-2},{5,0}}],Arrow[{{0,0},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,1},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,2},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,-1},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Arrow[{{0,-2},{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]}}],Disk[{5,5*Re[yDistance[order,5000,incidentAngle,wavelength]]},.1]},{order,-3,-1,1}],1]},{wavelength,{400,475,510,570,590,650}}]}],ImageSize->Large],{incidentAngle,-Pi/2,Pi/2,Pi/196}]

Manipulate[allOptions2[[i]],{i,48,196-48,1}]

Animate[allOptions2[[i]],{i,48,196-48,1}]

IV. Birefringence

Anisotropic crystals also diffract light and disperse it like diffraction gratings, but due to a property called birefringence, light is actually reflected and polarized in two different directions and two different wavelengths. Birefringence refers to the optical property of a material where its index of refraction depends on the propagation and polarization of light. In uniaxial crystals, the system we will be analyzing, the indices of refraction along two of the principal axes are equal and fixed as (Subscript[n, ordinary]). The index of refraction along the third principal axis is called (Subscript[n, extraordinary]) and varies with light polarization and angle. The z axis is known as the optic axis. The electric flux density, D, is given by  ϵ*E. In birefringent materials, epsilon - the electrical permittivity, is not a constant, but rather a tensor - a 3x3 matrix.  Usually this tensor is symmetric, so there are at most six coefficients defining the electric permittivity. If your coordinate axes is chosen along the principal axes of symmetry in the crystal, the diagnoal components of the matrix vanish and you only have three coefficients of electric permittivity. To simplify calculations, our coordinate system will always be along the principal axes. The permittivities (Subscript[ϵ, 1]), (Subscript[ϵ, 2]), (Subscript[ϵ, 3]) correspond to the refractive indices (Subscript[n, 1])=(Sqrt[((Subscript[ϵ, 1])/(Subscript[ϵ, 0]))]),  (Subscript[n, 1])=(Sqrt[((Subscript[ϵ, 2])/(Subscript[ϵ, 0]))]) ,   (Subscript[n, 1])=(Sqrt[((Subscript[ϵ, 3])/(Subscript[ϵ, 0]))]), where (Subscript[ϵ, 0]) is the permittivity of free space.

electricFluxDensity[i_]:=Sum[epsilon[[i]][[j]]*electricField[[j]],{j,1,i}]

epsilon={{epsilon11,epsilon12,epsilon13},{epsilon12,epsilon22,epsilon23},{epsilon13,epsilon23,epsilon33}};

electricField={E1,E2,E3}

epsilon0=8.854*10^-12

Electric impermeability tensor

The electric impermeability tensor, η = (Subscript[ϵ, 0])(ϵ^(-1)), is the permittivity of free space times the inverse of the electric permittivity tensor. The quadratic representation of the electric impermeability tensor, (∑_ij)(Subscript[η, ij])(Subscript[x, i])(Subscript[x, j])=1, where i,j = 1, 2, 3, is called the index ellipsoid and is used to determine the index of refraction of an anisotropic solid along an arbitrary propagation direction. The propagation direction’s normal modes’ indicies of refraction, (Subscript[n, a]) and  (Subscript[n, b]), will contain components from the principal indicides of refraction and will be determined by the intersection of the ellipsoid with a plane normal to the propagation direction that passes through the origin of the ellipsoid. The electric flux density, D, which will give the propagation direction of the normal modes is also determined through the index ellipsoid as shown in the figure.

eta=epsilon0*Inverse[epsilon/.{epsilon12->0,epsilon13->0,epsilon23->0}]

position={x1,x2,x3};

N1=Sqrt[epsilon11/epsilon0]

N2=Sqrt[epsilon22/epsilon0]

N3=Sqrt[epsilon33/epsilon0]

indexEllipsoid=(Sum[eta[[i]][[j]]*position[[i]]*position[[j]],{i,1,3,1},{j,1,3,1}]/.{epsilon0/epsilon11->(1/n1^2),epsilon0/epsilon22->(1/n2^2),epsilon0/epsilon33->(1/n3^2)})==1

x3Sols=x3/.Solve[((x1^2)/(n1^2))+((x2^2)/(n2^2))+((x3^2)/(n3^2))==1,x3]

normalVector={a,b,c}

eqOfPlane=normalVector[[1]](x1-0)+normalVector[[2]] (x2-0)+normalVector[[3]] (x3-0)==0

plane=x3/.Solve[eqOfPlane,x3]

Plot3D[{x3Sols[[1]],x3Sols[[2]],plane}/.{n1->1.2,n2->1.3,n3->2.7,a->1,b->2,c->1},{x1,-2,2},{x2,-2,2},PlotStyle->{Directive[Opacity[0.6],Blue]},ImageSize->Large]

normalVector={a,b,c}

eqOfPlane=normalVector[[1]](x1-0)+normalVector[[2]] (x2-0)+normalVector[[3]] (x3-0)==0

plane=x3/.Solve[eqOfPlane,x3]

crossSectionTop=Solve[{plane==x3Sols[[1]]},{x1,x2}]

Show[Plot[{plane/.crossSectionTop[[1]]/.{n1->1.2,n2->1.3,n3->2.7,a->1,b->2,c->1},plane/.crossSectionTop[[2]]/.{n1->1.2,n2->1.3,n3->2.7,a->1,b->2,c->1}},{x1,-2,2}],Graphics[{Text[ Style["(Subscript[n, a]),(OverVector[(Subscript[D, a])])",14],{-1.1,1.9}],{Dashed,Line[{{0,0},{-.92,1.6}}]},Text[ Style["(Subscript[n, b]),(OverVector[(Subscript[D, b])])",14],{1,1.2}],{Dashed,Line[{{0,0},{.8,1}}]}}]]

Uniaxial crystal

For a uniaxial crystal, (Subscript[n, 1])=(Subscript[n, 2])=(Subscript[n, 0]), the ordinary axis, and (Subscript[n, 3])=(Subscript[n, e]), the extraordinary axis. For a wave whose direction of travel, (OverHat[u]), forms an angle θ with the optic axis (either direction of (Subscript[n, 0]), the index ellipse gives normal modes that have refreactive indicies (Subscript[n, b]) = (Subscript[n, 0]) and n(θ), given by: (1/((n^2)(θ)))=(((cos(θ))^2)/(Subscript[n, 0]^2))+((sin(((theta))^2))/(Subscript[n, e]^2)). The first mode is called ordinary mode, and the second mode the extraordinary mode. The electric flux density, D, of the ordinary mode is normal to the plane defined by the optic axis - the z axis and the direction of propagation, k, and is perpendicular to the electric field, E. The extraordinary wave, however, has an electric flux density D that is normal to k, but E is not parallel to D because of its uniqe dependence on angle.

indexOfRefraction[theta_,nOrdinary_,nExtraordinary_]:=Sqrt[1/(((Cos[theta])^2)/nOrdinary+((Sin[theta])^2)/nExtraordinary)]

Snell’s Law

Light is always dispersed in a medium if it crosses from an index of refraction that is lower to an index of refraction that is higher. As we are modeling air to crystal interfaces, the index of refraction of the uniaxial crystal will always be greater than that of air. The energy of the electromagnetic wave must be conserved, so its frequency does not change. However, because the phase velocity within the material depends on the dispersion relation w(k) and is no longer the speed of light (it is actually slower than the speed of light), the wavelength of light changes. Birefringent materials do not only change the wavelength of light propagating through them but they also act as polarizers, since the ordinary and extraordinary modes are orthogonal to each other. Thus, light is also polarized along the two different indices of refraction present in birefringent materials. We can use the index of refraction to calculate the new wavelength propagating through the material (((Subscript[λ, 0])/n)) as dictated by the index of refraction for uniaxial crystals and Snell’s law (Subscript[n, 1])sin((Subscript[θ, 1]))=(Subscript[n, 2])sin((Subscript[θ, 2])).

thetaOrdinary[thetaIncident_,nOrdinary_]:=ArcSin[Sin[thetaIncident]/nOrdinary]

thetaExtraordinary[thetaIncident_,nOrdinary_,thetaOptical_,nExtraordinary_]:=thetaExtra/.(NSolve[Sin[thetaIncident]==indexOfRefraction[thetaOptical+thetaExtra,nOrdinary,nExtraordinary]*Sin[thetaExtra],thetaExtra,Reals]/.C[1]->0)[[1]]

lambda[lambdaNot_,n_]:=lambdaNot/n

yPositionOfLambda[theta_]:=1*Tan[theta]

colorLambda[wavelength_]:=Which[wavelength<475,Purple,475≤ wavelength<510,Blue,510≤ wavelength<570,Green,570≤wavelength<590,Yellow,590≤wavelength≤650,Orange,650≤ wavelength,Red]

"thetaExtraordinaryValues=Table[thetaExtraordinary[thetaIncident,1.5,Pi/3,2.7],{thetaIncident,-Pi/3,Pi/3,Pi/12}]"

flattenedExtraordinaryThetas={-2.3427641759479108`,-2.5860103642927736`,-2.3651358378114242`,-2.5863829897647`,-2.5250048842930477`,-2.6947000199979803`,-2.5269482276645903`,-2.696315902597508`,-2.7202123201553032`,-2.8229063326052244`,-2.7217439669923817`,-2.830377432475499`,-2.9250980023850413`,-2.9701562565159345`,-2.9316215188675185`,-2.9805796947321794`,0.`,0.`,0.`,0.`,0.21649465120475242`,0.17143639707385996`,0.20997113472227394`,0.16101295885761588`,0.42138033343448933`,0.3186863209845689`,0.41984868659741226`,0.31121522111429345`,0.6165877692967453`,0.4468926335918126`,0.6146444259252029`,0.44527675099228525`,0.7988284776418826`,0.5555822892970195`,0.7764568157783688`,0.5552096638250931`};

someOptions=Table[{thetaIncident,lambdaIncident,Graphics[{colorLambda[lambdaIncident],Line[{{0,0},{1,0}}],colorLambda[lambda[lambdaIncident,1.3]],Line[{{1,0},{2,yPositionOfLambda[thetaOrdinary[thetaIncident,1.3]]}}],colorLambda[lambda[lambdaIncident,indexOfRefraction[thetaExtraordinary[thetaIncident,1.3,Pi/3,1.5]+Pi/3,1.3,1.5]]],Line[{{1,0},{2,yPositionOfLambda[thetaExtraordinary[thetaIncident,1.3,Pi/3,1.5]]}}]}]},{thetaIncident,-Pi/3,Pi/3,Pi/12},{lambdaIncident,400,700,50}]

someOptions=Flatten[{{{-(π/3),400,},{-(π/3),450,},{-(π/3),500,},{-(π/3),550,},{-(π/3),600,},{-(π/3),650,},{-(π/3),700,}},{{-(π/4),400,},{-(π/4),450,},{-(π/4),500,},{-(π/4),550,},{-(π/4),600,},{-(π/4),650,},{-(π/4),700,}},{{-(π/6),400,},{-(π/6),450,},{-(π/6),500,},{-(π/6),550,},{-(π/6),600,},{-(π/6),650,},{-(π/6),700,}},{{-(π/12),400,},{-(π/12),450,},{-(π/12),500,},{-(π/12),550,},{-(π/12),600,},{-(π/12),650,},{-(π/12),700,}},{{0,400,},{0,450,},{0,500,},{0,550,},{0,600,},{0,650,},{0,700,}},{{(π/12),400,},{(π/12),450,},{(π/12),500,},{(π/12),550,},{(π/12),600,},{(π/12),650,},{(π/12),700,}},{{(π/6),400,},{(π/6),450,},{(π/6),500,},{(π/6),550,},{(π/6),600,},{(π/6),650,},{(π/6),700,}},{{(π/4),400,},{(π/4),450,},{(π/4),500,},{(π/4),550,},{(π/4),600,},{(π/4),650,},{(π/4),700,}},{{(π/3),400,},{(π/3),450,},{(π/3),500,},{(π/3),550,},{(π/3),600,},{(π/3),650,},{(π/3),700,}}},1]

Manipulate[someOptions[[i]][[3]],{i,28,63,1}]