/***********************************************************************
 *
 * File: rlkreis.c (11-May-2001)
 *                 (05-May-2006)
 *                 (17-May-2007)
 *                 (27-Mar-2008)
 *
 * Einschwingverhalten eines RL-Kreises
 *
 *   L*dI/dt + R*I = U,     I(t_0)=0
 *
 * bei Anregung mit
 *
 *          U_max*cos(w*t)   fuer cos(w*t) > 0
 *   U(t) =
 *          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.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() {

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

  const float pi=M_PI;
  
  FILE *fout;
  int it,itmax;
  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/(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("rlkreis.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;it<=itmax;it++) {

    t1=t+0.5*dt;            /* Hilfsschritt */
    xi1=xi+0.5*dt*f(xi,t);

    tn=t0+it*dt;            /* Vollschritt */
    xin=xi+dt*f(xi1,t1);            

    t=tn;
    xi=xin;
    fprintf(fout,"%10.5f %14.7e\n",t,xi);
  }

  fclose(fout);

  return 0;

}

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