Waves (comments to: moloney@nextwork.rose-hulman.edu) [*.ms = maple, *.xls = excel]

MAPLE text, input lines and worksheet downloads

In this animation, we imagine ourselves at x=0 with the 'beats' moving past us.

restart; with(plots):

c:=2;

omega:=c*k;

psi:=sin(k*x-omega*t);

y1:=subs(k=k1,psi); # one wavelength

y2:=subs(k=k2,psi); # another wavelength

ytotl:=y1+y2;

values:={k1=5,k2=45/10}; NN:=25; dtime:=6;

h:=animate(subs(values,ytotl), x=-2..8,t=0..dtime,color=red,frames=NN,numpoints=150):

g:=animate(subs(values,y1), x=-2..8,t=0..dtime,color=black,frames=NN,numpoints=150):

i:=animate(subs(values,y2), x=-2..8,t=0..dtime,color=blue, frames=NN,numpoints=150):

j:=display([h,g]): display([i,j]);

Download beats.mws [Use 'Save File']

'Steering a beam' in Multiple Source Interference

8/25/95 after BMWidth.ms Learning objectives: Translation of equation for
2D wave to off-origin source points Physical intuition regarding interference

Show a radiated beam being 'steered' by changing the phase shift between adjacent sources MAPLE: animate3d

Remarks: Takes some minutes to execute on most machines. Array of sources extends along the x-axis, centered at the origin

Wavelength, spacing, x, and y could be in meters k0 = wave number = 2 Pi/ wavelength spacing = distance between sources xmax, ymax used to set plot range dimensions

restart; with(plots):

wavelength:=1.0; k0:=628/(100*wavelength);

NrSources:=4; spacing:=6/10; Atotal:=0;

for j from 1 to NrSources do m:=(2*j-5)/2: # m = {-3/2 , -1/2, 1/2, 3/2}
for four sources

r:=((x-m*spacing)^2+y^2)^0.5: # distance from source to (x,y)

Atotal:=Atotal+cos(k0*r-m*phi): od:

Atotal:=A[1]+A[2]+...+A[NrSources]: # amplitude pattern at t=0 due to all four sources

Animate the effect of changing source separations

xmax:=4: ymax:=4: # Set plot area scale

The orientation is theta=0, phi =0 (phi = any), so the x-axis goes down
as we look from the z-axis

(Takes several minutes on most machines.)

animate3d(Atotal, x=-xmax..xmax, y=-ymax/5..9*ymax/5,phi=0/10..14/10, grid=[35,35],style=patchnogrid, shading=zgreyscale, scaling=constrained, orientation=[0,Pi/4], tickmarks=[0,0,0],frames=10);

Download bmsteer.mws [Use 'Save File']

**8/22/95 circvec.ms CW and CCW rotating vectors added together, animated.
Circular polarizations give linear polarization at x=0. But different wavelengths
give a phase shift when x>0 so that we get a rotation of the plane of
polarization.**

**restart;with(plots):
kval:=2*Pi/wavelength;
k[2]:=subs(wavelength=11/10,kval);
k[1]:=subs(wavelength=1,kval);**

**x:=1/80; clockwise rotating wave **

**Qx1:=5.0*cos(k[1]*x+3*t+phi[1]);
Qy1:=5.0*sin(k[1]*x+3*t+phi[1]); **

**counterclockwise rotating wave **

**Qx2:=5.0*cos(k[2]*x-3*t+phi[2]);
Qy2:=5.0*sin(k[2]*x-3*t+phi[2]);**

**NN:=20; phi[1]:=0; phi[2]:=0;**

**for i from 0 to NN do
t:=i/10;
g.(i):=plot([0,0,Qx1,Qy1],color='red',thickness=3):
h.(i):=plot([0,0,Qx2,Qy2],color='blue',thickness=3):
j.(i):=plot([0,0,Qx1+Qx2,Qy1+Qy2],thickness=5):
s.(i):=display({g.(i),h.(i),j.(i)}): od: **

**animation of CW, CCW, and their sum **

**display([s.(0..NN)],insequence=true); **

Download circvec.ms [Use 'Save File']

PHYSICS RESOURCE PACKETS PROJECT File: CIRCPOL.MS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 3

A ‘transverse' wave traveling along the z-axis may be represented as y(z,t) = A cos (kz - wt + phi). The frequency w of the wave normally stays the same as the wave passes through different media, but the velocity v = w/k changes, so k will be different in different media. A circularly-polarized wave is made up of two linearly polarized waves which are 90 degrees out of phase and also perpendicular to one another in space (e. g., one oscillates in x-dir, one oscillates in y-dir). a) Construct a circularly polarized wave whose displacement rotates counterclockwise in the x,y plane. Test the correctness of your wave by plotting at two successive times fairly close together, or by using the animate command in Maple, demonstrating the CCW rotation. b) Do the same for a clockwise-rotating circularly polarized wave, keeping the amplitude the same. Do a couple of plots to show the correctness of the function. c) Add the solutions together from parts a) and b). Do a couple of plots, or an animation, to show that the resultant wave is a ‘linearly' polarized wave.

