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

#define MAX_FILENAME_LENGTH 1024

char infile_name[MAX_FILENAME_LENGTH];
char noisefile_name[MAX_FILENAME_LENGTH];
char parfile_name[MAX_FILENAME_LENGTH];
char outfile_name[MAX_FILENAME_LENGTH];
FILE *infile;
FILE *parfile;
FILE *noisefile;
FILE *outfile;
FILE *logfile;
int n_proc;
int gauss_quality;
int exit_code;                  /* exit code for error identification */
int infile_length;
int parfile_length;
int noisefile_length;
int outfile_length;
int want_log;
double amp_grid;
double freq_grid;
double phase_grid;

double *jd;                   /* time values array */
double *amp;                  /* amplitude values array */
double *freq;                 /* frequency values array */
double *phase;                /* phase values array */
double *mag;                  /* magnitude values array */
double *rms;                  /* noise levels array */
double **noisy;               /* light curve plus noise */
double **newamp;
double **newfreq;
double **newphase;
double **var_amp;
double **var_freq;
double **var_phase;

int count_file(char *fname);
void read_parfile();
void read_infile();
void read_noisefile();
int check_outfile();
void least_squares();
double sum_of_squares(int noiseindex);
void update_outfile();

/*****************************************************************************/
/******************* Welcome to the start of the program! ********************/
/*****************************************************************************/

