/*****************************************************************************
 *                                                                           *
 *               This script reads in a VIBEIGVEC file produced              *
 *               by CPMD and writes the output in the form that              *
 *                         can be used in Mathematica.                       *
 *                                                                           *
 *                           Veronika Brazdova 2002                          *
 *                           vb@chemie.hu-berlin.de                          *
 *                                                                           *
 *****************************************************************************/

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

#define BUFFER 150

/***********************************main**************************************/

/// Return codes for function main:
//  0 successful run
//  1 no command line arguments or file not found
//  2 not enough memory
//  3 file not in a correct format

int main(int argc, char *argv[])
{
  int    i, j, k, tmp, found, no_atoms, no_modes;
  double at_coords[3];
  double **modes, *freq;
  char   line[BUFFER], at_symbol[3];
  int    at_number;
  FILE   *fp;


  // argv[0] = name-of-the-executable
  // argv[1] = number of modes
  
  /// Handle the command line arguments
  if (argc < 2)
   {
     fprintf(stderr,"usage: vib2math no_normal_modes > file.tab\n");
     return 1;
   }
  
  no_modes = atoi(argv[1]);
  if (!(fp = fopen("VIBEIGVEC","r")))
    {
      fprintf(stderr,"Vib2math: cannot open file VIBEIGVEC\n");
      return 1;
    }

  /// Memory allocation for the array which will contain the
  //  normal mode vectors. Because we don't know how many
  //  columns there are at the end of the VIBEIGVEC file,
  //  we'll allocate a bit more.

  modes = (double **) malloc(sizeof(double *) * (no_modes+8));
  if (modes == NULL) /* not enough memory */
  {
    fprintf(stderr,"Not enough memory.\n");
    return 2;
  }	    

  for (i = 0; i < no_modes; i++ )
  {
    modes[i] = (double *) malloc(sizeof(double) * no_modes);
    if (modes[i] == NULL) 
    {
      fprintf(stderr,"Not enough memory.\n");
      return 2;
    }
  }

  freq = (double *) malloc(sizeof(double) * no_modes);
  if (freq == NULL) 
  {
    fprintf(stderr,"Not enough memory.\n");
    return 2;
  }	    

  /// Read in the frequencies and the normal vectors 
  
  found = 0;
  while (fgets(line,BUFFER,fp))
  {
    if (strstr(line,"NEW SET"))
    {
      for (i = 0; i < (no_modes-1)/8+1; i++)
      {		  
        if (fgets(line,BUFFER,fp) == NULL)
        {
          fprintf(stderr,"Incorrect format of file VIBEIGVEC.\n");
          return 3;
        }
        if (fgets(line,BUFFER,fp) == NULL)
        {
          fprintf(stderr,"Incorrect format of file VIBEIGVEC.\n");
          return 3;
        }
	if (fgets(line,BUFFER,fp) == NULL)
          return 3;
        tmp = i*8;
        for (j = 0; j < no_modes; j++)
        {
          if (fgets(line,BUFFER,fp) == NULL)
            return 3;
          sscanf(line,"%lf %lf %lf %lf %lf %lf %lf %lf",&modes[tmp][j],\
            &modes[tmp+1][j],&modes[tmp+2][j],\
            &modes[tmp+3][j],&modes[tmp+4][j],\
            &modes[tmp+5][j],&modes[tmp+6][j],\
            &modes[tmp+7][j]);
        }      

      }
      found = 1;	
      for (j = 0; j < no_modes; j++)
      {
        printf(" f%i = {\n",j+1);
        for (i = 0; i < no_modes - 1 ; i++)
          printf(" %10.6lf,\n",modes[j][i]);
        printf(" %10.6lf\n}\n",modes[j][no_modes-1]);
      }
    }
  }  
  fclose(fp);

  if (!found)
    {
      fprintf(stderr,\
        "The file VIBEIGVEC does not contain projected frequencies.\n");
      return 3;
    }
}  /// main