/****************************************************************************
                  INTERNATIONAL AVS CENTER
	(This disclaimer must remain at the top of all files)

WARRANTY DISCLAIMER

This module and the files associated with it are distributed free of charge.
It is placed in the public domain and permission is granted for anyone to use,
duplicate, modify, and redistribute it unless otherwise noted.  Some modules
may be copyrighted.  You agree to abide by the conditions also included in
the AVS Licensing Agreement, version 1.0, located in the main module
directory located at the International AVS Center ftp site and to include
the AVS Licensing Agreement when you distribute any files downloaded from 
that site.

The International AVS Center, MCNC, the AVS Consortium and the individual
submitting the module and files associated with said module provide absolutely
NO WARRANTY OF ANY KIND with respect to this software.  The entire risk as to
the quality and performance of this software is with the user.  IN NO EVENT
WILL The International AVS Center, MCNC, the AVS Consortium and the individual
submitting the module and files associated with said module BE LIABLE TO
ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING,
WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY
SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES.

This AVS module and associated files are public domain software unless
otherwise noted.  Permission is hereby granted to do whatever you like with
it, subject to the conditions that may exist in copyrighted materials. Should
you wish to make a contribution toward the improvement, modification, or
general performance of this module, please send us your comments:  why you
liked or disliked it, how you use it, and most important, how it helps your
work. We will receive your comments at avs@ncsc.org.

Please send AVS module bug reports to avs@ncsc.org.

******************************************************************************/
/* --------------------------------------------------------------------------- */
/* /usr/avs/examples/read_image.c was used as an example of writing modules    */
/* shaker v5    Jan 1992     */
/*       Janet L. Jensen     */
/*       James O. Jensen     */ 
/*       Ed Bender           */
/* ------------------------------------------------------- */
#include <stdio.h>
/* IAC CODE CHANGE : #include <math.h> */
#include <avs/avs_math.h>
#include <avs/avs.h>
#include <avs/field.h>
#include <avs/flow.h>
#include "/usr/avs/include/geom.h"
#include "shaker.h"

/* Used for AVSmessage() */
static char file_version[] = "@(#)shaker.c, Ed says don't call him";

AVSinit_modules()
{
    int shaker();

    AVSmodule_from_desc(shaker);
}

/* -----------------------------------------*/
/* This is the description function         */
/* -----------------------------------------*/

static shaker()
{
    int param, param2, param3;
    int shaker_compute();
    
    AVSset_module_name("shaker", MODULE_DATA);

    AVScreate_output_port("output geometry", "geom");

    AVSadd_float_parameter("move", 0.00, -1.0, 1.0);

    param = AVSadd_parameter("Read Ascii File", "string",0,0 ,".vib");
    AVSconnect_widget(param, "browser");

    param2 = AVSadd_parameter("MODE", "integer",1,1 ,MAX_FREQS);
    param3 = AVSadd_float_parameter("FREQUENCY", 1,1,5000 );
    AVSconnect_widget(param3, "typein_real");


    AVSset_compute_proc(shaker_compute);
}

/* ---------------------------- */
/* This is the compute function */
/* ---------------------------- */

