Fancy3.mws

**vector stuff**

**op**

**electric field plot**

**envelope detector**

**driven oscillator**

E&M Packet Problem ABW solution Potential and field plots of three point charges Learning objectives: Translation of standard 1/r potential to off-origin points Physical intuition regarding potential Identify equilibrium points Relate potential plot and electric field vectors MAPLE: plot3d, grad, fieldplot commands |

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

arbitrary units are used throughout in Coulomb's law

**Q1:=1.;Q2:=1;Q3:=-1:**
# define relative magnitude of three charges

Set plot range; use non-even values to avoid singular points in grids

**xmax:=0.755; ymax:=0.755: **#plot domain
-xmax<x<xmax and-ymax<y<ymax

**coul1:= Q1/((x-0.5)^2+y^2)^0.5;** # potential due to 1st charge only

**coul2:=Q2/((x+.5)^2+y^2)^0.5; **# potential due to 2nd charge only

**coul3:=Q3/((x)^2+(y-.7)^2)^0.5; **# potential due to 3rd charge only

**coul123:=coul1+coul2+coul3; **# potential due to all three charges

3d countour plot of potential function

**plot3d(coul123, x=-xmax..xmax, y=-ymax/2..2*ymax,
grid=[49,49],style=patchcontour,contours=10,tickmarks=[1,1,100],view=-20..20);**

Get electric field from negative gradient of potential

(Note: "-grad(U)" doesn't seem to work, so using "grad(-U)"

**vect:=grad(-coul123, vector([x,y]));**

Create new vector field of unit vectors E/|E|

**vect2:=evalm(vect/sqrt(vect[1]^2+vect[2]^2));
fieldplot([vect2[1],vect2[2]],x=-xmax..xmax,y=-ymax..ymax+1/2,title=`unit
vectors showing direct'n of E field `);**

Compress dynamic range of E-vector lengths by taking ln(1+|E|)

log is taken twice to put all lengths within range of fieldplot command

**vect3:=evalm(vect2*ln(ln(sqrt(vect[1]^2+vect[2]^2)+1)+1));
fieldplot([vect3[1],vect3[2]],x=-xmax..xmax,y=-ymax..ymax+1/2,title=`electric
field vectors`);**

*Plot the phase and amplitude response of a driven,
damped oscillator as a function of frequency.*

This problem illustrates dealing with complex numbers, and the use of the 'op' command to parse an output expression.

It is interesting that very heavy damping (b=10, k=10, m=1) crushes the response away from zero frequency.

**restart;with(plots):
assume(k,real);assume(m,real);assume(A,real);
assume(omega>0);assume(t,real);assume(b,real);assume(phi,real);**

**eqn:=A*exp(I*(omega*t+phi))-k*x(t)-b*diff(x(t),t) = m*diff(x(t),t,t);
s:=dsolve(eqn,x(t));** # (takes a little while)

**x:=rhs(s);
values:={A=1,m=1,k=10,b=1/10};
xs:=subs(values,x);**

**steady:=op(1,xs);** # keep only the oscillating terms for steady
state response

**ev:=evalc(subs(t=0,phi=0,steady));**

Phase of the response as a function of frequency. The conventional definition is the phase of the driver minus the phase of the driven system. The relative amplitude response will be taken as exp(-I alph).

**assume(alph,real);
cosal:=op(1,ev); # cos(alph)
sinal:=-op(2,ev)/I; # sin(alph)
alph:=arctan(sinal/cosal);
plot(alph,omega=0..10,title=`response phase vs. omega`);**

The driver is in phase with the driven system at very low frequencies, and pulls ahead until it is Pi/2 ahead at resonance, and goes on to be Pi ahead at very high frequencies.

Amplitude response as a function of frequency

**xcc:=subs(I=-I,steady);** # get complex conjugate

**xmag:=simplify(steady*xcc);** # get magnitude squared

**xev:=simplify(evalc(xmag));
plot(sqrt(xev),omega=0..10,title=`Response amplitude vs freq.`);**

In AM (amplitude-modulated) radio, the amplitude of a high-frequency
signal is modulated by the lower-frequency audio information. AM radios
use '**envelope detectors' involving RC circuits** in order to remove
the high-frequency (the 'carrier wave') and leave the audio signal.

This problem gives you a simple modulated signal. It also gives you the 'rectified' signal we would see after passing through a diode (which lets current flow in one direction only.).

You are to pass the signal through a series RC circuit, first through a resistor R and then a capacitor C to ground. The 'envelope' signal is to be the voltage across the capacitor.

Letting v be the incoming signal, the equation to be solved is v -IR -q/C = 0 .

You must replace the current I by a function of q, the charge on the capacitor. This (call it eq1) will be the generic equation to be solved for q the charge on the capacitor as a function of time. By substituting values of R and C into eq1, you get eq2. This new equation is to be solved numerically for q (the charge on the capacitor) as a function of time.

a) Get the rectified signal first, by using the Heaviside function. The input signal is

v := sin(96 Pi t) *sin(12000 Pi t).

b) Set up the differential equation, using q(t) for the charge on the capacitor C . C is connected to ground and to the resistor R . R has the input voltage on one side and C on the other. [The current will be dq/dt]

c) Substitute values into the DE to make a new DE which will be solved numerically. Start off by setting RC equal to 1/10 the period of the 'carrier frequency'.

d) Solve the new DE numerically. Plot (use odeplot) and see what this gives you for an envelope function Compare this to the plot of the modulated function.

**restart; with(plots):
th:=6000*2*Pi*t; # c**arrier frequency is
6000 hz

**plot(v,t=0..0.015,numpoints=400, title=`plot
of modulated signal`); **

From here on down, the students should be doing the work.** **

**eq1:=v-R*diff(q(t),t)-q(t)/C=0;
eq2:=subs(R=45,C=1/100000,eq1); **

RC is 0.45 ms. carrier pd= 0.17 ms. modulation pd .= 20 ms.

**sol:=dsolve({eq2,q(0)=0}, q(t),numeric); odeplot(sol,[t,q(t)*100000/1],0..15/1000,numpoints=400,title=`signal
across capacitor`);**

Note that we don't get all of the signal back, nor is the shape a perfect fit to the envelope.

End of Fancy3.mws