/**********************************************************************
 *
 * File: frlkreis.c (28-May-2008)
 *
 * Strom in einem RL-Kreis
 *
 *   L*dI/dt + R*I = U
 *
 * bei periodischer Anregung mit positiven Cosinus-Halbwellen
 * (Periode 1(s))
 *
 * Loesung mit Fourier-Ansatz ueber dem Intervall [-0.5(s),0.5(s)]
 *
 * nt...............Anzahl Punkte auf t-Achse
 * nw...............Anzahl Terme in Fourier-Entwicklung (w=omega)
 * R................Widerstand
 * L................Induktivitaet
 * U_max............Spannung (Amplitude)
 *
 * "frlkreis.dat"...Ausgabefile
 *
 **********************************************************************/

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

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

int main() {

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

  FILE *fp;
  int it,iw,nt,nw;
  double r,su,t,umax,xl;
  double *a;
  double complex si,z;
  double complex *c;

  printf("       nt=");
  scanf("%d",&nt);
  printf("       nw=");
  scanf("%d",&nw);
  printf("   R[Ohm]=");
  scanf("%lf",&r);
  printf("    L[Hy]=");
  scanf("%lf",&xl);
  printf(" U_max[V]=");
  scanf("%lf",&umax);

  a=(double*)malloc((nw+1)*sizeof(double));
  c=(double complex*)malloc((nw+1)*sizeof(double complex));

  a[0]=(2.0/M_PI)*umax;
  a[1]=0.5*umax;
  for(iw=2;iw<=nw;iw++) {
    if(iw%4==2) {
      a[iw]=(2.0/(M_PI*(iw*iw-1)))*umax;
    } else if(iw%4==0) {
      a[iw]=-(2.0/(M_PI*(iw*iw-1)))*umax;
    } else {
      a[iw]=0.0;
    }
  }

  for(iw=0;iw<=nw;iw++) {
    z=r+(2.0*M_PI*iw)*I;
    c[iw]=0.5*a[iw]/z;
  }

  fp=fopen("frlkreis.dat","w");

  fprintf(fp,"#\n");
  fprintf(fp,"#   t(s)       U[V]           I[A]\n");
  fprintf(fp,"#\n");

  for(it=0;it<=nt;it++) {

    t=-0.5+((double)it)/nt;

    su=0.5*a[0];
    for(iw=1;iw<=nw;iw++) {
      su=su+a[iw]*cos(2.0*M_PI*iw*t);
    }

    si=c[0];
    for(iw=1;iw<=nw;iw++) {
      si=si+(c[iw]*cexp(2.0*M_PI*iw*t*I)+\
	     conj(c[iw])*cexp(-2.0*M_PI*iw*t*I));
    }

    fprintf(fp," %8.5f %14.7e %14.7e\n",t,su,creal(si));

  }

  fclose(fp);

  return 0;

}

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