static
shaker_compute(output, move, filename, mode)
    GEOMedit_list *output;
    float  *move;
    char *filename;
    int mode;
{
   int i;

   static int natm;
   int atomic_numb[MAX_ATOMS];
   float atom_center_gs[MAX_ATOMS][COORD];
   float atom_center[MAX_ATOMS][COORD];   /* ground state geom from gaussian90 */
   float atom_radius[MAX_ATOMS];
   float atom_color[MAX_ATOMS][COORD];    /* COORD holds r g b */

   static int   nfreq;
   float freq[MAX_FREQS];
   float dx[MAX_FREQS][MAX_ATOMS],  /* vector changes from gaussian90 */
         dy[MAX_FREQS][MAX_ATOMS], 
         dz[MAX_FREQS][MAX_ATOMS];

/* these variables are for makebonds and get_bonds */
   static int nbonds;                   /* number of bonds */
   static int vertices;                 /* number of vertices, 2 per bond */
   float   *bonds[MAX_BONDS][COORD];  /* vertex pairs  assigned in make_bonds */
   float   tbonds[MAX_BONDS][COORD];  /* temporary holder for passing */
   
/* these variables are for avs GEOM calls */
   float ext[6];
   char *objname;
   GEOMobj *atom_obj, *bond_obj;

/* this is for getting the ascii text file that jim creates from the gaussian .chk file */
/* with his fortran read program.                                                       */
   FILE *f1, *fopen();

/* this is for AVSmodify_parameter for collecting the modes */
   int flags;

   /*------------------------------------------------ */
   /* Read in the ascii data file.                    */
   /* It'm important that the ground state is re-read */
   /* each time because it gets clobber in get_move   */
   /* This is necessary, because the bond connectivity*/
   /* is pointed to by bonds. Then *bonds is put into */
   /* tbond[][] for passing to GEOMadd_disjoint_line. */
   /* This way the connectivity should not change.    */
   /* ----------------------------------------------- */

    /* No data file name */
    if (!filename)
      return(1);

    if(AVSparameter_changed("Read Ascii File")){
       /* Attempt to open file */
       if (!(f1 = fopen(filename, "r"))) {
           AVSmessage(file_version, AVS_Warning, AVSmodule, "shaker_compute",
                   "Ok", "Can't open data file %s", filename);
           return(0);
       }

       /*the file is closed in mover */

        mover(filename, f1, &natm, atomic_numb, atom_center_gs, &nfreq, freq, dx,dy,dz);
 
        flags = (AVS_VALUE | AVS_MAXVAL) | (AVS_MINVAL);
        AVSmodify_parameter("MODE",flags, 30, 1, nfreq);
        AVSmodify_float_parameter("FREQUENCY",flags, freq[nfreq-1], freq[0], freq[nfreq-1]);

       /* ---------------------------------------------------------------*/
       /* the bonds connectivity is based on the ground state geometry.*/
       /* bonds points to the two  ends of one bond in make_bonds().     */
       /* Then the actual                                                */  
       /* value of the vertices are then assigned to tbonds[][]          */
       /* before tbonds is passed to GEOMadd_disjoint_line.              */
       /* ---------------------------------------------------------------*/

        make_bonds(natm, atomic_numb, atom_center_gs, atom_center, &nbonds, bonds);

       /* set the radius to some constant arbitrary numb so balls don't change size */
       for(i=0; i<natm; i++){
         atom_radius[i] = .3;
       }
       /* make color sets the color of the atoms */
       make_atom_color(natm, atomic_numb, atom_color);
   }

        flags = AVS_VALUE ;
        AVSmodify_float_parameter("FREQUENCY",flags, freq[mode-1], freq[0], freq[nfreq-1]);

   /* calculate the number of bond vertices, 2 per bond */
   vertices = 2 * nbonds;

   get_move(mode, *move, natm, atom_center_gs, atom_center, dx, dy, dz);

   /* put the value of what bonds is pointing to an put into array that */
   /* GEOMadd_disjoint_line will like                                   */
   for(i=0; i<vertices; i++){
      tbonds[i][0] = *bonds[i][0];
      tbonds[i][1] = *bonds[i][1];
      tbonds[i][2] = *bonds[i][2];
   }

/* -------------------------------- */
/* avs GEOM calls                   */
/* -------------------------------- */


  *output = GEOMinit_edit_list(*output);

   atom_obj=GEOMcreate_sphere(GEOM_NULL,atom_center,atom_radius,NULL,NULL,natm,GEOM_COPY_DATA);
   GEOMset_extent(atom_obj);
   GEOMadd_float_colors(atom_obj,atom_color,natm,GEOM_COPY_DATA);

   objname = "spheres";
   GEOMedit_geometry(*output,objname,atom_obj);

   bond_obj = GEOMcreate_obj(GEOM_POLYTRI,NULL);
   GEOMadd_disjoint_line(bond_obj,tbonds,NULL,vertices,0);


   GEOMedit_geometry(*output,objname,atom_obj);
   GEOMedit_geometry(*output,objname,bond_obj);

   /* tells avs kernel it was successful */ 
    return(1);
}

/*---------------------------------------------------------------------------------- */
/*---------------------------------------------------------------------------------- */

make_atom_color(natm, atomic_numb, atom_color)
   int natm;
   int atomic_numb[MAX_ATOMS];
   float atom_color[MAX_ATOMS][COORD]; /* COORD holds r g b not xyz here */
{
   int i;

/* assign a color for a given atom */
   for(i=0; i<natm; i++){
      switch(atomic_numb[i]){
         case 1:  /* hydrogen is white */
            atom_color[i][0] = 1.0;
            atom_color[i][1] = 1.0;
            atom_color[i][2] = 1.0;
            break;
         case 6:  /* carbon is green */
            atom_color[i][0] = 0.0;
            atom_color[i][1] = 1.0;
            atom_color[i][2] = 0.0;
            break;

         case 8:  /* oxygen is red */
            atom_color[i][0] = 1.0;
            atom_color[i][1] = 0.0;
            atom_color[i][2] = 0.0;
            break;

         case 7:  /* nitrogen is blue */
            atom_color[i][0] = 0.0;
            atom_color[i][1] = 0.0;
            atom_color[i][2] = 1.0;
            break;
      }
   }
}

/*---------------------------------------------------------------------------------- */
/* mover() reads an ascii data file.  see enclosed data file for format.             */
/*---------------------------------------------------------------------------------- */