main(int argc, char *argv[])
{

  int i;                        /* index */
  int j;                        /* index */
  int k;                        /* index */
  double nval;                  /* value from random number generator */

/* Variable initialization */

  n_proc = -1;
  infile_length = -1;           /* no input file specified yet */
  parfile_length = -1;          /* no parameter file specified yet */
  noisefile_length = -1;        /* no noise file specified yet */
  outfile_length = -1;          /* no output file specified yet */
  gauss_quality = 10;           /* default iterations */
  amp_grid = .01;               /* default amplitude grid spacing */
  freq_grid = .0001;            /* default frequency grid spacing */
  phase_grid = .001;            /* default phase grid spacing */
  exit_code = 0;                /* no errors yet */
  want_log = 0;

  printf("EPSIM - Error Propagation SIMulator\n");
  printf("Version 1.0 beta\n");
  printf("**************************************");
  printf("****************************************\n");
  printf("by Peter Reegen\n");
  printf("Department of Astronomy\n");
  printf("University of Vienna\n");
  printf("Tuerkenschanzstrasse 17\n");
  printf("1180 Vienna, Austria\n");
  printf("Release date: April 1, 2003 (NO JOKE!)\n\n");
  printf("**************************************");
  printf("****************************************\n\n");

/* Scan command line arguments */

  while (--argc > 0)
    {
      if (strcmp(*++argv, "-I") == 0)           /* scan input file name */
	{
	  --argc;
	  sprintf(infile_name, "%s", *++argv);
	  printf("Input file:          %s\n", infile_name);
	  ++infile_length;
	}
      else if (strcmp(*argv, "-L") == 0)        /* create log file */
	  ++want_log;
      else if (strcmp(*argv, "-N") == 0)        /* scan noise file name */
	{
	  --argc;
	  sprintf(noisefile_name, "%s", *++argv);
	  printf("Noise file:          %s\n", noisefile_name);
	  ++noisefile_length;
	}
       else if (strcmp(*argv, "-O") == 0)       /* scan output file name */
	{
	  --argc;
	  sprintf(outfile_name, "%s", *++argv);
	  printf("Output file:         %s\n", outfile_name);
	  ++outfile_length;
	}
     else if (strcmp(*argv, "-P") == 0)         /* scan parameter file name */
	{
	  --argc;
	  sprintf(parfile_name, "%s", *++argv);
	  printf("Parameter file:      %s\n", parfile_name);
	  ++parfile_length;
	}
      else if (strcmp(*argv, "-#") == 0)        /* scan number of processes */
	{
	  --argc;
	  n_proc = atoi(*++argv);
	  printf("Number of processes: %lu\n", n_proc);
	}
       else if (strcmp(*argv, "-a") == 0)       /* scan amplitude grid */
	{
	  --argc;
	  amp_grid = strtod(*++argv, NULL);
	  printf("Amplitude grid:      %lg\n", amp_grid);
	}
       else if (strcmp(*argv, "-f") == 0)       /* scan frequency grid */
	{
	  --argc;
	  freq_grid = strtod(*++argv, NULL);
	  printf("Frequency grid:      %lg\n", freq_grid);
	}
       else if (strcmp(*argv, "-g") == 0)       /* scan Gaussian quality */
	{
	  --argc;
	  gauss_quality = atoi(*++argv);
	  printf("Gaussian quality:    %lu\n", gauss_quality);
	}
       else if (strcmp(*argv, "-p") == 0)       /* scan phase grid */
	{
	  --argc;
	  phase_grid = strtod(*++argv, NULL);
	  printf("Phase grid           %lg\n", phase_grid);
	}
     else
	{
	  printf("\nAbnormal program termination.\n");
	  printf("Exit code 1 - ");
	  printf ("Unknown command line argument \'%s\'.\n", *argv);
	  return 1;
	}
    }

  if (want_log) logfile = fopen("epsim.log", "w");

/* Ask for parameters missing in the command line */

  if (infile_length == -1) 
    {
      printf("Please enter (path and) name of input file: ");
      scanf("%s", infile_name);
    }
  if (noisefile_length == -1) 
    {
      printf("Please enter (path and) name of noise file: ");
      scanf("%s", noisefile_name);
    }
  if (outfile_length == -1) 
    {
      printf("Please enter (path and) name of parameter file: ");
      scanf("%s", parfile_name);
    }
  if (parfile_length == -1) 
    {
      printf("Please enter (path and) name of output file: ");
      scanf("%s", outfile_name);
    }
  while (n_proc < 0) 
    {
      printf("Please enter number of processes: ");
      scanf("%lu", &n_proc);
      if (n_proc < 0) printf("\nNumber of processes < 0. Please retry.\n\n");
    }

  if (n_proc == 0) n_proc = INT_MAX;
  printf("%s: counting number of entries...", infile_name);
  if ((infile_length = count_file(infile_name)) == -1)
    {
      printf("\nAbnormal program termination.\n");
      printf("Exit code 2 - Failed to count input file entries.\n");
      return 2;
    }
  printf("done.\n%s: %lu entries found.\n", infile_name, infile_length);

  printf("%s: counting number of entries...", noisefile_name);
  if ((noisefile_length = count_file(noisefile_name)) == -1)
    {
      printf("\nAbnormal program termination.\n");
      printf("Exit code 3 - Failed to count noise file entries.\n");
      return 3;
    }
  printf("done.\n%s: %lu entries found.\n", noisefile_name, noisefile_length);

  printf("%s: counting number of entries...", parfile_name);
  if ((parfile_length = count_file(parfile_name)) == -1)
    {
      printf("\nAbnormal program termination.\n");
      printf("Exit code 4 - Failed to count parameter file entries.");
      return 4;
    }
  else if (parfile_length / 3. > parfile_length / 3 + .1)
    {
      printf("\nAbnormal program termination.\n");
      printf("Exit code 5 - Format error in parameter file %s.\n", 
	     parfile_name);
      return 5;
    }
  else parfile_length /= 3;
  printf("done.\n%s: %lu entries found.\n", parfile_name, parfile_length);

  printf("Reading parameter file...");  
  amp = (double*)calloc(parfile_length, sizeof(double));
  freq = (double*)calloc(parfile_length, sizeof(double));
  phase = (double*)calloc(parfile_length, sizeof(double));
  newamp = (double**)calloc(noisefile_length, sizeof(double*));
  newfreq = (double**)calloc(noisefile_length, sizeof(double*));
  newphase = (double**)calloc(noisefile_length, sizeof(double*));
  noisy = (double**)calloc(noisefile_length, sizeof(double*));
  var_amp = (double**)calloc(noisefile_length, sizeof(double*));
  var_freq = (double**)calloc(noisefile_length, sizeof(double*));
  var_phase = (double**)calloc(noisefile_length, sizeof(double*));
  for (i = 0; i < noisefile_length; i++) 
    {
      noisy[i] = (double*)calloc(infile_length, sizeof(double));
      newamp[i] = (double*)calloc(parfile_length, sizeof(double));
      newfreq[i] = (double*)calloc(parfile_length, sizeof(double));
      newphase[i] = (double*)calloc(parfile_length, sizeof(double));
      var_amp[i] = (double*)calloc(parfile_length, sizeof(double));
      var_freq[i] = (double*)calloc(parfile_length, sizeof(double));
      var_phase[i] = (double*)calloc(parfile_length, sizeof(double));
      for (j = 0; j < parfile_length; j++)
	{
	  var_amp[i][j] = 0;
	  var_freq[i][j] = 0;
	  var_phase[i][j] = 0;
	}
    }

  read_parfile();
  printf("done.\n");
  
  printf("Reading input file and computing signal...");  
  jd = (double*)calloc(infile_length, sizeof(double));
  mag = (double*)calloc(infile_length, sizeof(double));
  read_infile();
  printf("done.\n");

  printf("Reading noise file...");  
  rms = (double*)calloc(noisefile_length, sizeof(double));
  read_noisefile();
  printf("done.\n");

  printf("Initializing random number generator using system time...");
  int fseed;
  fseed = time(NULL);
  printf("done.\nInitial value: %lu\n", time(NULL));

  if ((outfile_length = check_outfile()) == -6)
    {
      printf("\nAbnormal program termination.\n");
      printf("Exit code 6 - Failed to create output file %s.\n", 
	     outfile_name);
      return 6;
    }
  if (outfile_length == -1) return 0;

  outfile_length = parfile_length;

  printf("Entering main procedure.\n");
  printf("Processes finished: %16lu\n", 0);

  for (i = 0; i < n_proc; i++)
    {
      if (want_log) 
	{
	  fprintf(logfile, "Entering procedure %lu of %lu.\n", 
		  i + 1, n_proc);
	  fprintf(logfile, " Generating noise data.\n");
	}
      for (j = 0; j < infile_length; j++)
	{

/*****************************************************************************/
/************* Random number generator for Gaussian distribution *************/
/*****************************************************************************/

/*

This random number generator is supposed to produce normally distributed random
numbers. It uses the central limit theorem. The number of iterations is the 
gauss_quality parameter. The result is a normal distribution with expected 
value 0 and standard deviation 1;

*/

	  nval = 0;
	  for (k = 0; k < gauss_quality; k++) 
	    nval += (double)rand() / RAND_MAX;
	  nval = (nval / gauss_quality - .5) * sqrt(12. * gauss_quality);

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

/*

Now the light curve stored in the mag array has to be distorted by the noise 
evaluated above.

*/

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

	  for (k = 0; k < noisefile_length; k++)
	      noisy[k][j] = mag[j] + rms[k] *  nval;
	}

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

/*

The next step is to compute the least squares fit. The sum of squares of the 
residuals is evaluated. Then the minimum is searched for by means of a 
step-by-step approximation.

*/

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

      if (want_log) fprintf(logfile, " Entering least squares algorithm.\n");
      least_squares();
      if (want_log) fprintf(logfile, " Writing output file.\n");
      outfile = fopen(outfile_name, "w");
      for (j = 0; j < noisefile_length; j++) 
	{
	  if (i > 1)
	    {
	      fprintf(outfile, "NOISE LEVEL: %lg     ", rms[j]);
	      fprintf(outfile, "frequency grid: %lg / amp     ", 
		      freq_grid * rms[j] / (jd[infile_length - 1] - jd[0]));
	      fprintf(outfile, "amplitude grid: %lg      ", amp_grid * rms[j]);
	      fprintf(outfile, "phase grid: %lg / amp\n", phase_grid * rms[j]);
	      fprintf(outfile, "frequency                ");
	      fprintf(outfile, "rms frequency            ");
	      fprintf(outfile, "amplitude                ");
	      fprintf(outfile, "rms amplitude            ");
	      fprintf(outfile, "phase                    ");
	      fprintf(outfile, "rms phase                \n");
	      fprintf(outfile, "                         ");
	      fprintf(outfile, "                         ");
	      fprintf(outfile, "                         ");
	      fprintf(outfile, "                         ");
	      fprintf(outfile, "                         ");
	      fprintf(outfile, "                         \n");
	    }
	  for (k = 0; k < outfile_length; k++)
	    {
	      var_amp[j][k] += pow(newamp[j][k] - amp[k], 2);
	      var_freq[j][k] += pow(newfreq[j][k] - freq[k], 2);
	      var_phase[j][k] += pow(newphase[j][k] - phase[k], 2);
	      if (i > 1) 
		{
		  fprintf(outfile, "%24lg ", freq[k]);
		  fprintf(outfile, "%24lg ", sqrt(var_freq[j][k] / i));
		  fprintf(outfile, "%24lg ", amp[k]);
		  fprintf(outfile, "%24lg ", sqrt(var_amp[j][k]/ i));
		  fprintf(outfile, "%24lg ", phase[k]);
		  fprintf(outfile, "%24lg\n", sqrt(var_phase[j][k] / i));
		}
	    }
	}
      if (want_log)
	{
	  fprintf(logfile, "\nUpdated sums  of squares:\n");
	  for (j = 0; j < noisefile_length; j++)
	    {
	      fprintf(logfile, "Noise level: %lg\n", rms[j]); 
	      for (k = 0; k < outfile_length; k++)
		{
		  fprintf(logfile, "%24lg ", var_freq[j][k]);
		  fprintf(logfile, "%24lg ", var_amp[j][k]);
		  fprintf(logfile, "%24lg\n", var_phase[j][k]);
		}
	    }
	}
      fprintf(outfile, "\n%lu processes", i + 1);
      fclose(outfile);
      printf("Processes finished: %16lu\n", i + 1);
    }

  free(noisy);  
  free(rms);
  free(amp);
  free(freq);
  free(phase);
  free(jd);
  free(mag);
  free(newamp);
  free(newfreq);
  free(newphase);
  free(var_amp);
  free(var_freq);
  free(var_phase);

  printf("\nFinished.\n\n");
  printf("**************************************");
  printf("****************************************\n\n");
  printf("Thank you for using EPSIM.\n");
  printf("Questions or complaints?\n");
  printf("Please contact Peter Reegen (reegen@astro.univie.ac.at)\n");
  printf("Bye!\n\n");

  if (want_log) fclose(logfile);

  return exit_code;

}

