/* Molecular Vibrations Analyser - VibrAlyzer
 *
 * This programme fourier transforms input from ESPACK (temperature output)
 * in order to make the automated retrieval of vibrational frequencies possible.
*/

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

/** Main routine.
 * The routine needs a file name to be read as the temperature file, and also
 * a frequency range (start and steps). Standard one-dimensional fourier-trans-
 * formation via a simple discrete integration over the given values from the
 * file is performed and the result returned on stdout. 
 * \param argc parameter count
 * \param **argv array of parameter (array of chars)
 * \return error code
*/
int main(int argc, char **argv) 
{
	FILE *temperature_file;		// file with temperature values
	double *time_steps, *temperatures;	// contain data value pairs
	int counter;		// keeps track of data pairs array size
	double freq_start, freq_step;	// frequency start and step width (end determined by number of points in temp.file)
	char *filename;	// filename of temp.file
	char line[255];	// line buffer for parsing the temperature file
	int i,j;	// runtime variable
	double result, iresult;		// temporary result value for fourier transformation
	double frequency;		// current frequency during dumb O(N^2) integration
	double gauge;			// conversion to atomic units for time axis
				
	// Check for needed arguments
	if (argc < 4) {
		printf("Molecular Vibrations Analyser - VibrAlyzer\n\n");
		printf("Usage: %s <time gauge> <freq.step> <temperature file>\n", argv[0]);
		printf("\t<time gauge>\tConversion factor from step count to atomic time\n");
		printf("\t<freq.start>\tstart of frequency for fourier transform in atomic units\n");
		printf("\t<freq.step>\tstep width of frequency for fourier transform in atomic units\n");
		printf("\t<temperature file>\tfile with (time step, temperature)-pairs\n");
		exit(1);
	} else {
	  gauge = atof(argv[1]);
		freq_start = atof(argv[2]);
		freq_step = atof(argv[3]);
		filename = argv[4];
	}

	// read in file into buffer array
	temperature_file=fopen(filename, "r");
	if (temperature_file == NULL) { // check whether file could be opened
		printf("Could not open temperature file named '%s'!\n", filename);
		exit(255);
	}
  // check number of pairs
	counter=0;
	while (fgets(line,255, temperature_file)) {
		sscanf(line,"%lg %lg", &result, &iresult);
		counter++;
	}
  // ... allocate ...
  time_steps = malloc(counter*sizeof(double));
  temperatures = malloc(counter*sizeof(double));
  fclose(temperature_file);
  // ... and parse in
  temperature_file=fopen(filename, "r");
  counter=0;
  while (fgets(line,255, temperature_file)) {
    sscanf(line,"%lg %lg", &time_steps[counter-1], &temperatures[counter-1]);
    counter++;
  }
	// for debugging only: print read values
	//for(i=0;i<(counter-1);i++) {
	//	printf("%lg\t%lg\n",time_steps[i],temperatures[i]);
	//}
  
  // discretely integrate over desired frequency range
	frequency = freq_start;
	for(j=0;j<(counter-1);j++) {
		result = iresult = 0.;
		for(i=0;i<(counter-1);i++) {
			result += temperatures[i] * cos(frequency * time_steps[i]*gauge);
			iresult += temperatures[i] * sin(frequency * time_steps[i]*gauge);
		}
    // NOTE: It is by definition freq over 2*pi, however as the temperature curve counts double (there are two
    // standstills, one at the perihel one at the aphel!) we insert this factor to make the plots automatically
    // have the correct frequency! 
		printf("%lg\t%lg\t%lg\n",frequency/(2.*2.*M_PI),result/(counter-1),iresult/(counter-1));
		frequency += freq_step;
	}

	// dis'alloc and end
	exit(0);
}