/***********************************************************************
 *
 * File: fsaite.c (04-May-2001)
 *                (06-May-2005)
 *                (27-May-2005)
 *
 * Loesung der Wellengleichung fuer die gezupfte Saite
 *
 *   (d/dt)**2 A(x,t) = c**2 * (d/dx)**2 A(x,t)
 *
 * fuer das Intervall [0,1] und c=1 mit Fourier-Reihe
 *
 * nx.............Anzahl Punkte auf x-Achse
 * nk.............Anzahl Terme in Fourier-Entwicklung
 * "fsaite.dat"...Ausgabefile
 *
 **********************************************************************/

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

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

int main() {
  
/**********************************************************************/

  int ik,ix,nk,nx;
  float pi,s,t,x;
  float *b;
  FILE *fout;
  
  printf("nx=");
  scanf("%d",&nx);
  printf("nk=");
  scanf("%d",&nk);
  printf(" t=");
  scanf("%f",&t);
  
  pi=4.0*atan(1.0);
  
  b=(float*)malloc((nk+1)*sizeof(float));
  
  for(ik=1;ik<=nk;ik=ik+2) {
    if(ik%4==1) {
      b[ik]=8.0/((ik*pi)*(ik*pi));
    } else {
      b[ik]=-8.0/((ik*pi)*(ik*pi));
    }
  }
  
  fout=fopen("fsaite.dat","w");
  
  for(ix=0;ix<=nx;ix=ix+1) {
    
    x=((float)ix)/nx;
    s=0.0;
    for(ik=1;ik<=nk;ik=ik+2) {
      s=s+b[ik]*cos(ik*pi*t)*sin(ik*pi*x);
    }
    
    fprintf(fout," %8.5f %14.7e\n",x,s);
    
  }
  
  fclose(fout);
  
  return 0;
  
}

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