/*****************************************************************************/
/************************** The middle of the film ***************************/
/*****************************************************************************/

/*****************************************************************************/
/*************************** Preparing output file ***************************/
/*****************************************************************************/

int check_outfile()
{

  int i;                        /* index */
  int j;                        /* index */


  printf("Checking for output file...done.\n");
  if ((outfile = fopen(outfile_name, "r")) != NULL)
    {
      printf("Output file %s already exists!\n", outfile_name);
      printf("Please re-run with a different file name.\n");
      return -1;
    }
  if ((outfile = fopen(outfile_name, "w")) == NULL) return -6;
  for (j = 0; j < noisefile_length; j++)
    {
      fprintf(outfile, "NOISE LEVEL: %lg     ", rms[j]);
      fprintf(outfile, "frequency grid: %lg / amp    ", 
	      freq_grid * rms[j] / (jd[infile_length - 1] - jd[0]));
      fprintf(outfile, "amplitude grid: %lg     ", amp_grid * rms[j]);
      fprintf(outfile, "phase grid: %lg / amp\n", phase_grid * rms[j]);
      fprintf(outfile, "amplitude                ");
      fprintf(outfile, "rms amplitude            ");
      fprintf(outfile, "frequency                ");
      fprintf(outfile, "rms frequency            ");
      fprintf(outfile, "phase                    ");
      fprintf(outfile, "rms phase                \n");
      fprintf(outfile, "                         ");
      fprintf(outfile, "                         ");
      fprintf(outfile, "                         ");
      fprintf(outfile, "                         ");
      fprintf(outfile, "                         ");
      fprintf(outfile, "                         \n");
      for (i = 0; i < parfile_length; i++) 
	fprintf(outfile, 
		"%24lg %24lg %24lg %24lg %24lg %24lg \n",
		freq[i], 0., amp[i], 0., phase[i], 0.);
    }
  fprintf(outfile, "\n0 processes\n");
  fclose(outfile);

  return 0;
}