restart; with(linalg): with(plots):

Part A: Counterclockwise First, create two linearly polarized waves (traveling along the z axis) which are 90 degrees out of phase.

WaveXz:=A*cos(k*z-omega*t);

WaveYz:=A*cos(k*z-omega*t+Pi/2);

k:=omega/v; # Put k in terms of omega and v

Put the traveling wave in vector form .

WaveXYz:=vector([WaveXz,WaveYz,z]);

Substitute in values for the constants so the function can be plotted. Make 11 plot of the wave for t=0..10.

WaveXYzsub:=subs({omega=1,v=1,A=1},op(WaveXYz));

for t from 0 to 10 do h.t:=spacecurve(WaveXYzsub,z=0..20,axes=normal);
od

t:='t'; # change back to a variable. Display the 11 plots in sequence.

display([h.(0..10)],insequence=true,title=`Counterclockwise Rotation`);

Part B: Clockwise First, create two linearly polarized waves (traveling along the z axis) which are 90 degrees out of phase.

t:='t';

WaveXz2:=A*cos(k*z-omega*t);

WaveYz2:=A*cos(k*z-omega*t-Pi/2);

Put the traveling wave in vector form.

WaveXYz2:=vector([WaveXz2,WaveYz2,z]);

Substitute in values for the constants so the function can be plotted. Make 11 plots of the wave for t=0..10.

WaveXYzsub2:=subs({omega=1,v=1,A=1},op(WaveXYz2));

for t from 0 to 10 do j.t:=spacecurve(WaveXYzsub2,z=0..20,axes=normal);
od:

t:='t'; # change back to a variable, display 11 plots in sequence

display([j.(0..10)],insequence=true,title=`Clockwise Rotation`);

Part C: Summing up parts A and B.

Add together WaveXYz and WaveXYz2 to show that the result is linear.

ZwaveSum:=evalm(WaveXYzsub+WaveXYzsub2);

Make 11 plot of the wave for t=0..10.

for t from 0 to 10 do g.t:=spacecurve(ZwaveSum,z=0..20,axes=normal);
od:

t:='t'; # change t back to a variable. Display the 11 plots in sequence.

display([g.(0..10)],insequence=true,title=`Linearly polarized wave`);

Download circpol.mws [Use 'Save File']

18 Aug 95 (phase velocity is twice the group velocity).

This is a problem of dispersion of waves in deep sea waters. Let us only work with two waves and observe the beat pattern. We should be able to see that the phase velocity of the individual waves is much greater than the group velocity. The dispersion relation is:omega = sqrt(k*g), where omega is the angular frequency, k the wave vector, and g the acceleration due to gravity.

Remarks: Illustrates animation with specific colors assigned.

restart; with(plots):

omega:=sqrt(g*k); # write down the dispersion relation

psi:=sin(k*x-omega*t); #general form of a wave function for deep sea waves

y[1]:=subs(k=k1,psi); y[2]:=subs(k=k2,psi);

vp:=omega/k; # phase velocity

vg:=diff(omega,k); # group velocity

ratio:=vp/vg; # phase velocity over group velocity

evalf(subs(k=10,g=98/10,ratio));

y_total:=y[1]+y[2]; values:={k1=5,k2=9/2,g=98/10};

h:=animate(subs(values,y_total),x=-2..8,t=0..3,color=red, frames=30,numpoints=100):

g:=animate(subs(values,y[1]), x=-2..8,t=0..3,color=black,frames=30,numpoints=100):

i:=animate(subs(values,y[2]), x=-2..8,t=0..3,color=blue, frames=30,numpoints=100):

j:=display([h,g]): display([i,j]);

Notice that the individual peaks move ahead of the 'group' peak. By single stepping through the animation you can see a trailing peak come into the picture in the second or third frame, and overtake the big peak by the end of the animation.

Download deepsea.mws [Use 'Save File']

20 Aug 95 student 'starter' version

