/* square_wave_gen.c  using Fourier series         */
/* a1=4/pi  a3=a1/3  a5=a1/5 ...                   */
/* xn = sum i=1, 3, 5, ... i<N  ai*cos(2*pi*i*n/N) */
/*  n=0..N-1                 -3, -7, -11           */

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

#define Pi 3.1415926535897932384626433832795028841971

int main(int argc, char *argv[])
{
  int i, n, N;
  double sumc, sums, a1;
  int sgn;

  printf("square_wave_gen.c \n");
  N = 8;
  a1 = 4.0/Pi;
  for(N=8; N<=64; N=2*N)
  {
    printf("%3d points, cos %d up, 0, %d down, 0, %d up, 0, %d down ... \n",
	   N, N/4, N/2-1, N/2-1, N/2-1);
    printf("%3d terms,  sin 0, %d up, 0, %d down, 0, %d up, 0, %d down ... \n",
	   N, N/2-1, N/2-1, N/2-1, N/2-1);
    for(n=0; n<=2*N; n++)
    {
      sumc = 0.0;
      sums = 0.0;
      sgn = 1;
      for(i=1; i<2*N; i=i+2) /* could be i<50 always for 3 digit */
                             /* 2*N+3 for 0.99 1.01 -0.99 -1.01  */
      {
        sumc = sumc + (a1/(double)(i*sgn))*
                    cos((2.0*Pi*(double)i*(double)n)/(double)N);
        sgn = -sgn;
        sums = sums + (a1/(double)(i))*
                    sin((2.0*Pi*(double)i*(double)n)/(double)N);
      }
      printf("n=%3d, val_cos=%9.4f   val_sin=%9.4f \n", n, sumc, sums);
    }
    printf("\n");
  }
  return 0;
}