/*
** Known bugs and things to do :-) :
**
**	 - Left/right scrolling is buggy - especially after some zoom operations... *sigh*
**	 - scale_factor lower bound logic does not work properly!
**	 - z-function to reset current parameter set
**     	 - Use date_string in create_daily_plot to denote traces
**	(- implement remove_flaws)
**	 - Implement ticks on coordinate systems (esp. for spectrum data!) to enhace readability
**	 - Interactive with zero cut off frequency results in a crash!
*/

/*
** The seismic_toolkit, stk for short, is a small C program which is quite useful for producing hourly and daily plots, 
** perform low pass filter operations and to display measured data, if desired.
**
** 25/26-FEB-2003, B. Ulmann fecit.
*/

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

/*
** Local includes
*/

#include "misc.h"                                      /* String length settings, etc. */
#include "read_ini_file.h"                             /* Function prototypes for the read_ini_file routine collection */
#include "rfftw.h"                                     /* Include file for real FFT transformation package */

#ifdef VMS
# include "gltk.h"                                     /* Include file for Mesa graphics library */
#else
# include "gl.h"
# include "GL/glut.h"
#endif

/*
** Miscellaneous constant definitions
*/

#define INI_FILE_NAME "stk.ini"                        /* Name of the initialisation file to be read during startup */
#define MAX_INPUT_FILES 24                             /* Maximum number of input files to be processed */
#define MAX_HOURLY_CHANNELS 8                          /* Maximum number of channels to be displayed for hourly plots */
#define PROGRAM_NAME "Seismic toolkit"                 /* Name of the running program for use as window title */

#define LINE_LENGTH 37                                 /* Line length for postscript plots */
#define X0 3                                           /* X start offset for daily postscript plots */
#define Y0 0                                           /* Y start offset for daily postscript plots */
#define DY_DAILY .85
#define DY_HOURLY 2.5
#define EPSILON .001

#define STK$ACTION_INTERACTIVE 1                       /* Possible action codes */
#define STK$ACTION_FILTER      2
#define STK$ACTION_DAILY       3
#define STK$ACTION_HOURLY      4

#define DEBUG

/*
** Typedefs
*/

typedef struct                                         /* Structure which holds all necessary global control parameters */
{
  unsigned int first_line_contains_date,               /* TRUE -> The first line of each input file contains the start date */
    samples_per_second,                                /* Samples per second in each input file */
    seconds_per_file,                                  /* Seconds of data in each input file */
    action,                                            /* Selected action to be performed */
    channel,                                           /* Number of the channel to be processed */
    number_of_input_files,                             /* How many input files are to be read */
    number_of_samples,                                 /* Total number of samples to be processed */
    number_of_channels,                                /* Number of channels to be displayed in an hourly plot */
    verbose,                                           /* TRUE -> print verbose messages during program run */
    daily_skip_factor,                                 /* Determines how many samples are to be skipped for a daily display */
    hourly_skip_factor,                                /* Determines how many samples are to be skipped for an hourly display */
    max_samples_per_screen,                            /* Maximum number of samples to be displayed in the output window */
    tune,                                              /* Possible values: FFTW_ESTIMATE, FFTW_MEASURE */
    use_wisdom,                                        /* Use the wisdom mechanism of the (R)FFTW routines */
    curve_height,                                      /* Y factor for interactive graphics display of data */
    display_coordinate_systems,                        /* TRUE -> Show coordinate systems in the output plot */
    scale_filtered_like_raw,                           /* TRUE -> Use the same scaling factors for raw and filtered data */
    window_x_position,                                 /* Initial x position of output window */
    window_y_position,                                 /* Initial y position of output window */
    window_x_size,                                     /* Initial x size of output window */
    window_y_size;                                     /* Initial y size of output window */

  int raw_data_offset,                                 /* Y offset for the display of raw data */
    filtered_data_offset,                              /* Y offset for the display of filtered data */
    spectrum_data_offset;                              /* Y offset for the display of spectrum data */

  double low_pass_frequency,                           /* Low pass frequency for the FFT filter */
    hourly_low_pass_frequency [MAX_HOURLY_CHANNELS],   /* Cutoff frequencies for hourly plots */
    daily_y_factor,                                    /* Y axis magnifictation for daily plots */
    scale_factor,                                      /* Scale factor for display of graphics */
    hourly_y_factor [MAX_HOURLY_CHANNELS],                    /* Y scale factors for hourly plots */
    initial_scale_factor,                              /* Inital scale factor, used for scaling the graphics */
    scale_factor_multiplier,                           /* Factor to multiply the scale factor with during zoom operations */
    frequency_increment,                               /* Change to f_cutoff for interactive display */
    x_translate,                                       /* Actual X translation of curve being displayed */
    y_translate,                                       /* Actual Y translation of curve being displayed */
    z_translate,                                       /* Actual Z translation of curve being displayed */
    x_translate_increment,                             /* Increment factor for X translation */
    y_translate_increment;                             /* Increment factor for X translation */

  char output_file_name [STRING_LENGTH],               /* Name of the output file */
    daily_title [STRING_LENGTH],                       /* Title which is to appear on the daily plots */
    hourly_title [STRING_LENGTH],                      /* Title which is to appear on the hourly plots */
    hourly_channel_description [STRING_LENGTH],        /* Channel description string for hourly plots */
    wisdom_file_name [STRING_LENGTH],                  /* Name of the wisdom file to be used by the (R)FFTW routines */
    input_file_name [MAX_INPUT_FILES] [STRING_LENGTH]; /* Array of pointers to the names of all input file names */

  rfftw_plan plan_real_to_complex,                     /* Plan for all RFFTW routines used in forward direction */
    plan_complex_to_real;                              /* Plan for all RFFTW routines used in backward direction */

  GLuint graphics_object;                              /* Mesa graphics library graphics object */
} control_data;

typedef struct
{
  fftw_real *data;                                     /* Pointer to array containing the raw data */
  double max_value,                                    /* Maximum value of the data read into the data array */
    average;                                           /* Average value of the data read into the data array */
} rfftw_data;

/*
** Global variables
*/

control_data gbl$control_data;                         /* Global data structure holding all necessary control values */

rfftw_data gbl$raw_data,                               /* Global structure holding the raw data and its max/avg values */
  gbl$filtered_data,                                   /* Global structure holding the filtered data and its max/avg values */
  gbl$spectrum_data;                                   /* Global structure holding the necessary data for a power spectrum */

/*
** show_graphics actually shows a graphics display previously created by a call to make_graphics_object.
*/

void show_graphics ()
{
#ifdef DEBUG
  printf ("show_graphics ()\n");
#endif

  glClear (GL_COLOR_BUFFER_BIT);

  glPushMatrix ();
  gluLookAt (gbl$control_data.x_translate, gbl$control_data.y_translate, gbl$control_data.z_translate,
             gbl$control_data.x_translate, gbl$control_data.y_translate, 0.,
             0., 1., 0.);
  glScalef (gbl$control_data.scale_factor, gbl$control_data.scale_factor, gbl$control_data.scale_factor);
  glCallList (gbl$control_data.graphics_object);
  glPopMatrix ();

  glFlush ();
}

/*
** reshape is used for all window operations incorporating a reshaping of the display area.
*/