/*****************************************************************************/
/******************** Counting the number of file entries ********************/
/*****************************************************************************/

int count_file(char *fname)
{

  FILE *whichfile;
  int file_length;
  double file_entry;

  file_length = 0;
  while ((whichfile = fopen(fname, "r")) == NULL)
    {
      printf("\nUnable to open file %s. Please retry.\n\n", fname);
      printf("Please enter (path and) name: ");
      scanf("%s", fname);
    }

  for (file_length = 0; 
       fscanf(whichfile, "%lf", &file_entry) != EOF; ++file_length);

  fclose (whichfile);
  return file_length;
}

/*****************************************************************************/
/****************** Scan amplitudes, frequencies and phases ******************/
/*****************************************************************************/

void read_parfile()
{

  int i;                        /* index */
  int j;                        /* index */

  parfile = fopen(parfile_name, "r");
  for (i = 0; i < parfile_length; i++)
    {
      fscanf(parfile, "%lf", &freq[i]);
      fscanf(parfile, "%lf", &amp[i]);
      fscanf(parfile, "%lf", &phase[i]);
      for (j = 0; j < noisefile_length; j++)
	{
	  newamp[j][i] = amp[i];
	  newfreq[j][i] = freq[i];
	  newphase[j][i] = phase[i];
	}
    }
  fclose(parfile);

  return;
}

