/*
 *                            COPYRIGHT
 *
 *  nsines - Noisy sum of sampled sine waves program.
 *  Copyright (C) 2007 Exstrom Laboratories LLC
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  A copy of the GNU General Public License is available on the internet at:
 *
 *  http://www.gnu.org/copyleft/gpl.html
 *
 *  or you can write to:
 *
 *  The Free Software Foundation, Inc.
 *  675 Mass Ave
 *  Cambridge, MA 02139, USA
 *
 *  You can contact Exstrom Laboratories LLC via Email at:
 *
 *  stefan(AT)exstrom.com
 *
 *  or you can write to:
 *
 *  Exstrom Laboratories LLC
 *  662 Nelson Park Dr.
 *  Longmont, CO 80503, USA
 *
 */

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

#define MAXLINESIZE 128

double frand( void );

int main( int argc, char *argv[] )
{
  int np; /* number of data points to generate */
  double nmax; /* maximum amplitude of noise */
  int ns; /* number of sine waves used to generate data */
  double *A;   /* array of amplitudes of sine waves used */
  double *w;   /* array of frequencies of sine waves used */
  double *phi; /* array of phases of sine waves used */
  int i;
  int n; /* point number */
  double f; /* value of current data point generated */
  double noise; /* random number [-nmax...+nmax] */
  unsigned int seed; /* seed for random number generator */
  FILE *fp;
  char linebuff[MAXLINESIZE];

  if( argc < 3 )
  {
    printf("\n");
    printf("nsines produces a noisy sum of ns sampled\n");
    printf("sine waves of the form:\n");
    printf("fn = A1*sin(n*w1+phi1) + A2*sin(n*w2+phi2) + ... + noise\n");
    printf("usage: nsines infile outfile\n");
    printf("  infile = input file with the following structure:\n");
    printf("    np: number of data points to generate\n");
    printf("    nmax: maximum amplitude of noise\n");
    printf("    seed: unsigned int [0-2147483648] for random # generator\n");
    printf("    ns: number of sine waves to use\n");
    printf("    Following this are three columns of ns rows:\n");
    printf("      Amplitude  Frequency[rad/s]  Phase[rad] \n");
    printf("      A1         w1                phi1\n");
    printf("      '2         '2                  '2\n");
    printf("      '          '                   '\n");
    printf("      'ns        'ns                 'ns\n");
    printf("      w & phi ranges are 0...1, interpreted as multiples of Pi\n");
    printf("      Amplitude range is any real number\n");
    printf("  outfile = file where output is written\n");
    return(-1);
  }

  fp = fopen(argv[1], "r");
  if( fp == 0 )
  {
    perror( "Unable to open input file" );
    return(-1);
  }

  printf( "\nReading input file: %s\n", argv[1] );

  /* skip header */
  for( i = 0; i < 10000; ++i ) /* Assuming header < 10000 lines */
  {
    fgets( linebuff, MAXLINESIZE, fp );
    if( linebuff[0] != '#' ) break;
    printf( "%s", linebuff );
  }

  /* get the first number that's already in linebuff */
  sscanf( linebuff, "%d", &np );

  fscanf( fp, "%lf", &nmax );
  fscanf( fp, "%d", &seed );
  fscanf( fp, "%d", &ns );

  A = (double *)malloc( ns * sizeof(double) );
  w = (double *)malloc( ns * sizeof(double) );
  phi=(double *)malloc( ns * sizeof(double) );

  for( i = 0; i < ns; ++i )
  {
    fscanf( fp, "%lf %lf %lf", &(A[i]), &(w[i]), &(phi[i]) );
    w[i] *= M_PI;
    phi[i] *= M_PI;
  }

  fclose( fp ); /* close input file */

  printf( "Output file = %s\n", argv[2] );

  fp = fopen(argv[2], "w");
  if( fp == 0 )
  {
    perror( "Unable to open output file" );
    return(-1);
  }
   
  printf( "Saving data to file: %s\n", argv[2] );
  fprintf( fp, "# %s  :program that created this file\n", argv[0] );
  fprintf( fp, "# %s  :input file used\n", argv[1] );
  fprintf( fp, "# %20d  :np, number of data points in this file\n", np );
  fprintf( fp, "# %20f  :nmax, maximum noise amplitude\n", nmax );
  fprintf( fp, "# %20d  :seed, random number generator seed\n", seed );
  fprintf( fp, "# %20d  :ns, number of sine waves used\n", ns );
  fprintf( fp, "# Num     Amplitude     Frequency[rad/s]     Phase[rad]\n" );
  for( i = 0; i < ns; ++i )
  {
    fprintf( fp, "# %d  %13f %13f        %13f\n", i+1, A[i], w[i], phi[i] );
  }
  
  srandom(seed);

  for( n = 0; n < np; ++n ) /* for each point to be generated */
  {
    f = 0.0;
    for( i = 0; i < ns; ++i ) /* add individual sine waves together */
    {
      f += A[i]*sin(n*w[i]+phi[i]);
    }
    noise = nmax * frand();
    f += noise;
    fprintf( fp, "%lf\n", f );
  }
  fclose( fp );
  free( A );
  free( w );
  free( phi );
  return( 0 );
} /* end of program */
/*********************************************************************/
double frand( void )
{
/* frand produces a random number in the range of -1.0 to +1.0 */

double f; /* final result */

f = 2.0 * (double)random() / (double)RAND_MAX - 1.0;
return ( f );

} /* end function frand */