/*********************************************************
*                    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_final and number of points:                     *
* 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         *
*                                                        *
* Press any key to continue                              *
* (For a better accuracy, increase number of points).    *
*                                                        *
* ------------------------------------------------------ *
* Reference: "Problem Solving with Fortran 90 By David   *
*             R. Brooks, Springer-Verlag New York, 1997" *
*                                                        *
*                    C++ Version By J-P Moreau, Paris.   *
*                           (www.jpmoreau.fr)            *
*********************************************************/
#include <stdio.h>
#include <math.h>
 

//------------------------------------
      double AofT(double x,double v,double D,double K,
		          double M, double F)  {
// Calculate acceleration for "mass and spring" problem.
        return -(-F+D*v+K*x)/M;
      }
//---------------------------------------------
      void MassAndSpring(double *x,double *v, double D,double K,
		                 double M,double F,double dt)  {
// Calculate motion for mass and spring problem with constant force term.
        double k1_x,k2_x,k3_x,k4_x,k1_v,k2_v,k3_v,k4_v;
// Runge-Kutta coefficients...
        k1_x=*v;
        k2_x=*v+AofT(*x,*v,D,K,M,F)*dt/2.0;
        k3_x=*v+AofT(*x,*v,D,K,M,F)*dt/2.0;
        k4_x=*v+AofT(*x,*v,D,K,M,F)*dt;
        k1_v=AofT(*x,*v,D,K,M,F);
        k2_v=AofT(*x+(*v*dt/2.0),*v+k1_v*dt/2.0,D,K,M,F);
        k3_v=AofT(*x+(*v*dt/2.0),*v+k2_v*dt/2.0,D,K,M,F);
        k4_v=AofT(*x+(*v*dt),*v+k3_v*dt,D,K,M,F);
// Propagate solution...
        *x=*x+(k1_x+2.0*k2_x+2.0*k3_x+k4_x)*dt/6.0;
        *v=*v+(k1_v+2.0*k2_v+2.0*k3_v+k4_v)*dt/6.0;
      } // MassAndSpring()


void main()  {
/********************************************************
* 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                              *
********************************************************/
      double i,q,L,C,R,V,t,dt,t_final;
      double k1_i,k2_i,k3_i,k4_i,pi=3.1415926535;
      int j,n;
      char choice;

// Choose circuit type... 

      printf("\n Specify [o]scillating or [n]o oscillating term... ");
      scanf("%c",&choice);
      if (choice=='o' || choice=='O') {
// Ld^2q/dt^2+Rdq/dt+q/C=V
		do {
          printf(" Give V or F, L or M, R or 1/K, C or D (4L/C-R^2 > 0):\n");
          printf(" "); scanf("%lf %lf %lf %lf",&V,&L,&R,&C);
          printf(" one period at t=%f s\n", 4*pi*L/sqrt(4.0*L/C-R*R));
          printf(" Give t_final and number of points:\n");
          printf(" "); scanf("%lf %d", &t_final,&n);
        } while (4.0*L/C-R*R <= 0.0);
        q=0.0; i=C*V*R/2.0/L; t=0.0;   //Initial values
        dt=t_final/n;
        printf("\n       time        q or d    i or speed  analytic q or d\n");
        for (j=1; j<=n; j++) {
          MassAndSpring(&q,&i,R,1.0/C,L,V,dt);
          t=t+dt;
          printf(" %12.6f%12.4f%12.4f%12.4f\n", t, q, i, C*V*(1.0-exp(-R*t/2.0/L)
			     *cos(sqrt(4.0*L/C-R*R)*t/2.0/L)));
        }
      }
      else if (choice=='n' || choice=='N') {
// Ldi/dt+Ri=V (no oscillating term)...
        printf(" Give V or F, L or M, R or 1/K:\n");
        printf(" "); scanf("%lf %lf %lf", &V,&L,&R);
        printf(" time constant at t=%f s\n", L/R);
        printf(" Give t_final and number of points:\n");
        printf(" "); scanf("%lf %d", &t_final, &n);
        i=0.0; t=0.0;   //Initial values
        dt=t_final/n;
        printf("\n       time      i or speed  analytic i\n");
        for (j=1; j<=n; j++) {
// Runge-Kutta coefficients...
          k1_i=(V-R*i             )/L;
          k2_i=(V-R*(i+k1_i*dt/2.0))/L;
          k3_i=(V-R*(i+k2_i*dt/2.0))/L;
          k4_i=(V-R*(i+k3_i*dt)   )/L;
// Propagate solution...
          i=i+(k1_i+2.0*k2_i+2.0*k3_i+k4_i)*dt/6.0;
          t=t+dt;
          printf(" %12.6f,%12.4f%12.4f\n", t, i, V/R*(1.0-exp(-R*t/L)));
        }
      }
      else
        printf(" No such choice.  Try again...\n");
      printf("\n");
}

// end of file circuit.cpp