PROBLEM STATEMENT. A source of radio waves in the sky is detected by a line of N equally-spaced detectors along the y-axis. When the source has a positive y-coordinate, the detector with the greatest y-value will be closest to the source. Signals from the source will reach this detector (detector # 1) first. Signals reaching the second detector will be out of phase with those in detector #1, and so on with the other detectors. The incoming waves look like y = A cos (kL -wt). The phase of each wave is the argument of cosine, namely phase = kL-wt. At any time t, every detector will have the same k, w, and t value, but the L's are all different, so there is a phase difference between the signals from the source. To get a big signal from the array, all detector signals must be in phase, so a phase delay D is introduced between detectors. When D is set correctly, it compensates for the phase difference due to different L values, and brings the signals back into phase. [D is subtracted from the phase of dectector #2, 2D is subtracted from the phase of detector #3, and so forth.]

In this problem, we calculate the resultant of all 'phasors' from the detectors, and plot its square vs. the phase delay D. The magnitude of the resultant is the amplitude of the total detector response, something we do not prove in this problem.

This (detect0.ms) is a 'starter' program, intended for students to learn from and build on. The source is initially placed at theta=0 (it is on the x-axis). The array output will be a maximum when the phase delay D is zero, because waves from the source are all arriving in phase. Some 'aliasing' can be seen in the detector output.

The radio source in this problem is located at a distance of 10^9 m from the origin, and its wavelength is 0.21 m. The detector spacing is 0.10 m, and initially there are N=7 detectors. Later in the problem A is set equal to 1.

Note: In 'do' loops, be sure to close the loop with the 'od' statement. Failing to do this often makes the program stop executing, and you have to close and then re-enter Maple.

PARTS OF THE PROBLEM :

a) Initial run-through. Run through the program detect0.ms, with N=7 detectors in the array. In the plot of detector response vs phase delay D, notice that the big response comes at D=0. Record the delay value when Rsquared first equals zero.

b) Effect of array size on ability to locate a source. Increase the number of detectors N to equal 15. Rerun and replot s1. Make a note of the D value at which the detector response first becomes zero. This value should be less than when you had 7 detector elements. Your larger array should be able to locate a source more precisely than before.

c) Path differences between detectors are quite constant. Now move the source in the sky to theta=0.600 radians, label its detector response s1. Substitute numerical values into the distances L[i], and call these LVal[i]. They are distances from the new source to the detector elements. Create DeltaL[i], which are differences between LVal[i], and display them. They should be rather constant. Create DDeltaL[i], which are differences between DeltaL[i], and display them. They should be tiny compared to a wavelength. Report the sizes of the path differences. Compare this to W sin theta, where W is the spacing between detectors.

restart; with(plots): assume(k>0,wavelength>0,t>0);

The number of digits is set to 14 because of the 10^9 m distance used later in the problem. It turns out you need at least 4 or 5 more digits available so that the calculated values of kL are accurate. There is a compromise between accuracy and speed. Using more digits improves accuracy, but calculations (and plots) can take a lot longer.

Digits:=14; k:=2*Pi/wavelength;

X:=R*cos(Theta); Y:=R*sin(Theta);

N:=7; m:=trunc(N/2);

Determine the distances L[i] from the source to each detector element

for i from 1 to N do L[i]:=sqrt( (Y-spacing*(N-i-m))^2
+ X^2) od:

R_x:=0: R_y:=0:

Let t=0 and get magnitude of resultant; first obtain x and y components of the resultant, then square.

for i from 1 to N do R_x:=R_x+A*cos(k*L[i]-D*(i-1)):
R_y:=R_y+A*sin(k*L[i]-D*(i-1)): od:

Rsquared:=R_x^2+R_y^2:

values:={R=1*10^9,A=1,spacing=0.10,wavelength=0.21};

Substitute problem values into the resultant vector squared, except for theta value

g:=subs(values,Rsquared):

Now create the response function from source #1 in the sky at theta=0

s1:=subs(Theta=0.000,g): plot(s1,D=0..2*Pi,numpoints=200);

Download detect0.mws [Use 'Save File']

20 Aug 95 faculty version, with extra stuff done [student 'starter' version is detect0.ms]

PROBLEM STATEMENT. A source of radio waves in the sky is detected by a line of N equally-spaced detectors along the y-axis. When the source has a positive x-coordinate, the detector with the greatest y-value will be closest to the source. Signals from the source will reach this detector (detector # 1) first. Signals reaching the second detector will be out of phase with those in detector #1, and so on with the other detectors. The incoming waves look like y = A cos (kL -wt). The phase of each wave is the argument of cosine, namely phase = kL-wt. At any time t, every detector will have the same k, w, and t value, but the L's are all different, so there is a phase difference between the signals from the source. To get a big signal from the array, all detector signals must be in phase, so a phase delay D is introduced between detectors. When D is set correctly, it compensates for the phase difference due to different L values, and brings the signals back into phase. [D is subtracted from the phase of dectector #2, 2D is subtracted from the phase of detector #3, and so forth.]

In this problem, we calculate the resultant of all 'phasors' from the detectors, and plot its square vs. the phase delay D. The magnitude of the resultant is the amplitude of the total detector response, something we do not prove in this problem.

detect0.ms is a 'starter' program, intended for students to learn from and build on. The source is initially placed at theta=0 (it is on the x-axis). The array output will be a maximum when the phase delay D is zero, because waves from the source are all arriving in phase.

The radio source in this problem is located at a distance of 10^9 m from the origin, and its wavelength is 0.21 m. The detector spacing is 0.10 m, and initially there are N=7 detectors. Later in the problem A is set equal to 1.

Note: In 'do' loops, be sure to close the loop with the 'od' statement. Failing to do this often makes the program stop executing, and you have to close and then re-enter Maple.

PARTS OF THE PROBLEM :

a) Initial run-through. Run through the program detect0.ms, with N=7 detectors in the array. In the plot of detector response vs phase delay D, notice that the big response comes at D=0. Record the delay value when Rsquared first equals zero.