static void reshape (int width, int height)
{
#ifdef DEBUG
  printf ("reshape (width = %d, height = %d)\n", width, height);
#endif

  glViewport (0, 0, (GLint) width, (GLint) height);
  glMatrixMode (GL_PROJECTION);
  glLoadIdentity ();
  glFrustum (-1., 1., -1., 1., 5., 15.);
  glMatrixMode (GL_MODELVIEW);
}

/*
** make_graphics_object creates the necessary data structures for a complete curve display which can be shown with show_graphics.
*/

static GLuint make_graphics_object ()
{
  unsigned int i,
    skip_factor_display;

  double x_stretch,
    y_stretch;
  GLuint list;

#ifdef DEBUG
  printf ("make_graphics_object ()\n");
#endif

  list = glGenLists (1);
  glNewList (list, GL_COMPILE);

  skip_factor_display = gbl$control_data.number_of_samples / gbl$control_data.max_samples_per_screen;

  /* Plot raw data */

  glBegin (GL_LINE_STRIP);
  glColor3f (0., 1., 1.);

  x_stretch = (double) gbl$control_data.window_x_size / (double) gbl$control_data.number_of_samples;
  y_stretch = gbl$control_data.curve_height / gbl$raw_data.max_value;

  for (i = 0; i < gbl$control_data.number_of_samples; i += skip_factor_display)
    glVertex3f ((double) i * x_stretch, (gbl$raw_data.data [i] - gbl$raw_data.average) * y_stretch + 
                                        gbl$control_data.raw_data_offset, 0.);

  glEnd ();

  /* Plot power spectrum of the raw data */

  glBegin (GL_LINE_STRIP);
  glColor3f (1., 0., 1.);

  x_stretch = (double) gbl$control_data.window_x_size / (double) (gbl$control_data.number_of_samples >> 1);
  y_stretch = gbl$control_data.curve_height / gbl$spectrum_data.max_value;

  for (i = 0; i < (gbl$control_data.number_of_samples >> 1); i += (skip_factor_display << 1))
    glVertex3f ((double) i * x_stretch, (gbl$spectrum_data.data [i] - gbl$spectrum_data.average) * y_stretch + 
                                        gbl$control_data.spectrum_data_offset, 0.);

  glEnd ();

  /* Plot filtered data */

  glBegin (GL_LINE_STRIP);
  glColor3f (1., 1., 0.);

  x_stretch = (double) gbl$control_data.window_x_size / (double) gbl$control_data.number_of_samples;
  if (gbl$control_data.scale_filtered_like_raw)
    y_stretch = gbl$control_data.curve_height / gbl$raw_data.max_value; 
  else
    y_stretch = gbl$control_data.curve_height / gbl$filtered_data.max_value;

  for (i = 0; i < gbl$control_data.number_of_samples; i += skip_factor_display)
    glVertex3f ((double) i * x_stretch, (gbl$filtered_data.data [i] - gbl$filtered_data.average) * y_stretch + 
                                        gbl$control_data.filtered_data_offset, 0.);

  glEnd ();

  /* Plot coordinate systems if desired */

  if (gbl$control_data.display_coordinate_systems)
  {
    /* Plot coordinate system for raw data */

    glBegin (GL_LINE_STRIP);
    glColor3f (1., 1., 1.);

    glVertex3f (0., (double) gbl$control_data.raw_data_offset, 0.);
    glVertex3f ((double) gbl$control_data.window_x_size, (double) gbl$control_data.raw_data_offset, 0.);

    glEnd ();

    /* Plot coordinate system for filtered data */

    glBegin (GL_LINE_STRIP);
    glColor3f (1., 1., 1.);

    glVertex3f (0., (double) gbl$control_data.filtered_data_offset, 0.);
    glVertex3f ((double) gbl$control_data.window_x_size, (double) gbl$control_data.filtered_data_offset, 0.);

    glEnd ();

    /* Plot coordinate system for spectrum data */

    glBegin (GL_LINE_STRIP);
    glColor3f (1., 1., 1.);
 
    glVertex3f (0., (double) gbl$control_data.spectrum_data_offset, 0.);
    glVertex3f ((double) gbl$control_data.window_x_size, (double) gbl$control_data.spectrum_data_offset, 0.);

    glEnd ();
  }

  glEndList ();

  return list;
}

/*
** write_filtered_data writes the actual filtered data set out to the desired output file. 
*/

void write_filtered_data ()
{
  unsigned int i;

  FILE *handle;

#ifdef DEBUG
  printf ("write_filtered_data ()\n");
#endif

  if (!(handle = fopen (gbl$control_data.output_file_name, "w")))
    printf ("write_filtered_data> Unable to open output file >>%s<< for filtered data!\n", gbl$control_data.output_file_name);
  else
  { 
    if (gbl$control_data.verbose)
      printf ("write_filtered_data> Write filtered data to output file >>%s<<.\n", gbl$control_data.output_file_name);

    for (i = 0; i < gbl$control_data.number_of_samples; i++)
      fprintf (handle, "%g\n", gbl$filtered_data.data [i]);

    fclose (handle);
  }
}

/*
** low_pass implements an FFT based low pass filter for real data using the RFFTW libraries.
*/

void low_pass (unsigned int number_of_samples, double low_pass_frequency)
{
  unsigned int count, 
    i;

  double value,
    sum;

#ifdef DEBUG
  printf ("low_pass ()\n");
#endif

  if (low_pass_frequency <= 0.)
  {
    printf ("low_pass> Function was called with a cut off frequency of 0Hz -> skip!\n");
    return;
  }

  if (gbl$control_data.verbose)
    printf ("low_pass> Perform forward FFT transformation.\n");

  rfftw_one (gbl$control_data.plan_real_to_complex, gbl$raw_data.data, gbl$spectrum_data.data);

  count =  number_of_samples * low_pass_frequency / gbl$control_data.samples_per_second;

  if (gbl$control_data.verbose)
    printf ("          Apply low pass filter with cutoff freqency %f Hz.\n\
          Keep %d coefficients.\n", low_pass_frequency, count);

  for (i = 0; i < number_of_samples; i++)
    if ((i >= count) && (i <= number_of_samples - count))
      gbl$spectrum_data.data [i] = (fftw_real) 0.;

  if (gbl$control_data.verbose)
    printf ("          Perform backward FFT transformation.\n");

  rfftw_one (gbl$control_data.plan_complex_to_real, gbl$spectrum_data.data, gbl$filtered_data.data);

  if (gbl$control_data.verbose)
    printf ("          Perform normalisation of filtered output data.\n");

  gbl$filtered_data.max_value = gbl$filtered_data.average = 0.;
  for (i = 0; i < number_of_samples; i++)
  {
    value = gbl$filtered_data.data [i] / (double) number_of_samples;
    sum += (gbl$filtered_data.data [i] = value);
    if (fabs (value) > gbl$filtered_data.max_value)
      gbl$filtered_data.max_value = fabs (value);
  }
  gbl$filtered_data.average = sum / (double) number_of_samples;

  if (gbl$control_data.verbose)
    printf ("low_pass> Maximum value: %f\n\
          Average value: %f\n", gbl$filtered_data.max_value, gbl$filtered_data.average);
}

/*
** power_spectrum calculates the power spectrum of a previously transformed raw data set.
*/

