source: util/VibrAlyzer.c@ 2f6525

/* 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, freq_end; // 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;  // 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], &temperatures[counter]);
    //fprintf(stderr, "%lg\t%lg\n", time_steps[counter], temperatures[counter]);
    counter++;
  }
  // auto-set good values for start and step based on step range
  if ((freq_start == -1) || (freq_step == -1)) {
    // we cannot detect frequencies above twice a time step and not below half the time range
    freq_start = 1./(fabs(time_steps[counter-1] - time_steps[0])*gauge);  
    freq_step = ( 1./(fabs(time_steps[1] - time_steps[0])*gauge) - 1./(fabs(time_steps[counter-1] - time_steps[0])*gauge) )/counter;
    fprintf(stderr, "Using %lg and frequency start and %lg as step size.\n", freq_start/(2.*2.*M_PI), freq_step/(2.*2.*M_PI));
  } else {
    freq_start *= (2.*2.*M_PI);
    freq_step *= (2.*2.*M_PI);  
  }

  // for debugging only: print read values
        //for(i=0;i<(counter-1);i++) {
        //      printf("%lg\t%lg\n",time_steps[i],temperatures[i]);
        //}
  
  printf("#frequency(a.u.)\tcos\tsin");
  // discretely integrate over desired frequency range
  freq_end = freq_start+freq_step*counter-1;
        for(frequency = freq_start;frequency < freq_end;frequency += freq_step) {
                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));
        }

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