/***********************************************************************
 *
 * File: rlckreis-rk4.c (24-May-2001)
 *                      (03-Jun-2002)
 *                      (27-May-2005)
 *
 * 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
 *
 *   U(t) = U_max*cos(w*t)   fuer cos(w*t) > 0
 *
 *          0                sonst
 *
 * Integration mit Runge-Kutta Verfahren 4. Ordnung
 *
 *   R....................Widerstand
 *   L....................Induktivitaet
 *   C....................Kapazitaet
 *   U_max................Spannung (Maximalamplitude)
 *   w....................Kreisfrequenz
 *   dt...................Zeitschritt
 *   t_0,t_max............Anfangs/Endzeitpunkt
 *   "rlckreis-rk4.dat"...Ausgabefile
 *
 **********************************************************************/

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

float c,r,umax,w,xl;

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

void f(float *dxidt,float *dqdt,float xi,float q,float t) {

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

  float u;
  
  if(cos(w*t)>0.0) {
    u=umax*cos(w*t);
  } else {
    u=0.0;
  }

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

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

void rk4(float *u1n,float *u2n,float tn,float u1,float u2,float t) {

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

  float dt,f11,f12,f21,f22,f31,f32,f41,f42;
  
  dt=tn-t;

  f(&f11,&f12,u1,u2,t);
  f(&f21,&f22,u1+0.5*dt*f11,u2+0.5*dt*f12,t+0.5*dt);
  f(&f31,&f32,u1+0.5*dt*f21,u2+0.5*dt*f22,t+0.5*dt);
  f(&f41,&f42,u1+dt*f31,u2+dt*f32,t+dt);

  *u1n=u1+dt*(f11+2.0*f21+2.0*f31+f41)/6.0;
  *u2n=u2+dt*(f12+2.0*f22+2.0*f32+f42)/6.0;
  
}

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

int main() {

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

  FILE *fout;
  int it;
  float dt,q,qn,t,t0,tmax,tn,xi,xin;

  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=8.0*atan(1.0)*w;
  printf("t_0,dt,t_max(2)=");
  scanf("%f,%f,%f",&t0,&dt,&tmax);

  fout=fopen("rlckreis-rk4.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;t<tmax;it=it+1) {
    
    tn=t0+it*dt;
    rk4(&xin,&qn,tn,xi,q,t);
    
    t=tn;
    xi=xin;
    q=qn;
    fprintf(fout," %10.5f %14.7e %14.7e\n",t,xi,q);

  }

  close(fout);
  
  return 0;

}

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