void power_spectrum ()
{
  unsigned int i;

  double sum;

#ifdef DEBUG
  printf ("power_spectrum ()\n");
#endif

  /* Recalculate the FFT transform of the raw data since it has been destroyed by the previous back transformation */

  rfftw_one (gbl$control_data.plan_real_to_complex, gbl$raw_data.data, gbl$spectrum_data.data);

  gbl$spectrum_data.data [0] *= gbl$spectrum_data.data [0]; /* DC component */

  /* Do not use the DC component since it will be too large for display */

  sum = gbl$spectrum_data.max_value = gbl$spectrum_data.average = 0.;

  for (i = 1; i < (gbl$control_data.number_of_samples >> 1); i++)
  {
    if (fabs (gbl$spectrum_data.data [i] = gbl$spectrum_data.data [i] * gbl$spectrum_data.data [i] +
                                           gbl$spectrum_data.data [gbl$control_data.number_of_samples - i - 1] * 
                                           gbl$spectrum_data.data [gbl$control_data.number_of_samples - i - 1]) 
        > gbl$spectrum_data.max_value)
      gbl$spectrum_data.max_value = fabs (gbl$spectrum_data.data [i]);

    sum += gbl$spectrum_data.data [i];
  }

  if (!(gbl$control_data.number_of_samples % 2))                        /* Even number of samples */
  {
    gbl$spectrum_data.data [gbl$control_data.number_of_samples >> 1] *= 
      gbl$spectrum_data.data [gbl$control_data.number_of_samples >> 1]; /* Nyquist frequency */

    if (fabs (gbl$spectrum_data.data [gbl$control_data.number_of_samples >> 1]) > gbl$spectrum_data.max_value)
      gbl$spectrum_data.max_value = fabs (gbl$spectrum_data.data [gbl$control_data.number_of_samples >> 1]);
    sum += gbl$spectrum_data.data [gbl$control_data.number_of_samples >> 1];
  }

  gbl$spectrum_data.average = sum / (gbl$control_data.number_of_samples >> 1);

  if (gbl$control_data.verbose)
    printf ("power_spectrum> Maximum value: %f\n\
                Average value: %f\n", gbl$spectrum_data.max_value, gbl$spectrum_data.average);
}

/*
** This function performs all filtering calculations and the calculation of the power spectrum which are needed for the
** -interactive and the -filter mode. It was necessary to get this out of the function filter () to be able to control
** it by the key events 'I' and 'D'.
*/

void filter_and_spectrum_calc ()
{
#ifdef DEBUG
  printf ("filter_and_spectrum_calc ()\n");
#endif

  if (gbl$control_data.low_pass_frequency != 0.)
    low_pass (gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file * gbl$control_data.number_of_input_files,
              gbl$control_data.low_pass_frequency);

  if (gbl$control_data.action == STK$ACTION_INTERACTIVE)
    power_spectrum ();
}

/*
** key_event is the main routine to process key events signalized to the tk-system.
*/

#ifdef VMS
static GLenum key_event (int key, GLenum mask)
#else
static void key_event (unsigned char key, int x, int y)
#endif
{
#ifdef DEBUG
# ifdef VMS
  printf ("key_event (key = %d)\n", key);
# else
  printf ("key_event (key = %d, x = %d, y = %d)\n", (int) key, x, y);
# endif
#endif

  switch (key)
  {
    case 'q':
#ifdef VMS
      tkQuit ();
#else
      exit (0);
#endif
      break;
    case 'r':           /* Shift right */
      gbl$control_data.x_translate += gbl$control_data.x_translate_increment;
      show_graphics ();
      break;
    case 'l':           /* Shift left */
      gbl$control_data.x_translate -= gbl$control_data.x_translate_increment;
      show_graphics ();
      break;
    case 'u':	        /* Shift up */
      gbl$control_data.y_translate += gbl$control_data.y_translate_increment;
      show_graphics ();
      break;
    case 'd':           /* Shift down */
      gbl$control_data.y_translate -= gbl$control_data.y_translate_increment;
      show_graphics ();
      break;
    case 'h':           /* Print help */
      printf ("\n\n\
+------------------------------+\n\
! l........Shift left          !\n\
! r........Shift right         !\n\
! u........Shift up            !\n\
! d........Shift down          !\n\
! h........Print this help     !\n\
! L........Larger              !\n\
! q........Quit                !\n\
! s........Smaller             !\n\
! I........Increment f_cutoff  !\n\
! D........Decrement f_cutoff  !\n\
! w........Write filtered data !\n\
+------------------------------+\n\
\n");
      break;
    case 'I': /* Increment cut off frequency */
      gbl$control_data.low_pass_frequency += gbl$control_data.frequency_increment;
      printf ("key_event> New low pass frequency: %f\n", gbl$control_data.low_pass_frequency);
      filter_and_spectrum_calc ();
      gbl$control_data.graphics_object = make_graphics_object ();
      show_graphics ();
      break;
    case 'D': /* Decrement cut off frequency */
      gbl$control_data.low_pass_frequency -= gbl$control_data.frequency_increment;
      if (gbl$control_data.low_pass_frequency <= EPSILON)
        gbl$control_data.low_pass_frequency = 0.;
      else
      {
        printf ("key_event> New low pass frequency: %f\n", gbl$control_data.low_pass_frequency);
        filter_and_spectrum_calc ();
        gbl$control_data.graphics_object = make_graphics_object ();
        show_graphics ();
      }
      break;
    case 'L':           /* Enlarge view */
      gbl$control_data.scale_factor *= gbl$control_data.scale_factor_multiplier;
      show_graphics ();
      break;
    case 'p':           /* Print out current parameter set */
      printf ("\n\n\
+----------------------------+\n\
! x_translate....%f           \n\
! y_translate....%f           \n\
! scale_factor...%f           \n\
+----------------------------+\n\
\n",
              gbl$control_data.x_translate, gbl$control_data.y_translate, gbl$control_data.scale_factor);
      break;
    case 's':           /* Make view smaller */
      gbl$control_data.scale_factor /= gbl$control_data.scale_factor_multiplier;
      if (gbl$control_data.scale_factor == 0)
      {
        printf ("key_event> scale_factor was set to zero. Reset to %f.\n", gbl$control_data.initial_scale_factor);
        gbl$control_data.scale_factor = gbl$control_data.initial_scale_factor;
      }
      show_graphics ();
      break;
    case 'w':
      write_filtered_data ();
      break;
    default:
      printf ("key_event> Unknown key event, skip event!\n");
      break;
  }

#ifdef VMS
  return GL_FALSE;
#endif
}

/*
** stk_create_plan creates actually two plans for the RFFTW library routines: One plan for forward FFTs and one for backward FFTs.
*/