/*****************************************************************************/
/********************* Scan time values from input file **********************/
/*****************************************************************************/

void read_infile()
{

  int i;                        /* index */
  int j;                        /* index */

  infile = fopen(infile_name, "r");
  printf("\n");
  for (i = 0; i < infile_length; i++)
    {
      fscanf(infile, "%lf", &jd[i]);
      mag[i] = 0;
      for (j = 0; j < parfile_length; j++) 
	mag[i] += amp[j] * sin(2 * M_PI * (freq[j] * jd[i] + phase[j]));
    }
  fclose(infile);

  return;
}

/*****************************************************************************/
/********************* Scan noise levels from noise file *********************/
/*****************************************************************************/

void read_noisefile()
{

  int i;                        /* index */

  noisefile = fopen(noisefile_name, "r");
  for (i = 0; i < noisefile_length; i++) fscanf(noisefile, "%lf", &rms[i]);
  fclose(noisefile);

  return;
}

/*****************************************************************************/
/************************** Least squares algorithm **************************/
/*****************************************************************************/

void least_squares()
{

  int whichpar;

/*

This value is used as an indicator for the type of parameter that has been 
modified. It is set -1 or +1 for amplitude, -2 or +2 for frequency, and -3 or 
+3 for phase, respectively.

*/

  int index_changed;

/*

This is the line number in the parameters file corresponding to the modified
parameter.

*/

  double reference;

/*

The result of the sum of squares for the current parameter setup is stored 
here.

*/

  double best_improvement;

/*

This variable is for the result of the sum of squares for the new parameter 
setup.

*/

  int i;                        /* index */
  int j;                        /* index */
  int k;                        /* index */

  double zero;

  for (i = 0; i < noisefile_length; i++)
    {
      if (want_log) fprintf(logfile, "\n Initialization...");
      for (j = 0; j < parfile_length; j++) 
	{
	  newamp[i][j] = amp[j];
	  newfreq[i][j] = freq[j];
	  newphase[i][j] = phase[j];
	}
      if (want_log) 
	{
	  fprintf(logfile, "done.\n");
	  fprintf(logfile, "\n  Least squares for noise level %lg:\n", rms[i]);
	}
      whichpar = 0;
      best_improvement = sum_of_squares(i);
      if (want_log) fprintf(logfile, "   RMS %lg.\n", 
	      sqrt(best_improvement / (infile_length - 1)));
      k = 0;
      do 
        {
	  if (want_log) fprintf(logfile, "   Starting iteration %lu.\n", ++k);
          reference = best_improvement;
          switch (whichpar)
	    {
	    case -3:
	      newphase[i][index_changed] -= phase_grid
		* rms[i] / amp[index_changed];
	      if (want_log) 
		{
		  fprintf(logfile, "    Changing phase %lu by %lg.\n",
			  index_changed + 1, -phase_grid
			  * rms[i] / amp[index_changed]);
		  fprintf(logfile, "    New phase: %lg\n", 
			  newphase[i][index_changed]);
		}
	      break;
	    case -2:
	      newfreq[i][index_changed] -= freq_grid * rms[i] / 
		(jd[infile_length - 1] - jd[0]) / amp[index_changed];
	      if (want_log) 
		{
		  fprintf(logfile, "    Changing frequency %lu by %lg.\n",
			  index_changed + 1, 
			  -freq_grid * rms[i] / amp[index_changed] /
			  (jd[infile_length - 1] - jd[0]));
		  fprintf(logfile, "    New frequency: %lg\n", 
			  newfreq[i][index_changed]);
		}
	      break;
	    case -1:
	      newamp[i][index_changed] -= amp_grid * rms[i]; 
	      if (want_log) 
		{
		  fprintf(logfile, "    Changing amplitude %lu by %lg.\n", 
			  index_changed + 1, -amp_grid * rms[i]);
		  fprintf(logfile, "    New amplitude: %lg\n", 
			  newamp[i][index_changed]);
		}
	      break;
	    case 1:
	      newamp[i][index_changed] += amp_grid * rms[i];
	      if (want_log) 
		{
		  fprintf(logfile, "    Changing amplitude %lu by %lg.\n", 
			  index_changed + 1, amp_grid * rms[i]);
		  fprintf(logfile, "    New amplitude: %lg\n", 
			  newamp[i][index_changed]);
		}
	      break;
	    case 2:
	      newfreq[i][index_changed] += freq_grid * rms[i] / 
		(jd[infile_length - 1] - jd[0]) / amp[index_changed];
	      if (want_log) 
		{
		  fprintf(logfile, "    Changing frequency %lu by %lg.\n", 
			  index_changed + 1, 
			  freq_grid * rms[i] / amp[index_changed] /
			  (jd[infile_length - 1] - jd[0]));
		  fprintf(logfile, "    New frequency: %lg\n", 
			  newfreq[i][index_changed]);
		}
	      break;
	    case 3:
	      newphase[i][index_changed] += phase_grid
		* rms[i] / amp[index_changed];
	      if (want_log) 
		{
		  fprintf(logfile, "    Changing phase %lu by %lg.\n", 
			  index_changed + 1, phase_grid
			  * rms[i] / amp[index_changed]);
		  fprintf(logfile, "    New phase: %lg\n", 
			  newphase[i][index_changed]);
		}
	      break;
	    default:
	      break;
	    }
	  for (j = 0; j < parfile_length; j++)
	    {
	      newamp[i][j] -= amp_grid * rms[i];
	      if ((zero = sum_of_squares(i)) < best_improvement)
		{
		  whichpar = -1;
		  index_changed = j;
		  best_improvement = zero;
		  if (want_log) 
		    {
		      fprintf(logfile,
			      "    Improvement found at component %lu.\n",
			      index_changed + 1);
		      fprintf(logfile, "    Preliminary amplitude: %24lg\n",
			      newamp[i][j]);
		      fprintf(logfile, "    Preliminary RMS:       %24lg\n",
			      sqrt(best_improvement / (infile_length - 1)));
		    }
		}
	      newamp[i][j] += 2 * amp_grid * rms[i];
	      if ((zero = sum_of_squares(i)) < best_improvement)
		{
		  whichpar = 1;
		  index_changed = j;
		  best_improvement = zero;
		  if (want_log) 
		    {
		      fprintf(logfile,
			      "    Improvement found at component %lu.\n",
			      index_changed + 1);
		      fprintf(logfile, "    Preliminary amplitude: %24lg\n",
			      newamp[i][j]);
		      fprintf(logfile, "    Preliminary RMS:       %24lg\n",
			      sqrt(best_improvement / (infile_length - 1)));
		    }
		}
	      newamp[i][j] -= amp_grid * rms[i];
	      newfreq[i][j] -= freq_grid * rms[i] / 
		(jd[infile_length - 1] - jd[0]) / amp[i];
	      if ((zero = sum_of_squares(i)) < best_improvement)
		{
		  whichpar = -2;
		  index_changed = j;
		  best_improvement = zero;
		  if (want_log) 
		    {
		      fprintf(logfile,
			      "    Improvement found at component %lu.\n",
			      index_changed + 1);
		      fprintf(logfile, "    Preliminary frequency: %24lg\n",
			      newfreq[i][j]);
		      fprintf(logfile, "    Preliminary RMS:       %24lg\n",
			      sqrt(best_improvement / (infile_length - 1)));
		    }
		}
	      newfreq[i][j] += 2 * freq_grid * rms[i] / 
		(jd[infile_length - 1] - jd[0]) / amp[i];
	      if ((zero = sum_of_squares(i)) < best_improvement)
		{
		  whichpar = 2;
		  index_changed = j;
		  best_improvement = zero;
		  if (want_log) 
		    {
		      fprintf(logfile,
			      "    Improvement found at component %lu.\n",
			      index_changed + 1);
		      fprintf(logfile, "    Preliminary frequency: %24lg\n",
			      newfreq[i][j]);
		      fprintf(logfile, "    Preliminary RMS:       %24lg\n",
			      sqrt(best_improvement / (infile_length - 1)));
		    }
		}
	      newfreq[i][j] -= freq_grid * rms[i] / 
		(jd[infile_length - 1] - jd[0]) / amp[i];
	      newphase[i][j] -= phase_grid * rms[i] / amp[index_changed];
	      if ((zero = sum_of_squares(i)) < best_improvement)
		{
		  whichpar = -3;
		  index_changed = j;
		  best_improvement = zero;
		  if (want_log) 
		    {
		      fprintf(logfile,
			      "    Improvement found at component %lu.\n",
			      index_changed + 1);
		      fprintf(logfile, "    Preliminary phase:     %24lg\n",
			      newphase[i][j]);
		      fprintf(logfile, "    Preliminary RMS:       %24lg\n",
			      sqrt(best_improvement / (infile_length - 1)));
		    }
		}
	      newphase[i][j] += 2 * phase_grid * rms[i] / amp[index_changed];
	      if ((zero = sum_of_squares(i)) < best_improvement)
		{
		  whichpar = 3;
		  index_changed = j;
		  best_improvement = zero;
		  if (want_log) 
		    {
		      fprintf(logfile,
			      "    Improvement found at component %lu.\n",
			      index_changed + 1);
		      fprintf(logfile, "    Preliminary phase:     %24lg\n",
			      newphase[i][j]);
		      fprintf(logfile, "    Preliminary RMS:       %24lg\n",
			      sqrt(best_improvement / (infile_length - 1)));
		    }
		}
	      newphase[i][j] -= phase_grid * rms[i] / amp[index_changed];
	    }
	  if (want_log) for (j = 0; j < parfile_length; j++)
	    {
	      fprintf(logfile, "\n   Result of iteration %8lu:\n", k);
	      fprintf
		(logfile, "   C: %4lu F: %24lg A: %24lg P: %24lg S: %24lg\n",
		 j + 1, newfreq[i][j], newamp[i][j], newphase[i][j],
		 best_improvement);
	    }
	}
      while(best_improvement < reference);
      if (want_log) 
	fprintf(logfile, "\n  Least squares for noise level %lg finished.\n",
	      rms[i]);
    }

  return;
}

/*****************************************************************************/
/********** This function has to attain a minimum for least squares **********/
/*****************************************************************************/

double sum_of_squares(int noiseindex)
{

  int i;                        /* index */
  int j;                        /* index */
  double result;
  double buf;

  result = 0;

  for (i = 0; i < infile_length; i++)
    {
      buf = 0;
      for (j = 0; j < parfile_length; j++)
	  buf += newamp[noiseindex][j] * sin(2 * M_PI * 
	     (newfreq[noiseindex][j] * jd[i] + newphase[noiseindex][j]));
      result += pow(noisy[noiseindex][i] - buf, 2);
    }

  return result;

}