/***********************************************************************
 *
 * File: heavi.c (04-May-2001)
 *               (06-May-2005)
 *               (27-May-2005)
 *
 * Fouriersynthese der Heaviside-Funktion
 *
 *   H(x) = 0 fuer x < 0
 *          1 fuer x > 0
 *
 * auf dem Intervall [-0.5,0.5]
 *
 * nx............Anzahl Punkte auf x-Achse
 * nk............Anzahl Terme in Fourier-Entwicklung
 * "heavi.dat"...Ausgabefile
 *
 **********************************************************************/

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

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

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

  const float xmax=0.5,xmin=-xmax;

  int ik,ix,nk,nx;
  float a0,dk,dx,pi,s,x;
  float *b;
  FILE *fout;
  
  printf("nx=");
  scanf("%d",&nx);
  printf("nk=");
  scanf("%d",&nk);
  
  dx=(xmax-xmin)/nx;
  pi=4.0*atan(1.0);
  dk=2.0*pi/(xmax-xmin);
  
  b=(float*)malloc((nk+1)*sizeof(float));

  a0=1.0;
  for(ik=1;ik<=nk;ik=ik+2) {
    b[ik]=(4.0/(ik*dk))/(xmax-xmin);
  }
  
  fout=fopen("heavi.dat","w");
  
  for(ix=0;ix<=nx;ix=ix+1) {
    
    x=xmin+ix*dx;
    s=0.5*a0;
    for(ik=1;ik<=nk;ik=ik+2) {
      s=s+b[ik]*sin(ik*dk*x);
    }
  
    fprintf(fout," %8.5f %14.7e\n",x,s);
    
  }
  
  fclose(fout);
  
  return 0;
  
}

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