b) Effect of array size on ability to locate a source. Increase the number of detectors N to equal 15. Rerun and replot s1. Make a note of the D value at which the detector response first becomes zero. This value should be less than when you had 7 detector elements. Your larger array should be able to locate a source more precisely than before.

c) Path differences between detectors are quite constant. Now move the source in the sky to theta=0.600 radians, label its detector response s1. Substitute numerical values into the distances L[i], and call these LVal[i]. They are distances from the new source to the detector elements. Create DeltaL[i], which are differences between LVal[i], and display them. They should be rather constant. Create DDeltaL[i], which are differences between DeltaL[i], and display them. They should be tiny compared to a wavelength. Report the sizes of the path differences. Compare this to W sin theta, where W is the spacing between detectors.

restart; assume(k>0,wavelength>0,t>0);

The number of digits is set to 14 because of the 10^9 m distance used later in the problem. It turns out you need at least 4 or 5 more digits available so that the calculated values of kL are accurate. There is a compromise between accuracy and speed. Using more digits improves accuracy, but calculations (and plots) can take a lot longer.

Digits:=14; k:=2*Pi/wavelength:

X:=R*cos(Theta): Y:=R*sin(Theta):

N:=7: m:=trunc(N/2):

DistBtwDets:=1/10: # distance between detectors

Determine the distances L[i] from the source to each detector element

for i from 1 to N do L[i]:=sqrt( (Y-spacing*(N-i-m))^2
+ X^2) od:

R_x:=0: R_y:=0:

Let t=0 and get magnitude of resultant; first obtain x and y components of the resultant, then square.

for i from 1 to N do

R_x:=R_x+A*cos(k*L[i]-D*(i-1)):

R_y:=R_y+A*sin(k*L[i]-D*(i-1)): od:

Rsquared:=R_x^2+R_y^2:

values:={R=1*10^9,A=1,spacing=DistBtwDets,wavelength=0.21};

Substitute problem values into the resultant vector squared, except for theta value

g:=subs(values,Rsquared):

Now create the response function from source #1 in the sky

s1:=subs(Theta=0.600,g):

plot(s1,D=0..2*Pi,numpoints=200);

for i from 1 to N do LVal[i]:=evalf(subs(values,Theta=0.600,L[i])); od:

for i from 1 to N-1 do DeltaL[i]:=LVal[i+1]-LVal[i]; od;

Wsint:=DistBtwDets*sin(0.600); # spacing sin(theta) , should be very close
to Delta-L values.

Remarks: Replace colons by semicolons to see actual values. The value of W sin theta is extremely close to the path differences calculated.

Download detect1.mws [Use 'Save File']

PHYSICS RESOURCE PACKETS PROJECT File: DOPPL.MS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 3

Assume a source radiating radially and draw 10 crests (circles) of the sound waves. Show the position of the wave crests for an object moving near the speed of sound, at the speed of sound, and faster than the speed of sound.

restart; with(plots):

Define the function "draw" so that it will return A for A>0 and 0 for A<0.

draw:=A->A*Heaviside(A);

The loop below will generate 11 parametric equations (h0,h1,h2,...,h10). h0 describes the motion of the first sound wave emitted, h1 describes the motion of the second wave and so on. The parametric equations are of the form [x(n),y(n),n=a..b]. We want to generate circles that radiate outward from the sound source at the speed of sound, so there is a cosine term in x(n), and a sine term in y(n). 343*j*(0.75) is the location of the center of the circle being drawn. In the animate command all 11 parametric equations are ploted on the same graph at the same time for t=0 to 11 seconds. The object moves along at a speed of 343*0.75 m/s and (for clarity) emits only one wavefront every second. The "draw" function controls the number and size of the circles drawn. When 0>t>1, the only time (t-j) is positive in the draw function is when j=0, so only the contents of h0 are visible in the animation. When 1>t>2, h0 is still drawn, but now the circle in h1 is also visible. The radius of the first circle is 343*(t-0), and the radius of the second circle is 343*(t-1). More of the circles will become visible as time progresses until they all are when t>10.

for j from 0 to 10 do h.j:=[343*draw(t-j)*cos(theta)+343*j*(0.75),343*draw(t-j)*sin(theta),theta=0..2*Pi]; od;

