/***********************************************************************
 *
 * File: rlkreis-rk2.c (11-May-2001)
 *                     (17-May-2001)
 *                     (06-May-2005)
 *                     (27-May-2005)
 *
 * Einschwingverhalten eines RL-Kreises
 *
 *   L*dI/dt + R*I = U,     I(t_0)=0
 *
 * wo
 *
 *   U(t) = U_max*cos(w*t)   fuer cos(w*t) > 0
 *
 *          0              sonst
 *
 * Integration mit Verfahren 2. Ordnung mit Hilfsschritt
 * (Runge-Kutta Verfahren 2. Ordnung)
 *
 * R...................Widerstand
 * L...................Induktivitaet
 * U_max...............Spannung (Amplitude)
 * w...................Kreisfrequenz
 * dt..................Zeitschritt
 * t_0.................Startzeitpunkt
 * t_max...............Endzeitpunkt
 * "rlkreis-rk2.dat"...Ausgabefile
 *
 **********************************************************************/

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

float r,w,umax,xl;

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

float f(float xi,float t) {

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

  if(cos(w*t)>=0.0) {
    return ((-r*xi+umax*cos(w*t))/xl);
  } else {
    return ((-r*xi)/xl);
  }

}

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

int main() {

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

  FILE *fout;
  int it;
  float dt,t,t0,t1,tmax,tn,xi,xi1,xin;

  printf("         R(Ohm)=");
  scanf("%f",&r);
  printf("          L(Hy)=");
  scanf("%f",&xl);
  printf("       U_max(V)=");
  scanf("%f",&umax);
  printf("     w/2pi(1/s)=");
  scanf("%f",&w);
  w=8.0*atan(1.0)*w;
  printf("t_0,dt,t_max(s)=");
  scanf("%f,%f,%f",&t0,&dt,&tmax);

  fout=fopen("rlkreis-rk2.dat","w");
  fprintf(fout,"#\n");
  fprintf(fout,"#   t(s)        I(A)\n");
  fprintf(fout,"#\n");

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

  for(it=1;t<=tmax;it=it+1) {
    t1=t+0.5*dt;
    xi1=xi+0.5*dt*f(xi,t);          /* Halbschritt */
    tn=t0+it*dt;
    xin=xi+dt*f(xi1,t1);            /* Vollschritt */
    t=tn;
    xi=xin;
    fprintf(fout,"%10.5f %12.5e\n",t,xi);
  }

  fclose(fout);

  return 0;

}

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