mover(filename,f1 ,p_natm, atomic_numb, atom_center_gs ,p_nfreq, freq, dx,dy,dz)
   char *filename;
   FILE *f1;                            /* file descriptor) */
   int *p_natm;                         /* pointer to  number of atoms natm in compute function */
   int atomic_numb[MAX_ATOMS];          /* ....atomic number of each atom */
   float atom_center_gs[MAX_ATOMS][COORD]; /* ground state x,y,z coord of each atom */
   int *p_nfreq;                         /* pointer to number of freqs (3n) in compute function */
   float freq[MAX_FREQS];               /* ....the value of each freq */
   float dx[MAX_FREQS][MAX_ATOMS],      /* ....displacement vectors */
         dy[MAX_FREQS][MAX_ATOMS],
         dz[MAX_FREQS][MAX_ATOMS];
{

   int i, j;

  /* ------------------------------------ */
  /* read number of atoms, atomic numbers */
  /* ------------------------------------ */

   fscanf(f1, "%d", p_natm);  /* the first number is the number of atoms */

   for(i=0; i<*p_natm; i++){   /* the next set of numbers are the atomic numbers */
      fscanf(f1, "%d", &atomic_numb[i]);
   }

  /* -------------------------------- */
  /* read the ground state geometries */
  /* -------------------------------- */


   for(i=0; i<*p_natm; i++){
      fscanf(f1, "%f %f %f", &atom_center_gs[i][0], /* x */
                             &atom_center_gs[i][1], /* y */
                             &atom_center_gs[i][2]); /* z */
   }

   /* -------------------------------- */
   /* read the frequency data          */
   /* -------------------------------- */

   /* the number of frequencies is 3 * number of atoms */
   /* the first six are usually not very interesting   */
   /* but jim says to put them in anyway               */

   *p_nfreq = 3 **p_natm;
   for(i=0; i<*p_nfreq; i++){
      fscanf(f1, "%f", &freq[i]);   /* the value of the frequency */
      for(j=0; j<*p_natm; j++){       /* read displacement vectors  */
         fscanf(f1, "%f %f %f", &dx[i][j], &dy[i][j], &dz[i][j]);
      }
   }
   fclose(f1);  /* close the file before returning to shaker_compute */
                /* file is only read in mover                        */
}

/*---------------------------------------------------------------------------------- */
/* make_bonds returns the number of vertices (2 per bond) and array bonds that       */
/*            have the connectivity                                                  */
/*---------------------------------------------------------------------------------- */
make_bonds(natm, atomic_numb, atom_center_gs, atom_center, p_nbonds, bond_list)
   int natm;
   int atomic_numb[MAX_ATOMS];
   float atom_center_gs[MAX_ATOMS][COORD];
   float atom_center[MAX_ATOMS][COORD];
   int *p_nbonds;
   float *bond_list[MAX_BONDS][COORD];
{
   int i, j;
   int bond_index;
   float distance, distance_squared;
   float bohr;

  *p_nbonds = 0; /* initialize bond counter */
   bohr = 3.02;   /* maximum bond length possible */
   bond_index =0;  /* for making the bond list indices */

   for(i=0; i<natm; i++){
      for(j=i+1; j<natm; j++){
         /* if both atoms are hydrogen, then skip */
         if(atomic_numb[i] != 1 ||  atomic_numb[j] != 1){
            /* then they aren't both hydrogens, go ahead to check bond */
            distance_squared = pow((atom_center_gs[i][0] - atom_center_gs[j][0]),2)
                              +pow((atom_center_gs[i][1] - atom_center_gs[j][1]),2)
                              +pow((atom_center_gs[i][2] - atom_center_gs[j][2]),2);
            distance = sqrt(distance_squared);

            if(distance < bohr){       /* then count bond and form bond */
               *p_nbonds = *p_nbonds + 1;

               bond_list[bond_index][0] = &atom_center[i][0];
               bond_list[bond_index][1] = &atom_center[i][1];
               bond_list[bond_index][2] = &atom_center[i][2];

               bond_list[bond_index+1][0] = &atom_center[j][0];
               bond_list[bond_index+1][1] = &atom_center[j][1];
               bond_list[bond_index+1][2] = &atom_center[j][2];

                bond_index = bond_index + 2;
            }
         }
      }
   }
}

/* ----------------------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------------------------- */

get_move(mode, fmove, natm, atom_center_gs, atom_center, dx, dy, dz)
   int mode;
   float fmove;
   int natm;
   float atom_center_gs[MAX_ATOMS][COORD];
   float atom_center[MAX_ATOMS][COORD];
   float dx[MAX_FREQS][MAX_ATOMS],
         dy[MAX_FREQS][MAX_ATOMS],
         dz[MAX_FREQS][MAX_ATOMS];
{
   int i;
   int set_freq;

   set_freq = mode-1;
   for(i=0; i<natm; i++){
      atom_center[i][0] = atom_center_gs[i][0] + (fmove * dx[set_freq][i]);
      atom_center[i][1] = atom_center_gs[i][1] + (fmove * dy[set_freq][i]);
      atom_center[i][2] = atom_center_gs[i][2] + (fmove * dz[set_freq][i]);
   }
}