Animate h0..h10 for t=0..11 s.

animate({h.(0..10)},t=0..11,title=`Mach 0.75`,color=red);

Now increase the speed of the source from 343*0.75 m/s to 343 m/s.

for j from 0 to 10 do h.j:=[343*draw(t-j)*cos(theta)+343*j*(1),343*draw(t-j)*sin(theta),theta=0..2*Pi]; od: animate({h.(0..10)},t=0..11,title=`Mach 1`,color=red);

Increase the speed to 2*343 m/s. (Mach 2)

for j from 0 to 10 do h.j:=[343*draw(t-j)*cos(theta)+343*j*(2),343*draw(t-j)*sin(theta),theta=0..2*Pi]; od: animate({h.(0..10)},t=0..11,title=`Mach 2`,color=red);

Accelerate the source from Mach 0.75 to Mach 1.25 in 10 seconds.

for j from 0 to 10 do h.j:=[343*draw(t-j)*cos(theta)+343*j*(0.75+j/20),343*draw(t-j)*sin(theta),theta=0..2*Pi]; od: animate({h.(0..10)},t=0..11,title=`Accelerating from Mach 0.75 to Mach 1.25`,color=red);

Instead of animating the wavefronts as a function of time,animate them as a function of speed.

for j from 1 to 11 do h.j:=[343*j*cos(theta)-343*j*(speed),343*j*sin(theta),theta=0..2*Pi]; od: animate({h.(1..11)},speed=0..2,title=`Shape of wavefronts from Mach 0 to Mach 2`,color=red);

PHYSICS RESOURCE PACKETS PROJECT File: FUZZYTV.MS

Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby

http://www.rose-hulman.edu/~moloney Software: Maple V Release 3

TV sets receiving signals out of the air without benefit of cable are sometimes subject to periodic jumping of the picture when an airplane flies through or near the line-of-sight to the TV station. The line of sight distance from station to the TV set is 10 miles, and the wavelength is 9.5 feet. The TV station is at the origin (0,0,0) and the TV set is at (10,0,0) mi. An airplane flies perpendicular to this line at an altitude of 1/2 mile at a speed of 150 mi/hr. At t=0, the plane is closest to the TV set. At that time, its coordinates are (9.0, 0, 0.5) mi. Find the first 3 times after t=0 when waves bouncing off the plane will be out of phase with waves arriving directly at the set from the station. (At this point, the TV picture should be all messed up.)

restart; with(linalg): with(student): assume(t>0);

plane:=vector([9,150*t,1/2]); # The location of the plane
in miles. t is in hours.

station:=vector([0,0,0]); #
The location of the TV station (miles)

TVset:=vector([10,0,0]); #
The location of the TV set (miles)

dist1:=distance(station,plane); #
Distance from the station to the plane.

dist2:=distance(TVset,plane);

distreflect:=dist1+dist2; #
Total path length of reflected wave

distdirect:=10; #
Distance from the TV station direct to the TV set

difference:=distreflect-distdirect; # Path
difference between the direct path and the reflected path

Since a phase shift of 180 deg occurs when the waves are reflected from the plane, distructive interference will occur when the path difference is a multiple of 9.5 feet. Convert 9.5 feet into miles to solve the equation for time.

sol:=solve(difference=n*(95/10)/5280,t);

Pick the positive value for time.

Time:=sol[1];

Find the path difference in feet at t=0

dist0:=evalf(subs(t=0,difference))*5280;

Find the value of n at t=0.

dist0/9.5;

destructive interference would first occur when n = 74.

Find the first three times in seconds when the reflected signal is 180 degrees out of phase with the direct signal

for n from 74 to 76 do time.(n-73)=evalf(Time*3600); od;

The 'fuzzy' signals are occurring every 0.7 seconds or so

8/19/95 All plots in nice shape. Scale always 0..2*Pi.

Interference pattern for equispaced two, three, four, five, fifteen and twenty-five point sources. Shows the narrowing of interference peaks with more sources at constant source spacing.

restart; assume(y,real);

Definition of the different terms in the interference pattern: a (n) stands for the amplitude of the sources, d is the, spacing between the sources and is expressed in micrometers, L is the wavelength of light used (micrometers), Z is the distance of observation, and y is the corordinate that is perpendicular to the propagation of light.

d := 6.33; L := 0.633; Z := 10; t:= y;

Definition of the amplitude of the light scattered from the point sources at a plane perpendicular to the propagation of light. The point source at the center is a(1). The even numbered source lie above a(1) and the odd numbered sources lie below a(1). The distances of the even numbered sources are positive and the odd numbered sources are negative. For example source a(2) above a(1) is at a distance 'd' and the source a(3) below a(1) is at '-d'. Similarly source a(4) is at a distance '2d' and source a(5) is at '-2d' and so on.

