DECLARE FUNCTION ATN2! (x!, y!) DECLARE SUB Pager () pi = 4! * ATN(1!) 10 PRINT CHR$(12); 20 PRINT "Forced Harmonic Oscillator:" 22 PRINT "Do you want to go directly to the calculations? (Y or N)"; 23 INPUT ANS$ 26 IF ANS$ = "Y" OR ANS$ = "y" OR ANS$ = "YES" OR ANS$ = "yes" THEN GOTO 1200 30 PRINT 40 PRINT "We will study the motion of a FORCED harmonic Oscillator." 50 PRINT "It is a 1D system. The displacement from equilbrium is called" 60 PRINT "x. It may be a translation or a rotation or a voltage, etc." 70 PRINT "If we select to use the terminology of a spring attached to" 80 PRINT "a mass m in a resistive media where the damping force, F(v)" 90 PRINT "is -cv, then the unforced differential equation is:" 100 PRINT 110 PRINT " m x'' + c x' + k x = 0" 120 PRINT "where ' is the time differentiation operator d/dt." 130 CALL Pager 140 PRINT "Next we apply an external driving force F(ext)." 150 PRINT "If that force is periodic, it can be decomposed via Fourier" 160 PRINT "techniques into a frequency spectrum of many sinusoids. Thus" 170 PRINT "our first study of damping will deal with one external driving" 180 PRINT "frequency w. Thus F(ext) = F(0)exp(iwt)." 190 PRINT "Then the differential equation is:" 200 PRINT 210 PRINT " m x'' + c x' + k x = F(0)exp(iwt)" 220 PRINT "This equation has a transient response & a steady state solution." 230 PRINT "We will deal only with the latter." 231 PRINT 240 PRINT "We assume a solution like x(t)=Aexp(i(wt-p))" 250 PRINT "where A is the amplitude ( a function of w) and p is a phase." 260 PRINT "Substitution into the DEQ, yields" 270 PRINT " -mw^2 +iwca +kA = F(0) exp(ip) = F(0)(cos(p)+i sin(p))" 290 PRINT 300 PRINT "Equating REAL and IMAGINARY parts yields two equations." 310 CALL Pager 320 PRINT "These are:" 330 PRINT " A(k-mw^2) =F(0) cos(p) and cwA = F(0) sin(p)" 350 PRINT "Solution is straight forward and yields:" 370 PRINT " cw" 380 PRINT " tan(p) = -------------" 390 PRINT " k - m w^2" 400 PRINT "and" 410 PRINT " F(0)" 420 PRINT " A(w) = -------------------" 430 PRINT " [ (k-mw^2)^2 + c^2 w^2 ] ^(1/2)" 440 PRINT "These can be rearranged. Recall w0^2 = k/m and let g=c/(2m)" 450 PRINT "then we get:" 460 PRINT " 2gw" 470 PRINT " tan(p) = -------------" 480 PRINT " w0^2-w^2" 490 PRINT "and" 500 PRINT " F(0)/m" 510 PRINT " A(w) = ---------------" 520 PRINT " D(w)" 530 PRINT "where D(w)^2 = [ (w0^2-w^2)^2 + 4 g^2 w^2]." 550 CALL Pager 600 PRINT "D(w) is a RESONANCE DEMONATOR. That is, it can nearly vanish for" 610 PRINT "some values of the applied angular frequency w. When this is the" 620 PRINT "situation A(w), the amplitude of vibration, can become large. A" 630 PRINT "RESONANCE occurs. The frequency that this occurs (if it does), is" 640 PRINT "called the RESONANCE FREQUENCY wr. " 650 PRINT "This frequency occurs where w--> wr, e.g., when A(w) is a maximum." 660 PRINT "Thus we set dA(w)/dw = 0, and solve for w (then called wr). Thus" 670 PRINT "we find " 680 PRINT " wr = SQRT( w0^2-2g^2)" 690 PRINT "For WEAK DAMPING, g= w0/sqrt(2)."; 730 PRINT "The AMPLITUDE AT THE RESONANCE PEAK is call Amax. It is found by" 740 PRINT "back substituting wr into the A(w) equation. It is given by:" 750 PRINT " F(0)/m F(0)" 760 PRINT " A max = -------------------- = --------------" 770 PRINT " 2gsqrt( w0^2-g^2) c wd" 780 PRINT 790 CALL Pager 800 PRINT "The sharpness of the resonace peak is of great interest. In" 810 PRINT "applications in mechanical and civil engineering we wish it to" 820 PRINT "be broad. In many examples in electrical engineering and atomic" 830 PRINT "physics, we wish it to be very small. Here, we consider the case" 840 PRINT "of WEAK DAMPING. [ g<< w0 ] Here, w0^2-w^2 = about 2(w0(w0-w))" 850 PRINT "and gw = about g w0. Then A max is approximately equal to:" 860 PRINT " A max g" 870 PRINT " A (w) = ---------------------------" 880 PRINT " Sqrt[ (w0-w)^2 + g^2 ]" 890 PRINT "When w=w0+/-g , then A^2= 0.5 * Amax^2. Note 0.5kA^2 is Potential" 900 PRINT "Energy in undamped motion, so A^2 is proportional to energy." 910 PRINT "Thus g(gamma)=c/2m is a frequency and the curve falls to its" 920 PRINT "HALF HEIGHT when w-> w0+/- g. Thus 2g = FULL WIDTH at HALF HEIGHT." 930 PRINT "and g = Half Width at HALF HEIGHT (HWHH)." 940 PRINT 950 CALL Pager 1010 PRINT "Another way to look at the resonance is in terms if the Q factor." 1020 PRINT "Q is called the Quality factor. It is used a lot in circuits and" 1030 PRINT "acoustical resonators. Q is defined as the ratio of 2 frequencies." 1040 PRINT "and thus is dimensionless. Q = wr/2g, the ratio of the resonant" 1050 PRINT "frequency to that of the FULL WIDTH @ HALF MAXIMUM of the resonance" 1060 PRINT "response curve, e.g. the amplitude curve. The fractional width of " 1070 PRINT "resonance peak is the ratio delta w over width & is about" 1080 PRINT "equal to 1/Q."; 1090 PRINT " Q's of 1000 are common for quartz oscillators. For specially " 1100 PRINT "constructed laser cavities they can be made to be 10^11 !" 1110 CALL Pager 1120 PRINT "THE PHASE ANGLE:" 1130 PRINT " THE PHASE DIFFERENCE between the driving function and the result-" 1140 PRINT "ing oscillation is called p above. At low values of w, p remains " 1150 PRINT "small. As w -> wr, the Tan^-1 goes to 90 degrees, and as w-> big," 1160 PRINT "w aproaches 180 degrees." 1170 CALL Pager 1200 PRINT "Now we will do a simulation of these details numerically and" 1205 DIM x(500), y(500), P(500) 1210 PRINT "write a plotting file for your use. We will still assume that " 1220 PRINT "the model is a spring system. You could input details of another" 1230 PRINT "1D system here, IF YOU KNOW THE ANALOGS." 1240 PRINT 1250 PRINT " We will need to know for the STREADY STATE SOLUTION SEVERAL" 1260 PRINT "parameters. They are in S. I. units." 1270 PRINT "Input the mass in Kg?"; 1280 INPUT M 1290 PRINT "Input the spring constant?"; 1300 INPUT K 1310 PRINT "Input the damping constant c?"; 1320 INPUT C 1330 PRINT "Input the driving force F(0) in Newtons?"; 1340 INPUT F0 1350 PRINT "We will vary the applied frequency from w0-5g to w0+5g." 1360 PRINT 1370 CALL Pager 1380 W0 = SQR(K / M) 1390 G = C / (2 * M) 1400 WR = SQR(W0 ^ 2 - 2 * G ^ 2) 1410 WD = SQR(W0 ^ 2 - G ^ 2) 1420 PRINT "For thes values of M, k and C, we find:" 1430 PRINT "UNDAMPED ANGULAR FREQUENCY W0 ="; W0 1440 PRINT " DAMPED ANGULAR FREQUENCY Wd ="; WD 1450 PRINT "RESONANT ANGULAR FREQUENCY Wr ="; WR 1460 PRINT 1470 IF G < .1 * W0 THEN PRINT " WEAK DAMPING CASE" 1480 AMAX = F0 / (C * WD) 1490 PRINT "The MAXIMUM DISPLACEMENT ="; AMAX 1500 PRINT "The g = "; G; 1520 PRINT " If weak damping it is the HWHH." 1530 Q = WD / (2 * G) 1532 PRINT " Q FACTOR ="; Q 1540 TOP = F0 / M 1545 COE = 4 * G * G 1550 W02 = W0 * W0 1560 I = 0 1600 FOR W = W0 - 5 * G TO W0 + 5 * G STEP .05 * G 1610 I = I + 1 1620 x(I) = W 1622 W2 = W * W 1670 BOT = (W02 - W2) ^ 2 + COE * W2 1680 BOT = SQR(BOT) 1690 y(I) = TOP / BOT REM need an ATN2 function REM and a covert to degrees function REM atn2(x,y) returen the proper sign and quadrant of a angle 1692 P(I) = ATN2(2 * G * W, (W02 - W2)) * 180! / pi 1700 NEXT W 1702 NPTS = I 1705 PRINT "We will now write files for use by graphics programs." 2100 REM TEKTRONIX OUTPUT HERE WRITE LABELS AND THE DATA IN PLTTEK.DAT 2110 OPEN "PLTTEK.DAT" FOR OUTPUT AS #3 2112 PRINT #3, 2 2120 PRINT #3, 1 2122 PRINT #3, NPTS 2130 PRINT #3, " AMPLITUDE" 2140 PRINT #3, "POSITION " 2150 PRINT #3, " M="; M; " C="; C; " K="; K; " WR="; WR 2155 PRINT #3, "DRIVEN DAMPED HARMONIC OSCILLATOR" 2160 FOR J = 1 TO NPTS 2170 PRINT #3, x(J), y(J) 2172 NEXT J 2174 PRINT #3, 1 2176 PRINT #3, NPTS 2181 PRINT #3, " PHASE DIFFERENCE DEGREES" 2183 PRINT #3, "ANGULAR FREQUENCY (RADS/S)" 2184 PRINT #3, " M="; M; " C="; C; " K="; K; " WR="; WR 2185 PRINT #3, "DRIVEN DAMPED HARMONIC OSCILLATOR" 2190 FOR I = 1 TO NPTS 2192 PRINT #3, x(I), P(I) 2194 NEXT I CLOSE #3 2200 PRINT "Now run WPLOT." CHAIN "TEK2WPLT.EXE" 2800 STOP FUNCTION ATN2 (x, y) pi = 4! * ATN(1!) R = SQR(x * x + y * y) SA = y / R CA = x / R REM TWO ARGUMENT ARCTAN REM CA=COS(A);SA=SIN(A);RETURNS A IN PROPER QUADRANT IF ABS(CA) < 1E-08 THEN 3781 3761 a3 = ATN(SA / CA) GOTO 3821 3781 IF SA > 0 THEN a3 = .5 * pi 3801 IF SA < 0 THEN a3 = 1.5 * pi 3821 IF CA < 0 THEN 3881 3841 IF CA > 0 AND SA < 0 THEN 3901 3861 GOTO 3921 3881 a3 = a3 + pi ATN2 = a3 EXIT FUNCTION 3901 a3 = a3 + 2 * pi 3921 ATN2 = a3 END FUNCTION SUB Pager 18600 PRINT , "Push RETURN to CONTINUE, S to STOP."; 18700 INPUT Q$ 18720 IF Q$ = "S" OR Q$ = "STOP" OR Q$ = "s" THEN END 18800 PRINT CHR$(12) END SUB