void stk_create_plan (unsigned int number_of_samples)
{
  unsigned int wisdom_ok;

  FILE *wisdom_file;

#ifdef DEBUG
  printf ("stk_create_plan (number_of_samples = %d)\n", number_of_samples);
#endif

  if (!gbl$control_data.use_wisdom) /* Do not use wisdom */
  {
    if (gbl$control_data.verbose)
      printf ("stk_create_plan> Create forward plan without using wisdom.\n");

    gbl$control_data.plan_real_to_complex = rfftw_create_plan (number_of_samples, FFTW_REAL_TO_COMPLEX, gbl$control_data.tune);

    if (gbl$control_data.verbose)
      printf ("stk_create_plan> Create backward plan without using wisdom.\n");

    gbl$control_data.plan_complex_to_real = rfftw_create_plan (number_of_samples, FFTW_COMPLEX_TO_REAL, gbl$control_data.tune);
  }
  else                              /* Use wisdom to speed up plan creation */
  {
    wisdom_ok = TRUE;
    wisdom_file = fopen (gbl$control_data.wisdom_file_name, "r");
    if (fftw_import_wisdom_from_file (wisdom_file) == FFTW_FAILURE)
    {
      wisdom_ok = FALSE;
      printf ("stk_create_plan> Error reading wisdom file >>%s<<!\n", gbl$control_data.wisdom_file_name);
    }
    fclose (wisdom_file);

    if (gbl$control_data.verbose)
      printf ("stk_create_plan> Create forward plan with using wisdom.\n");

    gbl$control_data.plan_real_to_complex = rfftw_create_plan (number_of_samples, FFTW_REAL_TO_COMPLEX, 
                                                               gbl$control_data.tune || FFTW_USE_WISDOM);
    if (gbl$control_data.verbose)
      printf ("stk_create_plan> Create backward plan with using wisdom.\n");

    gbl$control_data.plan_complex_to_real = rfftw_create_plan (number_of_samples, FFTW_COMPLEX_TO_REAL, 
                                                               gbl$control_data.tune || FFTW_USE_WISDOM);

    if (!wisdom_ok) /* The wisdom file could not be read so try creating one now */
    {
      if (gbl$control_data.verbose)
        printf ("stk_create_plan> Write new wisdom file >>%s<<.\n", gbl$control_data.wisdom_file_name);

      if (!(wisdom_file = fopen (gbl$control_data.wisdom_file_name, "w")))
        printf ("stk_create_plan> Unable to create new wisdom file >>%s<<!\n", gbl$control_data.wisdom_file_name);
      else
      {
        fftw_export_wisdom_to_file (wisdom_file);
        fclose (wisdom_file);
      }
    }
  }
}

/*
** allocate_memory allocates three buffers of memory - one for raw, filtered and spectrum data each. number_of_buffers contains 
** the number of files to be processed simultaneously.
*/

int allocate_memory (unsigned int number_of_buffers)
{
  unsigned int memory_needed;

#ifdef DEBUG
  printf ("allocate_memory (number_of_buffers = %d)\n", number_of_buffers);
#endif

  if (gbl$control_data.verbose)
    printf ("allocate_memory> Allocate memory for %d buffers.\n", number_of_buffers);

  memory_needed = gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file * number_of_buffers * sizeof (fftw_real);

  if (!(gbl$raw_data.data = malloc (memory_needed)))
  {
    printf ("allocate_memory> Unable to allocate %d bytes of memory for input data array!\n", memory_needed);
    return FALSE;
  }

  if (!(gbl$filtered_data.data = malloc (memory_needed)))
  {
    printf ("allocate_memory> Unable to allocate %d bytes of memory for output data array!\n", memory_needed);
    return FALSE;
  }

  if (!(gbl$spectrum_data.data = malloc (memory_needed)))
  {
    printf ("allocate_memory> Unable to allocate %d bytes of memory for spectrum data array!\n", memory_needed);
    return FALSE;
  }

  return TRUE;
}

/*
** create_daily_plot is the main routine which performs all necessary actions to create a postscript file containing 
** a complete daily plot of seismic data after applying a low pass filter (if desired) to each input file.
*/

void create_daily_plot ()
{
  unsigned int file_nr,
    sample_nr,
    j;

  double x_pos,
    y0,
    dx,
    value,
    sum;

  char date_string [STRING_LENGTH],
    line [STRING_LENGTH],
    *row;

  FILE *output_handle,
    *input_handle;

#ifdef DEBUG
  printf ("create_daily_plot ()\n");
#endif

  /* Perform initial checks */

  if (gbl$control_data.samples_per_second % gbl$control_data.daily_skip_factor)
  {
    printf ("create_daily_plot> [daily_plot]skip_factor must be a divisor of [file]samples_per_second!\n");
    return;
  }

  if (gbl$control_data.low_pass_frequency != 0)
    stk_create_plan (gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file);

  *date_string = (char) 0;

  if (gbl$control_data.verbose)
    printf ("create_daily_plot> Open output file >>%s<<\n", gbl$control_data.output_file_name);

  if (!(output_handle = fopen (gbl$control_data.output_file_name, "w")))
  {
    printf ("create_daily_plot> Unable to open output file >>%s<<!\n", gbl$control_data.output_file_name);
    return;
  }

  /* Allocate memory for the fourier transform */

  if (!allocate_memory (1))
  {
    printf ("create_daily_plot> Unable to allocate memory! Abort!\n");
    fclose (output_handle);
    return;
  }

  /* Create postscript preliminaries */

  fprintf (output_handle, "\
72 2.54 div dup scale\n\
.01 setlinewidth\n\
-90 rotate\n\
-25 0 translate\n\
90 rotate\n\
0 -3 translate\n\
0.05 setlinewidth\n\
/Times-Roman findfont\n\
1.5 scalefont setfont\n\
newpath 3 2 moveto (%s) show\n\
0.015 setlinewidth\n\
/Times-Roman findfont\n\
.75 scalefont setfont\n\
newpath 3 1 moveto (Sample rate: %d, %d lines of \n\
%d seconds of data each.) show\n\
-90 rotate\n\
0.01 setlinewidth\n\
newpath\n",
           gbl$control_data.daily_title,
           gbl$control_data.samples_per_second / gbl$control_data.daily_skip_factor, gbl$control_data.number_of_input_files, 
           gbl$control_data.seconds_per_file);

  dx = (double) LINE_LENGTH / 
       (double) (gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file / gbl$control_data.daily_skip_factor);

  /* Draw graphs */

  for (file_nr = 0; file_nr < gbl$control_data.number_of_input_files; file_nr++)
  {
    if (gbl$control_data.verbose)
      printf ("create_daily_plot> Process input file >>%s<<.\n", gbl$control_data.input_file_name [file_nr]);

    if (!(input_handle = fopen (gbl$control_data.input_file_name [file_nr], "r")))
    {
      printf ("create_daily_plot> Unable to open input file >>%s<<!\n", gbl$control_data.input_file_name [file_nr]);
      fclose (output_handle);

      free (gbl$raw_data.data);
      free (gbl$spectrum_data.data);
      free (gbl$filtered_data.data);

      return;
    }

    if (gbl$control_data.first_line_contains_date) /* Skip first line and remember the start date */
      fgets (date_string, STRING_LENGTH, input_handle);

    gbl$raw_data.max_value = sum = 0.;
    for (sample_nr = 0; sample_nr < gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file; sample_nr++)
    {
      if (feof (input_handle))
      {
        printf ("create_daily_plot> Unexpected end of file - abort!\n");

        fclose (input_handle);
        fclose (output_handle);

        free (gbl$raw_data.data);
        free (gbl$spectrum_data.data);
        free (gbl$filtered_data.data);

        return;
      }

      fgets (line, STRING_LENGTH, input_handle);

      row = strtok (line, " ");                             /* Get first token from line */
      for (j = 0; row && j < gbl$control_data.channel; j++) /* Skip channels if desired */
        row = strtok (0, " ");

      if (!row)
      {
        printf ("create_daily_plot> Unexpected end of line - channel number too large? Abort!\n");

        fclose (input_handle);
        fclose (output_handle);

        free (gbl$raw_data.data);
        free (gbl$spectrum_data.data);
        free (gbl$filtered_data.data);

        return;
      }

      sscanf (row, "%f", &value);
      sum += value;
      if (fabs (gbl$raw_data.data [sample_nr] = value) > gbl$raw_data.max_value) /* Determine new maximum */
        gbl$raw_data.max_value = fabs (value);
    }
    gbl$raw_data.average = sum / (double) sample_nr;

    fclose (input_handle);

    if (gbl$control_data.verbose)
      printf ("                   Maximum value: %f\n\
                   Average value: %f\n", gbl$raw_data.max_value, gbl$raw_data.average);

    if (gbl$control_data.low_pass_frequency != 0)
      low_pass (gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file, gbl$control_data.low_pass_frequency);

    fprintf (output_handle, "%d %d moveto\n", 
             (int) (y0 = (double) (Y0 + file_nr * DY_DAILY)), X0);

    for (sample_nr = 0, x_pos = (double) X0 + dx; 
         sample_nr < gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file;
         sample_nr += gbl$control_data.daily_skip_factor, x_pos += dx)
      if (gbl$control_data.low_pass_frequency != 0.)
        fprintf (output_handle, "%f %f lineto\n", 
                 y0 + (gbl$filtered_data.data [sample_nr] - gbl$raw_data.average) * (double) gbl$control_data.daily_y_factor, 
                 x_pos);
      else
        fprintf (output_handle, "%f %f lineto\n", 
                 y0 + (gbl$raw_data.data [sample_nr] - gbl$raw_data.average) * (double) gbl$control_data.daily_y_factor, x_pos);

    if (gbl$control_data.verbose)
      printf ("                   Line %d plotted.\n", file_nr);
  }

  /* Close output file */

  fprintf (output_handle, "stroke\nshowpage\n");
  fclose (output_handle);

  free (gbl$raw_data.data);
  free (gbl$spectrum_data.data);
  free (gbl$filtered_data.data);
}