Amplitude 15 point sources below source a(1):

for n from 1 to 15 do a(2*n-1) := 2*exp(-I* (n-1)*(Pi*d*t)/(L*Z)); od;

Amplitude of the 15 point sources above source a(1):

for n from 1 to 15 do a(2*n) := 2*exp(I* (n)*(Pi*d*t)/(L*Z)); od;

In the Young's double slit interferometer we have two point sources. Therefore the total amplitude of sources a(1) and a(2) or a(3) is calculated.

f3:=(abs(sum(a(n1), n1 = 1..2))); f4:=simplify(f3^2);

plot(f4, y = 0..6,title=`TwoSources`);

Three point sources we take sources s(1), s(2) and s(3).

f5:= abs(sum(a(n1), n1 = 1..3))^2; simplify(f5); plot(", y = 0..6);

Four,five,15, and 25 point sources

f6:= abs(sum(a(n1), n1 = 1..4))^2; simplify(f6);

plot(", y = 0..6,title=`Four Sources`);

f7:= sum(a(n1), n1 = 1..5)^2; f8:=simplify(f7); plot(f8, y = 0..6);

f9:= sum(a(n1), n1 = 1..15)^2; f10:=simplify(f9); plot(f10, y = 0..6, numpoints
= 300);

f11:= abs(sum(a(n1), n1 = 1..25))^2; simplify(f11); plot(", y = 0..6,numpoints
=500);

Download interfer.mws [Use 'Save File']

PHYSICS RESOURCE PACKETS PROJECT File: LICENSE.MS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 3

This problem has light rays coming parallel to the axis and refracting inside a sphere then coming back out. The idea is to show that for a range of index of refraction one will get rays coming back out which are mostly parallel to the axis, so making a spherical bead a good approximation to a retro-reflector. Small beads of this type are sold in 3M tape to go on license plates and other applications.

Perry shows n=1.85 below, but one would want to have a way to compare n=1.85 to other indices of refraction for effectiveness.

restart; with(plots):

n:=1.85; # index of refraction

for i from 0.05 to 1 by 0.05 do # Position of horizontal lines

L1[i]:=i; # Equation of incoming ray.

Find where incoming ray hits the circle.

sol[i]:=fsolve({y1[i]=L1[i],x1[i]=cos(theta1[i]),y1[i]=sin(theta1[i])},{theta1[i],x1[i],y1[i]},{theta1[i]=-Pi/2..Pi/2});

assign(sol[i]);

Find the angle of the refracted ray.

snell[i]:=1*sin(theta1[i])=n*sin(theta2[i]);

theta2[i]:=solve(snell[i],theta2[i]);

Find the slope angle of the refracted ray.

thetaL2[i]:=theta1[i]-theta2[i];

L2[i]:=tan(thetaL2[i])*(x-x1[i])+y1[i]; # Equation of the line for the
refracted ray.

Find where the refracted ray intersects the circle.

sol[i]:=fsolve({y2[i]=subs(x=x2[i],L2[i]),x2[i]=cos(theta3[i]),y2[i]=sin(theta3[i])},{theta3[i],x2[i],y2[i]},{theta3[i]=Pi/2..3*Pi/2});
assign(sol[i]); # make proper assignments of parts of the solution

thetaL3[i]:=thetaL2[i]-2*theta2[i]; # Find the slope angle of the reflected
ray.

L3[i]:=tan(thetaL3[i])*(x-x2[i])+y2[i]; # Equation of the line for the
reflected ray.

Find where the reflected ray intersects the circle.

sol[i]:=fsolve({y3[i]=subs(x=x3[i],L3[i]),x3[i]=cos(theta4[i]),y3[i]=sin(theta4[i])},{theta4[i],x3[i],y3[i]},{theta4[i]=Pi..2*Pi}); assign(sol[i]);

thetaL4[i]:=theta4[i]+theta1[i]; # Find the slope
angle of the reflected ray leaving the circle.

L4[i]:=tan(thetaL4[i])*(x-x3[i])+y3[i]; # Equation of the line for the
ray leaving the circle.

od:

Make plots for the incomming and outgoing rays.

for i from 0.05 to 1 by 0.05 do

f.(trunc(i*20)):=plot({[x,L1[i],x=x1[i]..2],[x,L2[i],x=x2[i]..x1[i]]},color=blue,scaling=constrained):
h.(trunc(i*20)):=plot({[x,L3[i],x=x2[i]..x3[i]],[x,L4[i],x=x3[i]..2]},color=green,scaling=constrained):
od:

g:=plot([cos(theta),sin(theta),theta=0..2*Pi],color=black): # Make a plot of the circle.

Display all the plots on the same graph.

display({g,f.(1..20),h.(1..20)});

Download license.mws [Use 'Save File']

PHYSICS RESOURCE PACKETS PROJECT File: MIRAGESM.MS

Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby

http://www.rose-hulman.edu/~moloney Software: Maple V Release 3

**Background and Theory.** A mirage will occur when surface air is
hotter than the air above it. Then a flat road surface will appear to give
glassy reflections of objects near the horizon. The index of refraction
increases with height, so we will say

n(y) = a + by = a (1+e y),

where a, b, and e = b/a are positive constants. Snell's law, n sin(theta) = constant, applies to all rays.

Call the point (x0,y0) the point where a ray is traveling parallel to the ground, and the angle theta in Snell's law is 90 degrees. In Snell's law anywhere on the path, sin(theta) = dx (dx^2 + dy^2)^(-1/2) = (1+(dy/dx)^2)^(-1/2). y0 is the lowest point on the path; all other y values will be greater. Snell's law now reads

n sin(theta) = a (1+ey)(1+(dy/dx)^2)^(-1/2) = constant = n(y0) sin(Pi/2) = n(y0) = a(1+ey0)

Rearranging this gives (y + 1/e)/(y0 + 1/e) = (1+(dy/dx)^2)^(1/2). This equation has a straightforward solution, when we make the substitution (y + 1/e)/(y0 + 1/e) = cosh(z), which leads to dy/dx = sinh(z), and finally to x=(y0 + 1/e) z. The full equation of a ray is then:

y + 1/e=(y0 + 1/e) cosh((x-x0)/(y0+1/e)), where (x0,y0) is the lowest point on the path.

We wish to launch several rays at different angles from a fixed point (xp,yp) and plot their paths. From Snell's law we find a different minimum height y0p for each ray: n(yp) sin(theta-p) = n(y0p), or (1+e yp) sin(theta-p) = 1+e y0p, or y0p=(1/e)((1+e yp) sin(theta-p)-1). Next, we find x0p from yp+1/e = (y0p+1/e) cosh((xp-x0p)/(y0p+1/e)), or 1/sin(theta-p) = cosh((xp-x0p)/(yp+1/e)sin(theta-p)).

Now we are ready to plot the path starting from (xp,yp).

a) Let H = 20 m, a = 1.003, b = 1.14607*10^(-5), and A = 1.5 degrees. Plot the path of a ray launched under these conditions. b) Find the location along the x-axis where the marage would be visible for a viewing height of 1.3 m to 1.4 m.

