REM
REM     GE Copeland
REM     March 23, 1999
REM     Old Dominion University
REM     Dept of Physics
REM Updated 3/26/2001

PRINT CHR$(12);
pi = 4 * ATN(1!)
PRINT "Series RLC Underdamped Circuit"
PRINT "Input C in Farads amd L in Henries."
INPUT C, L
 
   Rmax = SQR(4! * L / C)
1  PRINT " Max value of R for underdamped oscillation = "; Rmax; " Ohms"
   PRINT " Input R in Ohms";
   INPUT R
IF R >= Rmax THEN GOTO 1
   wp = SQR(1! / (L * C) - R ^ 2 / (4! * L ^ 2))
   wfree = SQR(1! / (L * C))
   PRINT "Angular frequency for underdamped oscillation = "; wp; "radians per second."
   PRINT "in ordinary frequency = "; wp / (2 * pi); " Hz"
   PRINT "Free Oscillation angular frequency ="; wfree; "rad/s = "; wfree / (2 * pi); " Hz"
   tau = R / (2! / L)
   tau = 1! / tau

   PRINT "Eponential decay constant ="; tau; "seconds."
   PRINT
   PRINT "INITIAL CONDITIONS"
   PRINT "Input the initial charge on the capacitor in coulombs";
   INPUT q0
  
   PRINT "Input the initial current in the circuit in Amperes";
   INPUT I0

   A = SQR(q0 ^ 2 + (1! / wp ^ 2) * (I0 + q0 / tau) ^ 2)
   phi = ATN(-(I0 + q0 / tau) / (wp * q0))
   PRINT "We find the Amplitude A ="; A; "Coulombs"
   PRINT "and the phase angle = "; phi; " radians ="; phi * 180 / pi; " degrees"
   PRINT "_______________"
   PRINT " Now we generate the data file for plotting."
   DIM t(1000), q(1000), curr(1000)
   n = 0
   FOR time = 0! TO 3! * tau STEP tau / 250!
   n = n + 1
 t(n) = time
 curr(n) = -(A / tau) * EXP(-time / tau) * COS(wp * time + phi) - wp * A * EXP(-time / tau) * SIN(wp * time + phi)
 q(n) = A * EXP(-time / tau) * COS(wp * time + phi)
   NEXT time

2  PRINT "Select either the 1. Charge or the 2. Current to plot vs time. Input 1 or 2"
  INPUT p$
  IF p$ = "1" OR p$ = "2" THEN
  GOTO 3
  ELSE
  GOTO 2
  END IF

3  IF LEFT$(p$, 1) = "1" THEN
  ip = 1
  ELSE
  ip = 2
  END IF

     OPEN "plttek.dat" FOR OUTPUT AS #2
     PRINT #2, 1
     PRINT #2, 1
     PRINT #2, n
     IF ip = 1 THEN
     PRINT #2, "Charge <Coulombs>"
     ELSE
     PRINT #2, "Current <amps>"
     END IF

     PRINT #2, "Time in seconds"
     PRINT #2, "R = "; R; " Ohms L = "; L; "H & C = "; C; " F "
     PRINT #2, "Wp ="; wp; " radians/s"
 FOR i = 1 TO n
     IF ip = 1 THEN
     PRINT #2, t(i), q(i)
     ELSE
     PRINT #2, t(i), curr(i)
     END IF
 NEXT i
   CLOSE #2
   CHAIN "Tek2WPLT.bas"
       END