/***********************************************************************
 *
 * File: rlckreis.c (24-May-2001)
 *                  (27-May-2005)
 *                  (27-Mar-2008)
 *
 * Einschwingverhalten eines RLC-Kreises
 *
 *   dI/dt = (U - R*I - Q/C)/L
 *
 *   dQ/dt = I
 *
 * mit Anfangsbedingung I(t_0)=0, Q(t_0)=0 und Spannungsquelle
 *
 *          U_max*cos(w*t)   fuer cos(w*t) > 0
 *   U(t) =
 *          0                sonst
 *
 * Integration mit Runge-Kutta Verfahren 4. Ordnung
 *
 *   R................Widerstand
 *   L................Induktivitaet
 *   C................Kapazitaet
 *   U_max............Spannung (Amplitude)
 *   w................Kreisfrequenz
 *   dt...............Zeitschritt
 *   t_0,t_max........Anfangs/Endzeitpunkt
 *
 *   "rlckreis.dat"...Ausgabefile
 *
 **********************************************************************/

#include <stdio.h>
#include <math.h>

#define N 2

float c,r,umax,w,xl;

/**********************************************************************/

void f(float dudt[],float u[],float t) {

/**********************************************************************/

  float dqdt,dxidt,q,v,xi;

  xi=u[0];
  q=u[1];
  
  if(cos(w*t)>0.0) {
    v=umax*cos(w*t);
  } else {
    v=0.0;
  }

  dxidt=(v-r*xi-q/c)/xl;
  dqdt=xi;

  dudt[0]=dxidt;
  dudt[1]=dqdt;

  return;
  
}

/**********************************************************************/

void rk4(float un[],float u[],float t,float dt) {

/**********************************************************************/

  int i;
  float f1[N],f2[N],f3[N],f4[N];

  f(f1,u,t);

  for(i=0;i<=N-1;i++) {
    un[i]=u[i]+0.5*dt*f1[i];
  }
  f(f2,un,t+0.5*dt);

  for(i=0;i<=N-1;i++) {
    un[i]=u[i]+0.5*dt*f2[i];
  }
  f(f3,un,t+0.5*dt);

  for(i=0;i<=N-1;i++) {
    un[i]=u[i]+dt*f3[i];
  }
  f(f4,un,t+dt);

  for(i=0;i<=N-1;i++) {
    un[i]=u[i]+dt*(f1[i]+2.0*(f2[i]+f3[i])+f4[i])/6.0;
  }

  return;
  
}

/**********************************************************************/

int main() {

/**********************************************************************/

  const float pi=M_PI;

  FILE *fout;
  int it,itmax;
  float dt,q,t,t0,tmax,xi;
  float u[2],un[2];

  printf("          R[Ohm]=");
  scanf("%f",&r);
  printf("           L[Hy]=");
  scanf("%f",&xl);
  printf("            C[F]=");
  scanf("%f",&c);
  printf("        U_max[V]=");
  scanf("%f",&umax);
  printf("   w/(2*pi)[1/s]=");
  scanf("%f",&w);
  w=2.0*pi*w;
  printf(" t_0,dt,t_max[s]=");
  scanf("%f,%f,%f",&t0,&dt,&tmax);

  itmax=(int)((tmax-t0)/dt+0.5);

  fout=fopen("rlckreis.dat","w");

  fprintf(fout,"#\n");
  fprintf(fout,"#    t[s]        I[A]           Q[C]\n");
  fprintf(fout,"#\n");
  
  t=t0;
  xi=0.0;
  q=0.0;

  fprintf(fout," %10.5f %14.7e %14.7e\n",t,xi,q);
  
  for(it=1;it<=itmax;it++) {

    u[0]=xi;
    u[1]=q;

    rk4(un,u,t,dt);
    
    t=t0+it*dt;
    xi=un[0];
    q=un[1];

    fprintf(fout," %10.5f %14.7e %14.7e\n",t,xi,q);

  }

  close(fout);
  
  return 0;

}

/**********************************************************************/