Part A: Path of Ray

Enter in the given parameters. AlphaP and ThetaP must be in radians. Don't put in the value for AlphaP for now so that it remains a variable in the final equation.

restart; with(plots):

xp:=0; yp:=30; a:=1.003; b:=1.14607*10^(-5);

e:=b/a; AlphaP:=(1.5*Pi/180); Thetap:=(Pi/2-AlphaP);

Now determine yop from the formula derived in the problem statement.

yop:=evalf((((1+e*yp)*sin(Thetap))-1)/e,20);

Now determine xop from the formula derived in the problem statement.

xop1:=yp+1/e=(yop+1/e)*cosh((xp-xop)/(yop+1/e));

Solve the above equation for xop.

xop2:=solve(xop1,xop);

Now enter in the formula for the path of the ray.

RayPath1:=y+1/e=(y0+1/e)*cosh((x-x0)/(y0+1/e));

Enter in the initial conditions for x0 and y0.

RayPath2:=subs({y0=yop,x0=xop2},RayPath1);

Solve RayPath2 for y.

RayPathy:=solve(RayPath2,y);

Now plot the path of the ray.

plot(RayPathy,x=0..3000,title=`Path of the Ray`);

Part B: Where is the mirage visible?

x1:=solve(RayPathy=1.4,x); # Find when the ray
is first visible.

x2:=solve(RayPathy=1.3,x); # Find when the ray is no longer visible.

Thus, the mirage is visible when the observer is 2768.5 m to 2786.4 m away from the object.

Download miragesm.mws [Use 'Save File']

PHYSICS RESOURCE PACKETS PROJECT File: MIRAGECM.MS Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney Software: Maple V Release 3

The problem statement is similar to miragesm.ms, except in this more complicated vesion, we will locate the position of the image in the mirage.

We wish to launch several rays at different angles from a fixed point (xp,yp) and plot their paths. From Snell's law we find a different minimum height y0p for each ray: n(yp) sin(theta-p) = n(y0p), or (1+e yp) sin(theta-p) = 1+e y0p, or y0p=(1/e)((1+e yp) sin(theta-p)-1). Next, we find xop from yp+1/e = (yop+1/e) cosh((xp-x0p)/(y0p+1/e)), or 1/sin(theta-p) = cosh((xp-x0p)/(yp+1/e)sin(theta-p)).

Now we are ready to plot the path starting from (xp,yp).