/*
** create_hourly_plot is the main routine which performs all necessary actions to create a postscript file containing 
** a complete hourly plot of seismic data.
*/

void create_hourly_plot ()
{
  unsigned int j,
    channel,
    sample_nr;

  double sum,
    value,
    y0,
    x_pos,
    dx;

  char date_string [STRING_LENGTH],
    line [STRING_LENGTH],
    *row;

  FILE *input_handle,
    *output_handle;

#ifdef DEBUG
  printf ("create_hourly_plot ()\n");
#endif

  /* Perform initial checks */

  if (gbl$control_data.samples_per_second % gbl$control_data.hourly_skip_factor)
  {
    printf ("create_hourly_plot> [hourly_plot]skip_factor must be a divisor of [file]samples_per_second!\n");
    return;
  }

  stk_create_plan (gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file);
  *date_string = (char) 0;

  if (gbl$control_data.verbose)
    printf ("create_hourly_plot> Open output file >>%s<<\n", gbl$control_data.output_file_name);

  if (!(output_handle = fopen (gbl$control_data.output_file_name, "w")))
  {
    printf ("create_hourly_plot> Unable to open output file >>%s<<!\n", gbl$control_data.output_file_name);
    return;
  }

  if (!allocate_memory (1)) /* Try to allocate memory for spectrum and filtered data */
  {
    printf ("create_hourly_plot> Unable to allocate FFTW buffers! Abort!\n");
    return;
  }

  /* Open input file to get the first line, if necessary */

  if (gbl$control_data.first_line_contains_date) /* Skip first line */
  {
    if (!(input_handle = fopen (gbl$control_data.input_file_name [0], "r")))
    {
      printf ("create_hourly_plot> Unable to open input file for reading datetime. Abort!\n");
      fclose (output_handle);

      free (gbl$raw_data.data);
      free (gbl$spectrum_data.data);
      free (gbl$filtered_data.data);

      return;
    }
    fgets (date_string, STRING_LENGTH, input_handle);
    fclose (input_handle);
  }

  /* Create postscript preliminaries */

  fprintf (output_handle, "\
72 2.54 div dup scale\n\
.01 setlinewidth\n\
-90 rotate\n\
-25 0 translate\n\
90 rotate\n\
0 -3 translate\n\
0.05 setlinewidth\n\
/Times-Roman findfont\n\
1.5 scalefont setfont\n\
newpath 3 2 moveto (%s %s) show\n\
0.015 setlinewidth\n\
/Times-Roman findfont\n\
.75 scalefont setfont\n\
newpath 3 1 moveto (Sample rate: %d, %d channels of \n\
%d seconds of data each.) show\n\
newpath 3 0 moveto (%s) show\n\
-90 rotate\n\
0.01 setlinewidth\n\
newpath\n",
           gbl$control_data.hourly_title, date_string,
           gbl$control_data.samples_per_second / gbl$control_data.hourly_skip_factor, gbl$control_data.number_of_channels, 
           gbl$control_data.seconds_per_file,
           gbl$control_data.hourly_channel_description);

  dx = (double) LINE_LENGTH / 
       (double) (gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file / gbl$control_data.hourly_skip_factor);

  /* To save memory, just loop gbl$control_data.number_of_channel times over the input file and process each channel */

  for (channel = 0; channel < gbl$control_data.number_of_channels; channel++)
  {
    if (gbl$control_data.verbose)
      printf ("create_hourly_plot> Read channel %d.\n", channel);

    if (!(input_handle = fopen (gbl$control_data.input_file_name [0], "r")))
    {
      printf ("create_hourly_plot> Unable to open input file for reading channel %d. Abort!\n", channel);
      fclose (output_handle);

      free (gbl$raw_data.data);
      free (gbl$spectrum_data.data);
      free (gbl$filtered_data.data);

      return;
    }

    if (gbl$control_data.first_line_contains_date) /* Skip first line */
      fgets (date_string, STRING_LENGTH, input_handle);

    gbl$raw_data.max_value = sum = 0.;
    for (sample_nr = 0; sample_nr < gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file; sample_nr++)
    {
      if (feof (input_handle))
      {
        printf ("create_hourly_plot> Unexpected end of file - abort!\n");

        fclose (input_handle);
        fclose (output_handle);

        free (gbl$raw_data.data);
        free (gbl$spectrum_data.data);
        free (gbl$filtered_data.data);

        return;
      }

      fgets (line, STRING_LENGTH, input_handle);

      row = strtok (line, " ");            /* Get first token from line */
      for (j = 0; row && j < channel; j++) /* Skip channels if desired */
        row = strtok (0, " ");

      if (!row)
      {
        printf ("create_hourly_plot> Unexpected end of line - channel number too large? Abort at sample number %d!\n", 
                sample_nr);

        fclose (input_handle);
        fclose (output_handle);

        free (gbl$raw_data.data);
        free (gbl$spectrum_data.data);
        free (gbl$filtered_data.data);

        return;
      }

      sscanf (row, "%f", &value);
      sum += value;
      if (fabs (gbl$raw_data.data [sample_nr] = value) > gbl$raw_data.max_value) /* Determine new maximum */
        gbl$raw_data.max_value = fabs (value);
    }
    gbl$raw_data.average = sum / (double) sample_nr;

    fclose (input_handle);

    if (gbl$control_data.verbose)
      printf ("                    Maximum value: %f\n\
                    Average value: %f\n", gbl$raw_data.max_value, gbl$raw_data.average);

    if (gbl$control_data.hourly_low_pass_frequency [channel] != 0.)
    {
      if (gbl$control_data.verbose)
        printf ("create_hourly_plot> Apply low pass filter to channel %d.\n", channel);

      low_pass (gbl$control_data.samples_per_second *
                gbl$control_data.seconds_per_file, gbl$control_data.hourly_low_pass_frequency [channel]);
    }

    fprintf (output_handle, "%d %d moveto\n", 
             (int) (y0 = (double) (Y0 + (channel + 2) * DY_HOURLY)), X0);

    for (sample_nr = 0, x_pos = (double) X0 + dx; 
         sample_nr < gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file;
         sample_nr += gbl$control_data.hourly_skip_factor, x_pos += dx)
      if (gbl$control_data.hourly_low_pass_frequency [channel] != 0.)
        fprintf (output_handle, "%f %f lineto\n", 
                 y0 + (gbl$filtered_data.data [sample_nr] - gbl$raw_data.average) * 
                 (double) gbl$control_data.hourly_y_factor [channel], 
                 x_pos);
      else
        fprintf (output_handle, "%f %f lineto\n", 
                 y0 + (gbl$raw_data.data [sample_nr] - gbl$raw_data.average) * (double) gbl$control_data.hourly_y_factor [channel], 
                 x_pos);

    if (gbl$control_data.verbose)
      printf ("                    Channel %d plotted.\n", channel);
  }

  /* Close output file */

  fprintf (output_handle, "stroke\nshowpage\n");
  fclose (output_handle);

  free (gbl$raw_data.data);
  free (gbl$spectrum_data.data);
  free (gbl$filtered_data.data);
}

