'**********************************************************
'*                    PROGRAM CIRCUIT                     *
'* Use Runge-Kutta method to solve a LRC circuit or equi- *
'* valent damped mass-spring problem.                     *
'* ------------------------------------------------------ *
'* SAMPLE RUN:                                            *
'*                                                        *
'* (Solve mass-spring system with:                        *
'*  F=75 N, M=50 kg, K=100 N (C=0.01), D=0.05)            *
'*                                                        *
'* Specify [o]scillating or [n]o oscillating term...      *
'* o                                                      *
'* Give V or F, L or M, R or D, C or 1/K (4L/C-R^2 > 0):  *
'* 75 50 0.05 0.01                                        *
'* one period at t=   4.442879    s                       *
'* Give t* 10 10                                          *
'*                                                        *
'*     time        q or d    i or speed  analytic q or d  *
'*   1.000000      0.7504      1.4993      0.6331         *
'*   2.000000      2.2485     -0.0017      1.4628         *
'*   3.000000      0.7483     -2.9955      1.0890         *
'*   4.000000     -2.2440      0.0055      0.1436         *
'*   5.000000      0.7555      5.9850      0.2223         *
'*   6.000000      6.7320     -0.0150      1.1913         *
'*   7.000000      0.7351    -11.9580      1.4147         *
'*   8.000000    -11.2021      0.0378      0.5163         *
'*   9.000000      0.7878     23.8922      0.0131         *
'*  10.000000     24.6302     -0.0915      0.7537         *
'*                                                        *
'* (For a better accuracy, increase number of points).    *
'*                                                        *
'*                  Basic version by J-P Moreau, Paris.   *
'*                           (www.jpmoreau.fr)            *
'**********************************************************
'Program Circuit;
DEFDBL A-H, O-Z
DEFINT I-N

'main program
'*********************************************************
'* Use Runge-Kutta method to solve LRC circuit problems. *
'* Variable equivalences with mass-and-spring problem:   *
'* V   => force F                                        *
'* q   => displacement x                                 *
'* i   => velocity v                                     *
'* L   => mass m                                         *
'* R   => damping constant D                             *
'* 1/C => spring constant K                              *
'********************************************************}
     
' xi,q,XL,C,R,v,t,dt,t,tfinal: REAL
' xk1,xk2,xk3,xk4: REAL
' j,n: INTEGER;
' choice$: CHAR

' Choose circuit type...
    CLS
    pi = 3.1415926535#
    PRINT
    INPUT " Specify [o]scillating or [n]o oscillating term... ", choice$
    IF choice$ = "o" OR choice$ = "O" THEN
' Ld^2q/dt^2+Rdq/dt+q/C=V
10    PRINT " Give V or F, L or M, R or 1/K, C or D (4L/C-R^2 > 0):"
      INPUT " ", vv, XL, R, C
      F = vv: XM = XL: XK = 1# / C: D = R
      PRINT " one period at t="; 4# * pi * XL / SQR(4# * XL / C - R * R); " s"
      PRINT " Give t_final and number of points:"
      INPUT " ", tfinal, n
      IF 4# * XL / C - R * R < 0# THEN GOTO 10
      x = 0#: v = C * vv * R / 2# / XL: t = 0# 'Initial values
      dt = tfinal / n
      PRINT
      PRINT " time    q or d      i or speed   analytic q or d"
      F1$ = "  ##  ####.######  ####.######   ####.######"
      FOR j = 1 TO n
        GOSUB 1000  'MassAndSpring(x,v,R,1./C,L,vv,dt)
        t = t + dt
        PRINT USING F1$; t; x; v; C * vv * (1# - EXP(-R * t / 2# / XL) * COS(SQR(4# * XL / C - R * R) * t / 2# / XL))
      NEXT j
    ELSEIF choice$ = "n" OR choice$ = "N" THEN
' Ldi/dt+Ri=V (no oscillating term)...
      PRINT " Give V or F, L or M, R or 1/K:"
      INPUT " ", v, XL, R
      PRINT
      PRINT " time constant at t="; XL / R; " s"
      PRINT " Give t_final and number of points:"
      INPUT " ", tfinal, n
      xi = 0#: t = 0#  'Initial values
      dt = tfinal / n
      PRINT
      PRINT " time   i or speed   analytic i"
      F$ = " ##  ####.######  ####.######"
      FOR j = 1 TO n
' Runge-Kutta coefficients...
        xk1 = (v - R * xi) / XL
        xk2 = (v - R * (xi + xk1 * dt / 2#)) / XL
        xk3 = (v - R * (xi + xk2 * dt / 2#)) / XL
        xk4 = (v - R * (xi + xk3 * dt)) / XL
' Propagate solution...
        xi = xi + (xk1 + 2# * xk2 + 2# * xk3 + xk4) * dt / 6#
        t = t + dt
        PRINT USING F$; t; xi; v / R * (1# - EXP(-R * t / XL))
      NEXT j
    ELSE
      PRINT " No such choice.  Try again..."
    END IF
END

'-------------------------------------------------------
500 'FUNCTION AofT(x,v,D,K,M,F:REAL):REAL;
' Calculate acceleration for "mass and spring" problem.
  AofT = -(-F + D * v1 + XK * x1) / XM
RETURN
'-------------------------------------------------------------------------
1000 'PROCEDURE MassAndSpring(VAR x,v:REAL; D,K,M,F,dt:REAL);
' Calculate motion for mass and spring problem with constant force term.
' xk1x,xk2x,xk3x,xk4x,xk1v,xk2v,xk3v,xk4v: REAL
' Runge-Kutta coefficients...
  xk1x = v
  x1 = x: v1 = v: GOSUB 500
  xk2x = v + AofT * dt / 2#
  xk3x = v + AofT * dt / 2#
  xk4x = v + AofT * dt
  xk1v = AofT
  x1 = x + v * dt / 2#: v1 = v + xk1v * dt / 2#
  GOSUB 500: xk2v = AofT
  v1 = v + xk2v * dt / 2#
  GOSUB 500: xk3v = AofT
  x1 = x + v * dt: v1 = v + xk3v * dt
  GOSUB 500: xk4v = AofT
' Propagate solution...
  x = x + (xk1x + 2# * xk2x + 2# * xk3x + xk4x) * dt / 6#
  v = v + (xk1v + 2# * xk2v + 2# * xk3v + xk4v) * dt / 6#
RETURN

'end of file circuit.bas