a) Let H = 20 m, a = 1.003, b = 1.14607*10^(-5), and A = 1.5 degrees. Plot the path of a ray launched under these conditions. Repeat for rays launched at 1.48 degrees and 1.49 degrees. b) Estimate the location of the image from you plot. This will involve drawing lines by hand on the plot, or using Maple or Mathematica to find tangent lines to the rays.

Part A: Path of Ray

Enter in the given parameters. AlphaP and ThetaP must be in radians. Don't put in the value for AlphaP for now so that it remains a variable in the final equation.

restart; with(plots): Digits:=20;

xp:=0; yp:=30; a:=1.003; b:=1.14607*10^(-5);

e:=b/a; AlphaP:=1.5*Pi/180; Thetap:=Pi/2-AlphaP;

Now determine yop from the formula derived in the problem statement.

yop:=evalf((((1+e*yp)*sin(Thetap))-1)/e);

Now determine xop from the formula derived in the problem statement.

xop1:=yp+1/e=(yop+1/e)*cosh((xp-xop)/(yop+1/e));

Solve the above equation for xop.

xop2:=solve(xop1,xop);

RayPath1:=y+1/e=(y0+1/e)*cosh((x-x0)/(y0+1/e));
# Now enter in the formula for the path of the ray.

RayPath2:=subs({y0=yop,x0=xop2},RayPath1); #
Enter in the initial conditions for x0 and y0 .

RayPathy:=solve(RayPath2,y); # Solve RayPath2
for y.

plot(RayPathy,x=0..3000,title=`Path of the Ray`);
# Now plot the path of the ray.

Part B: Vary Alpha

Now find the path of the ray when AlphaP is 1.49 degrees. All commands are the same, so only the final path length is printed.

AlphaP:=evalf(1.49*Pi/180): Thetap:=evalf(Pi/2-AlphaP):

yop:=(((1+e*yp)*sin(Thetap))-1)/e:

xop1:=yp+1/e=(yop+1/e)*cosh((xp-xop)/(yop+1/e));

xop2:=solve(xop1,xop);

RayPath1:=y+1/e=(y0+1/e)*cosh((x-x0)/(y0+1/e)):

RayPath2:=subs({y0=yop,x0=xop2},RayPath1):

RayPathya:=solve(RayPath2,y);

Now find the path of the ray for when AlphaP is 1.48 degrees.

AlphaP:=evalf(1.48*Pi/180): Thetap:=evalf(Pi/2-AlphaP):

yop:=(((1+e*yp)*sin(Thetap))-1)/e:

xop1:=yp+1/e=(yop+1/e)*cosh((xp-xop)/(yop+1/e));

xop2:=solve(xop1,xop);

RayPath1:=y+1/e=(y0+1/e)*cosh((x-x0)/(y0+1/e)):

RayPath2:=subs({y0=yop,x0=xop2},RayPath1):

RayPathyb:=solve(RayPath2,y);

Now plot the path of the three rays for AlphaP=1.5,1.49,1.48

plot({RayPathy,RayPathya,RayPathyb},x=2000..2800,title=`Rays for AlphaP=1.5,1.49,1.48`);

Part C: Determination of image location

Assume that the observer is at x=2770 and that the observer can see
all rays between 1.3 and 2.3 m high.

Where does the image appear to be located?

In order to determine the image's location, the tangent lines to the three paths must be calculated and then simultaneously solved.

At x=2770, the corresponding y-values are:

yvalue:=evalf(subs(x=2770,RayPathy));

yvaluea:=evalf(subs(x=2770,RayPathya));

yvalueb:=evalf(subs(x=2770,RayPathyb));

The slope at the points calculated is the derivative of the paths at the calculated point.

slope:=evalf(subs(x=2770,diff(RayPathy,x)));

slopea:=evalf(subs(x=2770,diff(RayPathya,x)));

slopeb:=evalf(subs(x=2770,diff(RayPathyb,x)));

The tangent lines at the points, in point-slope form, are:

line:=(y-yvalue)=slope*(x-2770);

linea:=(y-yvaluea)=slopea*(x-2770); l

ineb:=(y-yvalueb)=slopeb*(x-2770);

y1:=solve(line,y); y2:=solve(linea,y); y3:=solve(lineb,y);
# solve each equation for y

sol1:=solve({line,linea},{x,y}); # Find
the intersection of lines y1 and y2.

sol2:=solve({line,lineb},{x,y}); # Find the
intersection of lines y1 and y3

The three rays do not converge on the same point which indicates the the fuzzy appearance most mirages. Howerer, the rays are relatively close together, which allows the object to be seen.

Thus the image appears to be comming from approximately (-1.4,-13.846)

plot({y1,y2,y3,RayPathy,RayPathya,RayPathyb},x=-10..2800,y=-14..10,color=black,title=`Path of rays`);

Download miragecm.mws [Use 'Save File']