/*
** filter is the main routine which is called whenever an interactive filtering and display of data or a batch filter shall 
** be applyed to a number of input files.
*/

void filter (unsigned int action)
{
  unsigned int file_nr,
    sample_nr,
    i,
    j;

  double sum,
    value;

  char date_string [STRING_LENGTH],
    line [STRING_LENGTH],
    *row;

  FILE *input_handle;

#ifdef DEBUG
  printf ("filter (action = %d)\n", action);
#endif

  if (action != STK$ACTION_INTERACTIVE && action != STK$ACTION_FILTER)
  {
    printf ("filter> Routine has been called with an invalid action code!\n");
    return;
  }

  if (gbl$control_data.low_pass_frequency != 0)
    stk_create_plan (gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file * 
                     gbl$control_data.number_of_input_files);

  *date_string = (char) 0;

  /* Allocate memory for the fourier transform */

  gbl$control_data.number_of_samples = gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file * 
                                       gbl$control_data.number_of_input_files;
  if (!allocate_memory (gbl$control_data.number_of_input_files))
  {
    printf ("filter> Unable to allocate memory! Abort!\n");
    return;
  }

  gbl$raw_data.max_value = sum = 0.;
  sample_nr = 0;

  for (file_nr = 0; file_nr < gbl$control_data.number_of_input_files; file_nr++)
  {
    if (gbl$control_data.verbose)
      printf ("filter> Process input file >>%s<<.\n", gbl$control_data.input_file_name [file_nr]);

    if (!(input_handle = fopen (gbl$control_data.input_file_name [file_nr], "r")))
    {
      printf ("filter> Unable to open input file >>%s<<!\n", gbl$control_data.input_file_name [file_nr]);

      free (gbl$raw_data.data);
      free (gbl$spectrum_data.data);
      free (gbl$filtered_data.data);

      return;
    }

    if (gbl$control_data.first_line_contains_date) /* Skip first line and remember the start date */
      fgets (date_string, STRING_LENGTH, input_handle);

    for (i = 0; i < gbl$control_data.samples_per_second * gbl$control_data.seconds_per_file; i++)
    {
      if (feof (input_handle))
      {
        printf ("filter> Unexpected end of file - abort!\n");
        fclose (input_handle);

        free (gbl$raw_data.data);
        free (gbl$spectrum_data.data);
        free (gbl$filtered_data.data);

        return;
      }

      fgets (line, STRING_LENGTH, input_handle);

      row = strtok (line, " ");                             /* Get first token from line */
      for (j = 0; row && j < gbl$control_data.channel; j++) /* Skip channels if desired */
        row = strtok (0, " ");

      if (!row)
      {
        printf ("filter> Unexpected end of line - channel number too large? Abort!\n");
        fclose (input_handle);

        free (gbl$raw_data.data);
        free (gbl$spectrum_data.data);
        free (gbl$filtered_data.data);

        return;
      }

      value = atof (row);
      sum += value;
      if (fabs (gbl$raw_data.data [sample_nr++] = value) > gbl$raw_data.max_value) /* Determine max */
        gbl$raw_data.max_value = fabs (value);
    }
    fclose (input_handle);
  }
  gbl$raw_data.average = sum / (double) sample_nr;

  if (gbl$control_data.verbose)
    printf ("                   Maximum value: %f\n\
                   Average value: %f\n", gbl$raw_data.max_value, gbl$raw_data.average);

  filter_and_spectrum_calc ();

  if (action == STK$ACTION_INTERACTIVE) /* Open window only if running in interactive mode! */
  {
    /* Everything is done now, so create window for graphics display */

#ifdef VMS
    tkInitPosition (gbl$control_data.window_x_position, gbl$control_data.window_y_position, gbl$control_data.window_x_size, 
                    gbl$control_data.window_y_size);
    tkInitDisplayMode (TK_SINGLE | TK_DIRECT | TK_RGB);
    if (tkInitWindow (PROGRAM_NAME) == GL_FALSE)
      tkQuit ();
#else
    glutInitWindowPosition (gbl$control_data.window_x_position, gbl$control_data.window_y_position);
    glutInitWindowSize (gbl$control_data.window_x_size, gbl$control_data.window_y_size);

    glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);

    glutCreateWindow (PROGRAM_NAME);
#endif

    gbl$control_data.graphics_object = make_graphics_object ();

#ifdef VMS
    tkDisplayFunc (show_graphics);
    tkExposeFunc (reshape);
    tkReshapeFunc (reshape);
    tkKeyDownFunc (key_event);

    tkExec ();
#else
    glutDisplayFunc (show_graphics);
    glutReshapeFunc (reshape);
    glutKeyboardFunc (key_event);

    glutMainLoop ();
#endif
  }
  else
    write_filtered_data (); /* Running in -filter mode, so just write the resulting data set back */

  free (gbl$raw_data.data);
  free (gbl$spectrum_data.data);
  free (gbl$filtered_data.data);
}

/*
** Main program
*/

int main (int argc, char **argv)
{
  unsigned int status, 
    i;

  char scratch [STRING_LENGTH],
    parameter_name [STRING_LENGTH],
    mode [STRING_LENGTH];

  /* The actual code starts here */

  if (argc < 4) /* To few arguments supplied! */
  {
    printf ("Main> Usage:\n\
             stk -interactive <f_lowpass> <chnl> <out_file> <in_file> .. <in_file>\n\
                 -filter      <f_lowpass> <chnl> <out_file> <in_file> .. <in_file>\n\
                 -daily       <f_lowpass> <chnl> <out_file> <24 times in_file>\n\
                 -hourly                         <out_file> <in_file>\n");
    return -1;
  }

  /* Determine mode */

  if (!strcmp (*++argv, "-interactive"))
  {
    if (argc < 6)
    {
      printf ("Main> Interactive mode selected, but too few files supplied!\n");
      return -1;
    }

    gbl$control_data.action = STK$ACTION_INTERACTIVE;
    gbl$control_data.number_of_input_files = argc - 5; /* Remember number of input files */
  }
  else if (!strcmp (*argv, "-filter"))
  {
    if (argc < 6)
    {
      printf ("Main> Filter mode selected, but too few files supplied!\n");
      return -1;
    }

    gbl$control_data.action = STK$ACTION_FILTER;
    gbl$control_data.number_of_input_files = argc - 5; /* Remember number of input files */
  }
  else if (!strcmp (*argv, "-daily"))
  {
    if (argc != 29)
    {
      printf ("Main> Daily mode selected, but wrong number of files supplied!\n");
      return -1;
    }

    gbl$control_data.action = STK$ACTION_DAILY;
    gbl$control_data.number_of_input_files = argc - 5; /* Remember number of input files */
  }
  else if (!strcmp (*argv, "-hourly"))
  {
    if (argc != 4)
    {
      printf ("Main> Hourly mode selected, but too many files supplied!\n");
      return -1;
    }

    gbl$control_data.action = STK$ACTION_HOURLY;
    gbl$control_data.number_of_input_files = argc - 3; /* Remember number of input files */
  }
  else
  {
    printf ("Main> Illegal action type!\n\
        Valid actions are: -interactive, -filter, -daily and -hourly.\n");
    return -1;
  }

  strcpy (mode, (*argv) + 1);                                  /* Remember selected mode for later output */

  if (gbl$control_data.action != STK$ACTION_HOURLY)
  {
    gbl$control_data.low_pass_frequency = atof (*++argv);      /* Get low pass frequency */
    gbl$control_data.channel = atoi (*++argv);                   /* Get number of channel to be processed */
  }

  strcpy (gbl$control_data.output_file_name, *++argv);         /* Get output file name */

  for (i = 0; i < gbl$control_data.number_of_input_files; i++) /* Remember all input files */
  {
    strcpy (gbl$control_data.input_file_name [i], *++argv);
    if (!strchr (*argv, '.')) /* No extension found, so fake one */
      strcat (gbl$control_data.input_file_name [i], ".raw");
  }

  /* Process the initialisation file */

  if ((status = read_ini_file (INI_FILE_NAME)) != INI$SUCCESS)
  {
    printf ("Main> Unable to open Inifile >>%s<<!\n", INI_FILE_NAME);
    return -1;
  }

  /* Read entries from the [control]-section */

  get_ini_entry ("control", "verbose", "N", scratch);
  gbl$control_data.verbose = (scratch [0] == 'Y' || scratch [0] == 'y' || scratch [0] == 'J' ||
                              scratch [0] == 'j' || scratch [0] == '1');

  /* Read entries from the [file]-section */

  get_ini_entry ("file", "first_line_contains_date", "N", scratch);
  gbl$control_data.first_line_contains_date = (scratch [0] == 'Y' || scratch [0] == 'y' || scratch [0] == 'J' ||
                                               scratch [0] == 'j' || scratch [0] == '1');

  get_ini_entry ("file", "samples_per_second", "0", scratch);
  if ((gbl$control_data.samples_per_second = atoi (scratch)) < 1)
  {
    printf ("Main> [file]samples_per_second contains an illegal value or is missing!\n");
    return -1;
  }

  get_ini_entry ("file", "seconds_per_file", "0", scratch);
  if ((gbl$control_data.seconds_per_file = atoi (scratch)) < 1)
  {
    printf ("Main> [file]seconds_per_file contains an illegal value or is missing!\n");
    return -1;
  }

  /* Read entries from the [fftw]-section */

  get_ini_entry ("fftw", "tune", "N", scratch);
  if (scratch [0] == 'Y' || scratch [0] == 'y' || scratch [0] == 'J' || scratch [0] == 'j' || scratch [0] == '1')
    gbl$control_data.tune = FFTW_MEASURE;
  else
    gbl$control_data.tune = FFTW_ESTIMATE;

  get_ini_entry ("fftw", "wisdom_file_name", "", gbl$control_data.wisdom_file_name);
  gbl$control_data.use_wisdom = (strlen (gbl$control_data.wisdom_file_name) > 0) && gbl$control_data.tune;

  /* Read entries from the [daily_plot]-section */

  get_ini_entry ("daily_plot", "skip_factor", "1", scratch);
  gbl$control_data.daily_skip_factor = atoi (scratch);

  get_ini_entry ("daily_plot", "y_factor", "5.", scratch);
  gbl$control_data.daily_y_factor = atof (scratch);

  get_ini_entry ("daily_plot", "title", "Daily plot", gbl$control_data.daily_title);

  /* Read entries from the [hourly_plot]-section */

  get_ini_entry ("hourly_plot", "number_of_channels", "1", scratch);
  if ((gbl$control_data.number_of_channels = atoi (scratch)) > MAX_HOURLY_CHANNELS)
  {
    printf ("Main> Too many channels for hourly plot specified in INI file!\n");
    return -1;
  }

  for (i = 0; i < MAX_HOURLY_CHANNELS; i++)
  {
    sprintf (parameter_name, "y_factor_%d", i);
    get_ini_entry ("hourly_plot", parameter_name, "1.", scratch);
    gbl$control_data.hourly_y_factor [i] = atof (scratch);

    sprintf (parameter_name, "low_pass_frequency_%d", i);
    get_ini_entry ("hourly_plot", parameter_name, "1.", scratch);
    gbl$control_data.hourly_low_pass_frequency [i] = atof (scratch);
  }

  get_ini_entry ("hourly_plot", "title", "Hourly plot", gbl$control_data.hourly_title);

  get_ini_entry ("hourly_plot", "channel_description", "", gbl$control_data.hourly_channel_description);

  get_ini_entry ("hourly_plot", "skip_factor", "16", scratch);
  gbl$control_data.hourly_skip_factor = atoi (scratch);

  /* Read entries from the [display]-section */

  get_ini_entry ("display", "window_x_position", "0", scratch);
  gbl$control_data.window_x_position = (unsigned int) atoi (scratch);

  get_ini_entry ("display", "window_y_position", "0", scratch);
  gbl$control_data.window_y_position = (unsigned int) atoi (scratch);

  get_ini_entry ("display", "window_x_size", "1000", scratch);
  gbl$control_data.window_x_size = (unsigned int) atoi (scratch);

  get_ini_entry ("display", "window_y_size", "300", scratch);
  gbl$control_data.window_y_size = (unsigned int) atoi (scratch);

  get_ini_entry ("display", "raw_data_offset", "0", scratch);
  gbl$control_data.raw_data_offset = atoi (scratch);

  get_ini_entry ("display", "spectrum_data_offset", "-50", scratch);
  gbl$control_data.spectrum_data_offset = atoi (scratch);

  get_ini_entry ("display", "filtered_data_offset", "-100", scratch);
  gbl$control_data.filtered_data_offset = atoi (scratch);

  get_ini_entry ("display", "curve_height", "10", scratch);
  gbl$control_data.curve_height = (unsigned int) atoi (scratch);

  get_ini_entry ("display", "max_samples_per_screen", "10000", scratch);
  gbl$control_data.max_samples_per_screen = (unsigned int) atoi (scratch);

  get_ini_entry ("display", "display_coordinate_systems", "Y", scratch);
  gbl$control_data.display_coordinate_systems = (scratch [0] == 'Y' || scratch [0] == 'y' || scratch [0] == 'J' ||
                                                 scratch [0] == 'j' || scratch [0] == '1');

  get_ini_entry ("display", "scale_filtered_like_raw", "Y", scratch);
  gbl$control_data.scale_filtered_like_raw = (scratch [0] == 'Y' || scratch [0] == 'y' || scratch [0] == 'J' ||
                                              scratch [0] == 'j' || scratch [0] == '1');

  get_ini_entry ("display", "initial_x_translate", "2", scratch);
  gbl$control_data.x_translate = atof (scratch);;

  get_ini_entry ("display", "initial_y_translate", "-1", scratch);
  gbl$control_data.y_translate = atof (scratch);

  get_ini_entry ("display", "initial_z_translate", "10", scratch);
  gbl$control_data.z_translate = atof (scratch);

  get_ini_entry ("display", "x_translate_increment", "1.", scratch);
  gbl$control_data.x_translate_increment = atof (scratch);

  get_ini_entry ("display", "y_translate_increment", "1.", scratch);
  gbl$control_data.y_translate_increment = atof (scratch);

  get_ini_entry ("display", "scale_factor_multiplier", "2.", scratch);
  gbl$control_data.scale_factor_multiplier = atof (scratch);

  get_ini_entry ("display", "initial_scale_factor", ".003906", scratch);
  gbl$control_data.initial_scale_factor = gbl$control_data.scale_factor = atof (scratch);

  get_ini_entry ("display", "frequency_increment", ".01", scratch);
  gbl$control_data.frequency_increment = fabs (atof (scratch));

  if (gbl$control_data.verbose)
  {
    /* Dump all parameters read from the initialisation file and the command line */

    printf ("Main> Parameters read from the INI file >>%s<<:\n\n\
        [control]\n\
                verbose:                    %c\n\
        [file]\n\
                first_line_contains_date:   %c\n\
                samples_per_second:         %d\n\
                seconds_per_file:           %d\n\
        [fftw]\n\
                wisdom_file_name:           >>%s<< (use wisdom: %c)\n\
                tune:                       %c\n\
        [daily_plot]\n\
                daily_skip_factor:          %d\n\
                y_factor:                   %f\n\
                title:                      >>%s<<\n\
        [hourly_plot]\n\
                number_of_channels:         %d\n\
                y_factor_0:                 %f\n\
                y_factor_1:                 %f\n\
                y_factor_2:                 %f\n\
                y_factor_3:                 %f\n\
                y_factor_4:                 %f\n\
                y_factor_5:                 %f\n\
                y_factor_6:                 %f\n\
                y_factor_7:                 %f\n\
                low_pass_0:                 %f\n\
                low_pass_1:                 %f\n\
                low_pass_2:                 %f\n\
                low_pass_3:                 %f\n\
                low_pass_4:                 %f\n\
                low_pass_5:                 %f\n\
                low_pass_6:                 %f\n\
                low_pass_7:                 %f\n\
                title:                      >>%s<<\n\
                channel_description:        >>%s<<\n\
                skip_factor:                %d\n\
        [display]\n\
                window_x_position:          %d\n\
                window_y_position:          %d\n\
                window_x_size:              %d\n\
                window_y_size:              %d\n\
                raw_data_offset:            %d\n\
                filtered_data_offset:       %d\n\
                spectrum_data_offset:       %d\n\
                curve_height:               %d\n\
                max_samples_per_screen:     %d\n\
                display_coordinate_systems: %c\n\
                scale_filtered_like_raw:    %c\n\
                initial_x_translated:       %f\n\
                initial_y_translated:       %f\n\
                initial_z_translated:       %f\n\
                x_translate_increment:      %f\n\
                y_translate_increment:      %f\n\
                initial_scale_factor:       %f\n\
                scale_factor_multiplier:    %f\n\
                frequency_increment:        %f\n\
\n\
      Parameters read from the command line:\n\n\
                Mode:                       %s\n\
                Low pass frequency:         %f (apply low pass filter: %c)\n\
                Channel number:             %d\n\
                Output file name:           >>%s<<\n\
                Number of input files:      %d\n\
\n",
          INI_FILE_NAME,

          gbl$control_data.verbose ? 'Y' : 'N',

          gbl$control_data.first_line_contains_date ? 'Y' : 'N',
          gbl$control_data.samples_per_second,
          gbl$control_data.seconds_per_file,

          gbl$control_data.wisdom_file_name, gbl$control_data.use_wisdom ? 'Y' : 'N',
          gbl$control_data.tune == FFTW_MEASURE ? 'Y' : 'N',

          gbl$control_data.daily_skip_factor,
          gbl$control_data.daily_y_factor,
          gbl$control_data.daily_title,

          gbl$control_data.number_of_channels,
          gbl$control_data.hourly_y_factor [0], 
          gbl$control_data.hourly_y_factor [1], 
          gbl$control_data.hourly_y_factor [2], 
          gbl$control_data.hourly_y_factor [3], 
          gbl$control_data.hourly_y_factor [4], 
          gbl$control_data.hourly_y_factor [5], 
          gbl$control_data.hourly_y_factor [6], 
          gbl$control_data.hourly_y_factor [7], 
          gbl$control_data.hourly_low_pass_frequency [0],
          gbl$control_data.hourly_low_pass_frequency [1],
          gbl$control_data.hourly_low_pass_frequency [2],
          gbl$control_data.hourly_low_pass_frequency [3],
          gbl$control_data.hourly_low_pass_frequency [4],
          gbl$control_data.hourly_low_pass_frequency [5],
          gbl$control_data.hourly_low_pass_frequency [6],
          gbl$control_data.hourly_low_pass_frequency [7],
          gbl$control_data.hourly_title,
          gbl$control_data.hourly_channel_description,
          gbl$control_data.hourly_skip_factor,

          gbl$control_data.window_x_position,
          gbl$control_data.window_y_position,
          gbl$control_data.window_x_size,
          gbl$control_data.window_y_size,
          gbl$control_data.raw_data_offset,
          gbl$control_data.filtered_data_offset,
          gbl$control_data.spectrum_data_offset,
          gbl$control_data.curve_height,
          gbl$control_data.max_samples_per_screen,
          gbl$control_data.display_coordinate_systems ? 'Y' : 'N',
          gbl$control_data.scale_filtered_like_raw ? 'Y' : 'N',
          gbl$control_data.x_translate,
          gbl$control_data.y_translate,
          gbl$control_data.z_translate,
          gbl$control_data.x_translate_increment,
          gbl$control_data.y_translate_increment,
          gbl$control_data.initial_scale_factor,
          gbl$control_data.scale_factor_multiplier,
          gbl$control_data.frequency_increment,

          mode,
          gbl$control_data.low_pass_frequency, gbl$control_data.low_pass_frequency == 0 ? 'N' : 'Y',
          gbl$control_data.channel,
          gbl$control_data.output_file_name,
          gbl$control_data.number_of_input_files);
  }

  switch (gbl$control_data.action)
  {
    case STK$ACTION_INTERACTIVE:
      filter (STK$ACTION_INTERACTIVE);
      break;
    case STK$ACTION_FILTER:
      filter (STK$ACTION_FILTER);
      break;
    case STK$ACTION_DAILY:
      create_daily_plot ();
      break;
    case STK$ACTION_HOURLY:
      create_hourly_plot ();
      break;
    default:
      printf ("Main> This should never happen! :-) Illegal action type!\n");
      return -1;
  }

  return 0;
}
