trunk/ripple/readruf.c

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

// Structures to store our data:

// coords holds the data for a single tethered dipole:
struct coords{
  double pos[3]; // cartesian coords
  double theta;  // orientational angle relative to z axis
  double phi;    // orientational angle in x-y plane
  double mu;     // dipole strength
  char name[30]; // an identifier for the type of atom
};

// state holds the current "configuration" of the entire system
struct system {
  int    nAtoms;           // Number of Atoms in this configuration
  struct coords *r;        // The set of coordinates for all atoms
  double beta;             // beta = 1 /(kb*T)
  double strength;         // strength of the dipoles (Debye)
  double z0;               // default z axis position
  double theta0;           // default theta angle
  double kz;               // force constant for z displacement
  double ktheta;           // force constant for theta displacement
  int    nCycles;          // How many cycles to do in total
  int    iCycle;           // How many cycles have we done?
  int    nMoves;           // How many MC moves in each cycle
  int    nSample;          // How many cycles between samples
  double Hmat[2][2];       // The information about the size of the per. box
  double HmatI[2][2];      // The inverse box
  double energy;           // The current Energy
  double dtheta;           // maximum size of a theta move
  double deltaz;           // maximum size of a z move
  double deltaphi;         // maximum size of a phi move
  int    nAttempts;        // number of MC moves that have been attempted
  int    nAccepts;         // number of MC moves that have been accepted
  int    nx;               // number of unit cells in x direction
  int    ny;               // number of unit cells in y direction
  struct system *next;      // Next frame in the linked list
};

char *program_name;

int main(argc, argv)
     int argc;
     char *argv[];
{
  FILE *in_file;
  char in_name[500];
  char *eof_test, *foo;
  char read_buffer[1000];
  int lineCount = 0;
  double *mag, *newmag;
  int *present_in_old;
  double uxi, uyi, uzi, p1;
  double aLat, bLat;
  int cells;
  double sumZ, sumUx, sumUy, sumUz, sumP;
  double interpsum, value;
  int ninterp, px, py, newp;
  int i, j, nloops;
  int newx, newy, newindex, index;
  int new_i, new_j, new_index;
  int N, nframes;
  double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
  double maxfreqx, maxfreqy, maxfreq, dfreq;
  double dx, dy, dx1, dy1, pt1x, pt1y, pt2x, pt2y;
  int whichbin;
  int nbins, nx, ny;
  int *samples;
  double *bin;
  double sumz, sumz2, averz, averz2, sigmaz;

  int done;
  char current_flag;
  
  program_name = argv[0];

  for(i = 1; i < argc; i++){
    if(argv[i][0] == '-'){
      done = 0;
      j = 1;
      current_flag = argv[i][j];
      while( (current_flag != '\0') && ( !done ) ){
        switch(current_flag){
          case 'i':
            i++;
            strcpy( in_name, argv[i] );
            done = 1;
            break;
        }
      }
    }
  }

  struct system* state;
  struct system* temp_state;
  struct coords* r;

/*   FFTW_COMPLEX *in, *out; */
/*   fftwnd_plan p; */
        
/*   nbins = 100; */
/*   bin = (double *) malloc(sizeof(double) * nbins); */
/*   samples = (int *) malloc(sizeof(int) * nbins); */
/*   for (i=0; i < nbins; i++) { */
/*     bin[i] = 0.0; */
/*     samples[i] = 0; */
/*   } */
  
  in_file = fopen(in_name, "r");
  if(in_file == NULL){
    printf("Cannot open file \"%s\" for reading.\n", in_name);
    exit(8);
  }

  nframes = 0;
   
  // start reading the first frame

  eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
  lineCount++;

  state = (struct system *) malloc(sizeof(struct system));
  state->next = NULL;
  state->strength = 10.0;

  printf("averz\taverz2\tsigmaz\n");


  while(eof_test != NULL){
    
    nframes++;
    sumz = 0.0;
    sumz2 = 0.0;
    averz = 0.0;
    averz2 = 0.0;
    sigmaz = 0.0;
    (void)sscanf(read_buffer, "%d", &state->nAtoms);
    N = 2 * state->nAtoms;    

    state->r = (struct coords *)calloc(N, sizeof(struct coords));
    mag = (double *) malloc(sizeof(double) * N);
    newmag = (double *) malloc(sizeof(double) * N);
    present_in_old = (int *) malloc(sizeof(int) * N);
        
    // read and the comment line and grab the time and box dimensions
    
    eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
    lineCount++;
    if(eof_test == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    
    foo = strtok( read_buffer, " ,;\t\n" );
    (void)sscanf( read_buffer, "%d", &state->iCycle );

    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%d", &state->nx);

    nx = state->nx;

    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%d", &state->ny);
   
    ny = state->ny;

    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[0][0]);
    
    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[1][0]);
    
    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[0][1]);
    
    foo = strtok(NULL, " ,;\t\0");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&state->Hmat[1][1]);
    
    // Length of the two box vectors:

    dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
    dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));

    aLat = dx / (double)(state->nx);
    bLat = dy / (double)(state->ny);
 

    // FFT stuff depends on nx and ny, so delay allocation until we have
    // that information
    
/*     in = fftw_malloc(sizeof(FFTW_COMPLEX) * N); */
/*     out = fftw_malloc(sizeof(FFTW_COMPLEX) * N); */
/*     p =  fftw2d_create_plan(2*state->nx,  */
/*                          2*state->ny,  */
/*                          FFTW_FORWARD,  */
/*                          FFTW_ESTIMATE); */


/*     for (i = 0; i < N; i++) { */
/*       present_in_old[i] = 0; */
/*     } */

    for (i=0;i<2*state->nx;i=i+2){
      for(j=0;j<2*state->ny;j++){
        newy = j;
        if ((j % 2) == 0) {
          newx = i;
        } else {
          newx = i + 1;
        }
        newindex = newx*2*state->ny + newy;       

        eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
        lineCount++;
        if(eof_test == NULL){
          printf("error in reading file at line: %d\n", lineCount);
          exit(8);
        }
        
        foo = strtok(read_buffer, " ,;\t\0");
        (void)strcpy(state->r[newindex].name, foo); //copy the atom name
        
        // next we grab the positions
        
        foo = strtok(NULL, " ,;\t\0");
        if(foo == NULL){
          printf("error in reading postition x from %s\n"
                 "natoms  = %d, line = %d\n",
                 in_name, state->nAtoms, lineCount );
          exit(8);
        }
        (void)sscanf( foo, "%lf", &state->r[newindex].pos[0] );
      
        foo = strtok(NULL, " ,;\t\0");
        if(foo == NULL){
          printf("error in reading postition y from %s\n"
                 "natoms  = %d, line = %d\n",
                 in_name, state->nAtoms, lineCount );
          exit(8);
        }
        (void)sscanf( foo, "%lf", &state->r[newindex].pos[1] );
        
        foo = strtok(NULL, " ,;\t\0");
        if(foo == NULL){
          printf("error in reading postition z from %s\n"
                 "natoms  = %d, line = %d\n",
                 in_name, state->nAtoms, lineCount );
          exit(8);
        }
        (void)sscanf( foo, "%lf", &state->r[newindex].pos[2] );
        sumz += state->r[newindex].pos[2];
        sumz2 += pow(state->r[newindex].pos[2],2);
        
        foo = strtok(NULL, " ,;\t\0");
        if(foo == NULL){
          printf("error in reading angle phi from %s\n"
                 "natoms  = %d, line = %d\n",
                 in_name, state->nAtoms, lineCount );
          exit(8);
        }
        (void)sscanf( foo, "%lf", &state->r[newindex].phi );
        
        foo = strtok(NULL, " ,;\t\0");
        if(foo == NULL){
          printf("error in reading unit vector x from %s\n"
                 "natoms  = %d, line = %d\n",
                 in_name, state->nAtoms, lineCount );
          exit(8);
        }
        (void)sscanf( foo, "%lf", &uxi );
        
        foo = strtok(NULL, " ,;\t\0");
        if(foo == NULL){
          printf("error in reading unit vector y from %s\n"
                 "natoms  = %d, line = %d\n",
                 in_name, state->nAtoms, lineCount );
          exit(8);
        }
        (void)sscanf( foo, "%lf", &uyi );
        
        foo = strtok(NULL, " ,;\t\0");
        if(foo == NULL){
          printf("error in reading unit vector z from %s\n"
                 "natoms  = %d, line = %d\n",
                 in_name, state->nAtoms, lineCount );
          exit(8);
        }
        (void)sscanf( foo, "%lf", &uzi );
        
        state->r[newindex].theta = acos(uzi);
        
        // The one parameter not stored in the dump file is the dipole strength
        state->r[newindex].mu = state->strength;

        present_in_old[newindex] = 1;
      }
    }
    averz = sumz / state->nAtoms;
    averz2 = sumz2 / state->nAtoms;
    sigmaz = sqrt(averz2 - pow(averz,2));
    printf("%lf\t%lf\t%lf\n", averz, averz2, sigmaz);
    
/*     for (i=0; i< 2*state->nx; i++) { */
/*       for(j=0; j< 2*state->ny; j++) { */
/*      newindex = i*2*state->ny + j; */
/*      mag[newindex] = 0.0; */
/*      newmag[newindex] = 0.0; */
/*       } */
/*     } */
    
/*     for (i=0; i< 2*state->nx; i++) { */
/*       for(j=0; j< 2*state->ny; j++) { */
/*      newindex = i*2*state->ny + j; */
/*      if (present_in_old[newindex] == 0) { */
/*        // interpolate from surrounding points: */
          
/*        interpsum = 0.0; */
/*        ninterp = 0; */
          
/*        //point1 = bottom; */
          
/*        px = i - 1; */
/*        py = j; */
/*        newp = px*2*state->ny + py; */
/*        if (px >= 0) { */
/*          interpsum += state->r[newp].pos[2]; */
/*          ninterp++; */
/*          state->r[newindex].pos[1] = state->r[newp].pos[1]; */
/*        } */
          
/*        //point2 = top; */
          
/*        px = i + 1; */
/*        py = j; */
/*        newp = px*2*state->ny + py; */
/*        if (px < 2*state->nx) { */
/*          interpsum += state->r[newp].pos[2]; */
/*          ninterp++; */
/*          state->r[newindex].pos[1] = state->r[newp].pos[1]; */
/*        }  */
          
/*        //point3 = left; */
          
/*        px = i; */
/*        py = j - 1; */
/*        newp = px*2*state->ny + py; */
/*        if (py >= 0) { */
/*          interpsum += state->r[newp].pos[2]; */
/*          ninterp++; */
/*          state->r[newindex].pos[0] = state->r[newp].pos[0]; */
/*        }  */
          
/*        //point4 = right; */
          
/*        px = i; */
/*        py = j + 1; */
/*        newp = px*2*state->ny + py; */
/*        if (py < 2*state->ny) { */
/*          interpsum += state->r[newp].pos[2]; */
/*          ninterp++; */
/*          state->r[newindex].pos[0] = state->r[newp].pos[0]; */
/*        }  */
          
/*        value = interpsum / (double)ninterp; */
          
/*        state->r[newindex].pos[2] = value; */
/*      } */
/*       } */
/*     } */

    
/*     for (i=0; i < 2*state->nx; i++) { */
/*       for (j=0; j < 2*state->ny; j++) { */
/*      newindex = i*2*state->ny + j; */

/*      c_re(in[newindex]) = state->r[newindex].pos[2]; */
/*      c_im(in[newindex]) = 0.0; */
/*       }  */
/*     } */

    
/*     fftwnd(p, 1, in, 1, 0, out, 1, 0);  */
    
/*     for (i=0; i< 2*state->nx; i++) { */
/*       for(j=0; j< 2*state->ny; j++) { */
/*         newindex = i*2*state->ny + j; */
/*         mag[newindex] += sqrt(pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2)); */
/*       } */
/*     } */
    
            
    temp_state = state->next;
    state->next = NULL;
    
    if (temp_state != NULL) {
      free(temp_state->r);
      free(temp_state);
    }
    

/*     fftwnd_destroy_plan(p); */
/*     fftw_free(out); */
/*     fftw_free(in); */

    free(present_in_old);

    // Make a new frame
    
    temp_state = (struct system *) malloc(sizeof(struct system));
    temp_state->next = state;
    state = temp_state;
    eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
    lineCount++;
    
  }

  
  // This stuff is for printing out the FFT results directly for matlab:
  
/*   for (i=0; i < 2*nx; i++) { */
/*     for(j=0; j< 2*ny; j++) { */
/*       newindex = i*2*ny + j; */
/*       printf("%lf\t",  mag[newindex]/(double)nframes); */
/*     } */
/*     printf("\n"); */
/*   } */
  
  
  // This stuff is for stitching the four quadrants together: 
  
/*   for (i=0; i< (2*nx/2); i++) { */
/*     for(j=0; j< (2*ny/2); j++) { */
/*       index = i*2*ny + j; */
/*       new_i = i + (2*nx/2); */
/*       new_j = j + (2*ny/2); */
/*       new_index = new_i*2*ny + new_j; */
/*       newmag[new_index] = mag[index]; */
/*     } */
/*   } */
  
/*   for (i=(2*nx/2); i< 2*nx; i++) { */
/*     for(j=0; j< (2*ny/2); j++) { */
/*       index = i*2*ny + j; */
/*       new_i = i - (2*nx/2); */
/*       new_j = j + (2*ny/2); */
/*       new_index = new_i*2*ny + new_j; */
/*       newmag[new_index] = mag[index]; */
/*     } */
/*   } */
  
/*   for (i=0; i< (2*nx/2); i++) { */
/*     for(j=(2*ny/2); j< 2*ny; j++) { */
/*       index = i*2*ny + j; */
/*       new_i = i + (2*nx/2); */
/*       new_j = j - (2*ny/2); */
/*       new_index = new_i*2*ny + new_j; */
/*       newmag[new_index] = mag[index]; */
/*       } */
/*   } */
  
/*   for (i=(2*nx/2); i< 2*nx; i++) { */
/*     for(j=(2*ny/2); j< 2*ny; j++) { */
/*       index = i*2*ny + j; */
/*       new_i = i - (2*nx/2); */
/*       new_j = j - (2*ny/2); */
/*       new_index = new_i*2*ny + new_j; */
/*       newmag[new_index] = mag[index]; */
/*     } */
/*   } */
  
/*   maxfreqx = 2.0 / aLat; */
/*   maxfreqy = 2.0 / bLat; */
  
/*   //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy); */
  
/*   maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy); */
/*   dfreq = maxfreq/(double)(nbins-1); */
  
/*   //printf("%lf\n", dfreq); */
  
/*   zero_freq_x = nx;  */
/*   zero_freq_y = ny;  */
  
/*   for (i=0; i< 2*nx; i++) { */
/*     for(j=0; j< 2*ny; j++) { */
      
/*       freq_x = (double)(i - zero_freq_x)*maxfreqx / nx; */
/*       freq_y = (double)(j - zero_freq_y)*maxfreqy / ny; */
      
/*       freq = sqrt(freq_x*freq_x + freq_y*freq_y); */
      
/*       whichbin = (int) (freq / dfreq); */
/*       newindex = i*2*ny + j; */
/*       //     printf("%d %d %lf %lf\n", whichbin, newindex, freq, dfreq); */
/*       bin[whichbin] += mag[newindex]; */
/*       samples[whichbin]++; */
/*     } */
/*   } */
  
  
/*   for (i = 0; i < nbins; i++) {     */
/*     if (samples[i] > 0) { */
/*       printf("%lf\t%lf\n", i*dfreq, bin[i]/(double)samples[i]); */
/*     } */
/*   } */

    // This stuff is for stitching the four quadrants together: 
    
/*     for (i=0; i< (2*nx/2); i++) { */
/*       for(j=0; j< (2*ny/2); j++) { */
/*      index = i*2*ny + j; */
/*      new_i = i + (2*nx/2); */
/*      new_j = j + (2*ny/2); */
/*      new_index = new_i*2*ny + new_j; */
/*      newmag[new_index] = mag[index]; */
/*       } */
/*     } */
    
/*     for (i=(2*nx/2); i< 2*nx; i++) { */
/*       for(j=0; j< (2*ny/2); j++) { */
/*      index = i*2*ny + j; */
/*      new_i = i - (2*nx/2); */
/*      new_j = j + (2*ny/2); */
/*      new_index = new_i*2*ny + new_j; */
/*      newmag[new_index] = mag[index]; */
/*       } */
/*     } */
    
/*     for (i=0; i< (2*nx/2); i++) { */
/*       for(j=(2*ny/2); j< 2*ny; j++) { */
/*      index = i*2*ny + j; */
/*      new_i = i + (2*nx/2); */
/*      new_j = j - (2*ny/2); */
/*      new_index = new_i*2*ny + new_j; */
/*      newmag[new_index] = mag[index]; */
/*       } */
/*     } */
    
/*     for (i=(2*nx/2); i< 2*nx; i++) { */
/*       for(j=(2*ny/2); j< 2*ny; j++) { */
/*      index = i*2*ny + j; */
/*      new_i = i - (2*nx/2); */
/*      new_j = j - (2*ny/2); */
/*      new_index = new_i*2*ny + new_j; */
/*      newmag[new_index] = mag[index]; */
/*       } */
/*     } */

/*    for (i=0; i < 2*nx; i++) { */
/*      for(j=0; j< 2*ny; j++) { */
/*        newindex = i*2*ny + j; */
/*        printf("%lf\t",newmag[newindex]/(double)nframes); */
/*      } */
/*      printf("\n"); */
/*    } */
  
 // This stuff is for printing out the FFT results directly:
    
/*    for (i=0; i < 2*nx; i++) { */
/*      for(j=0; j< 2*ny; j++) { */
/*        newindex = i*2*ny + j; */
/*        printf("%lf\t", mag[newindex]/(double)nframes); */
/*      } */
/*      printf("\n"); */
/*    } */

  fclose(in_file);
  free(newmag);
  free(mag);
  free(samples);
  free(bin);
  free(state);
  return(0);
}

/*      } */
/*       for(i=0;i<120;i++) printf("%f\t%f\n", c_re(out[i]), c_im(out[i])); */
/*       fftw_destroy_plan(p); */
/*       fftw_free(in); */
/*       fftw_free(out); */
/*       averZ = sumZ / (double)nloops; */
/*       averUx = sumUx / (double)nloops; */
/*       averUy = sumUy / (double)nloops; */
/*       averUz = sumUz / (double)nloops; */
/*       averP = 1.5*sumP / (double)nloops-0.5; */
/*       averTheta = acos(averUz); */
/*       printf("nloops=%d\n",nloops); */
/*       printf("sumZ=%f\n",sumZ); */
/*       printf("average height is : %f\n",averZ); */
/*       printf("average ux is : %f\n",averUx); */
/*       printf("average uy is : %f\n",averUy); */
/*       printf("average uz is : %f\n",averUz); */
/*       printf("average angle is : %f\n",averTheta); */
/*       printf("average p is : %f\n",averP); */


/*       return 0; */
/*       }*/
Revision:	1185
Committed:	Fri May 21 20:07:10 2004 UTC (20 years, 1 month ago) by xsun
Content type:	text/plain
File size:	17242 byte(s)
Log Message:	this is the ripple codes
#	User	Rev	Content
1	xsun	1185	#include<stdio.h>
2			#include<string.h>
3			#include<stdlib.h>
4			#include<math.h>
5			#include<fftw.h>
6
7			// Structures to store our data:
8
9			// coords holds the data for a single tethered dipole:
10			struct coords{
11			double pos[3]; // cartesian coords
12			double theta; // orientational angle relative to z axis
13			double phi; // orientational angle in x-y plane
14			double mu; // dipole strength
15			char name[30]; // an identifier for the type of atom
16			};
17
18			// state holds the current "configuration" of the entire system
19			struct system {
20			int nAtoms; // Number of Atoms in this configuration
21			struct coords *r; // The set of coordinates for all atoms
22			double beta; // beta = 1 /(kb*T)
23			double strength; // strength of the dipoles (Debye)
24			double z0; // default z axis position
25			double theta0; // default theta angle
26			double kz; // force constant for z displacement
27			double ktheta; // force constant for theta displacement
28			int nCycles; // How many cycles to do in total
29			int iCycle; // How many cycles have we done?
30			int nMoves; // How many MC moves in each cycle
31			int nSample; // How many cycles between samples
32			double Hmat[2][2]; // The information about the size of the per. box
33			double HmatI[2][2]; // The inverse box
34			double energy; // The current Energy
35			double dtheta; // maximum size of a theta move
36			double deltaz; // maximum size of a z move
37			double deltaphi; // maximum size of a phi move
38			int nAttempts; // number of MC moves that have been attempted
39			int nAccepts; // number of MC moves that have been accepted
40			int nx; // number of unit cells in x direction
41			int ny; // number of unit cells in y direction
42			struct system *next; // Next frame in the linked list
43			};
44
45			char *program_name;
46
47			int main(argc, argv)
48			int argc;
49			char *argv[];
50			{
51			FILE *in_file;
52			char in_name[500];
53			char eof_test, foo;
54			char read_buffer[1000];
55			int lineCount = 0;
56			double mag, newmag;
57			int *present_in_old;
58			double uxi, uyi, uzi, p1;
59			double aLat, bLat;
60			int cells;
61			double sumZ, sumUx, sumUy, sumUz, sumP;
62			double interpsum, value;
63			int ninterp, px, py, newp;
64			int i, j, nloops;
65			int newx, newy, newindex, index;
66			int new_i, new_j, new_index;
67			int N, nframes;
68			double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
69			double maxfreqx, maxfreqy, maxfreq, dfreq;
70			double dx, dy, dx1, dy1, pt1x, pt1y, pt2x, pt2y;
71			int whichbin;
72			int nbins, nx, ny;
73			int *samples;
74			double *bin;
75			double sumz, sumz2, averz, averz2, sigmaz;
76
77			int done;
78			char current_flag;
79
80			program_name = argv[0];
81
82			for(i = 1; i < argc; i++){
83			if(argv[i][0] == '-'){
84			done = 0;
85			j = 1;
86			current_flag = argv[i][j];
87			while( (current_flag != '\0') && ( !done ) ){
88			switch(current_flag){
89			case 'i':
90			i++;
91			strcpy( in_name, argv[i] );
92			done = 1;
93			break;
94			}
95			}
96			}
97			}
98
99			struct system* state;
100			struct system* temp_state;
101			struct coords* r;
102
103			/* FFTW_COMPLEX in, out; */
104			/* fftwnd_plan p; */
105
106			/* nbins = 100; */
107			/* bin = (double ) malloc(sizeof(double) nbins); */
108			/* samples = (int ) malloc(sizeof(int) nbins); */
109			/* for (i=0; i < nbins; i++) { */
110			/* bin[i] = 0.0; */
111			/* samples[i] = 0; */
112			/* } */
113
114			in_file = fopen(in_name, "r");
115			if(in_file == NULL){
116			printf("Cannot open file \"%s\" for reading.\n", in_name);
117			exit(8);
118			}
119
120			nframes = 0;
121
122			// start reading the first frame
123
124			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
125			lineCount++;
126
127			state = (struct system *) malloc(sizeof(struct system));
128			state->next = NULL;
129			state->strength = 10.0;
130
131			printf("averz\taverz2\tsigmaz\n");
132
133
134			while(eof_test != NULL){
135
136			nframes++;
137			sumz = 0.0;
138			sumz2 = 0.0;
139			averz = 0.0;
140			averz2 = 0.0;
141			sigmaz = 0.0;
142			(void)sscanf(read_buffer, "%d", &state->nAtoms);
143			N = 2 * state->nAtoms;
144
145			state->r = (struct coords *)calloc(N, sizeof(struct coords));
146			mag = (double ) malloc(sizeof(double) N);
147			newmag = (double ) malloc(sizeof(double) N);
148			present_in_old = (int ) malloc(sizeof(int) N);
149
150			// read and the comment line and grab the time and box dimensions
151
152			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
153			lineCount++;
154			if(eof_test == NULL){
155			printf("error in reading file at line: %d\n", lineCount);
156			exit(8);
157			}
158
159			foo = strtok( read_buffer, " ,;\t\n" );
160			(void)sscanf( read_buffer, "%d", &state->iCycle );
161
162			foo = strtok(NULL, " ,;\t\0");
163			if(foo == NULL){
164			printf("error in reading file at line: %d\n", lineCount);
165			exit(8);
166			}
167			(void)sscanf(foo, "%d", &state->nx);
168
169			nx = state->nx;
170
171			foo = strtok(NULL, " ,;\t\0");
172			if(foo == NULL){
173			printf("error in reading file at line: %d\n", lineCount);
174			exit(8);
175			}
176			(void)sscanf(foo, "%d", &state->ny);
177
178			ny = state->ny;
179
180			foo = strtok(NULL, " ,;\t\0");
181			if(foo == NULL){
182			printf("error in reading file at line: %d\n", lineCount);
183			exit(8);
184			}
185			(void)sscanf(foo, "%lf",&state->Hmat[0][0]);
186
187			foo = strtok(NULL, " ,;\t\0");
188			if(foo == NULL){
189			printf("error in reading file at line: %d\n", lineCount);
190			exit(8);
191			}
192			(void)sscanf(foo, "%lf",&state->Hmat[1][0]);
193
194			foo = strtok(NULL, " ,;\t\0");
195			if(foo == NULL){
196			printf("error in reading file at line: %d\n", lineCount);
197			exit(8);
198			}
199			(void)sscanf(foo, "%lf",&state->Hmat[0][1]);
200
201			foo = strtok(NULL, " ,;\t\0");
202			if(foo == NULL){
203			printf("error in reading file at line: %d\n", lineCount);
204			exit(8);
205			}
206			(void)sscanf(foo, "%lf",&state->Hmat[1][1]);
207
208			// Length of the two box vectors:
209
210			dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
211			dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));
212
213			aLat = dx / (double)(state->nx);
214			bLat = dy / (double)(state->ny);
215
216
217			// FFT stuff depends on nx and ny, so delay allocation until we have
218			// that information
219
220			/* in = fftw_malloc(sizeof(FFTW_COMPLEX) * N); */
221			/* out = fftw_malloc(sizeof(FFTW_COMPLEX) * N); */
222			/* p = fftw2d_create_plan(2state->nx, /
223			/* 2state->ny, /
224			/* FFTW_FORWARD, */
225			/* FFTW_ESTIMATE); */
226
227
228			/* for (i = 0; i < N; i++) { */
229			/* present_in_old[i] = 0; */
230			/* } */
231
232			for (i=0;i<2*state->nx;i=i+2){
233			for(j=0;j<2*state->ny;j++){
234			newy = j;
235			if ((j % 2) == 0) {
236			newx = i;
237			} else {
238			newx = i + 1;
239			}
240			newindex = newx2state->ny + newy;
241
242			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
243			lineCount++;
244			if(eof_test == NULL){
245			printf("error in reading file at line: %d\n", lineCount);
246			exit(8);
247			}
248
249			foo = strtok(read_buffer, " ,;\t\0");
250			(void)strcpy(state->r[newindex].name, foo); //copy the atom name
251
252			// next we grab the positions
253
254			foo = strtok(NULL, " ,;\t\0");
255			if(foo == NULL){
256			printf("error in reading postition x from %s\n"
257			"natoms = %d, line = %d\n",
258			in_name, state->nAtoms, lineCount );
259			exit(8);
260			}
261			(void)sscanf( foo, "%lf", &state->r[newindex].pos[0] );
262
263			foo = strtok(NULL, " ,;\t\0");
264			if(foo == NULL){
265			printf("error in reading postition y from %s\n"
266			"natoms = %d, line = %d\n",
267			in_name, state->nAtoms, lineCount );
268			exit(8);
269			}
270			(void)sscanf( foo, "%lf", &state->r[newindex].pos[1] );
271
272			foo = strtok(NULL, " ,;\t\0");
273			if(foo == NULL){
274			printf("error in reading postition z from %s\n"
275			"natoms = %d, line = %d\n",
276			in_name, state->nAtoms, lineCount );
277			exit(8);
278			}
279			(void)sscanf( foo, "%lf", &state->r[newindex].pos[2] );
280			sumz += state->r[newindex].pos[2];
281			sumz2 += pow(state->r[newindex].pos[2],2);
282
283			foo = strtok(NULL, " ,;\t\0");
284			if(foo == NULL){
285			printf("error in reading angle phi from %s\n"
286			"natoms = %d, line = %d\n",
287			in_name, state->nAtoms, lineCount );
288			exit(8);
289			}
290			(void)sscanf( foo, "%lf", &state->r[newindex].phi );
291
292			foo = strtok(NULL, " ,;\t\0");
293			if(foo == NULL){
294			printf("error in reading unit vector x from %s\n"
295			"natoms = %d, line = %d\n",
296			in_name, state->nAtoms, lineCount );
297			exit(8);
298			}
299			(void)sscanf( foo, "%lf", &uxi );
300
301			foo = strtok(NULL, " ,;\t\0");
302			if(foo == NULL){
303			printf("error in reading unit vector y from %s\n"
304			"natoms = %d, line = %d\n",
305			in_name, state->nAtoms, lineCount );
306			exit(8);
307			}
308			(void)sscanf( foo, "%lf", &uyi );
309
310			foo = strtok(NULL, " ,;\t\0");
311			if(foo == NULL){
312			printf("error in reading unit vector z from %s\n"
313			"natoms = %d, line = %d\n",
314			in_name, state->nAtoms, lineCount );
315			exit(8);
316			}
317			(void)sscanf( foo, "%lf", &uzi );
318
319			state->r[newindex].theta = acos(uzi);
320
321			// The one parameter not stored in the dump file is the dipole strength
322			state->r[newindex].mu = state->strength;
323
324			present_in_old[newindex] = 1;
325			}
326			}
327			averz = sumz / state->nAtoms;
328			averz2 = sumz2 / state->nAtoms;
329			sigmaz = sqrt(averz2 - pow(averz,2));
330			printf("%lf\t%lf\t%lf\n", averz, averz2, sigmaz);
331
332			/* for (i=0; i< 2state->nx; i++) { /
333			/* for(j=0; j< 2state->ny; j++) { /
334			/* newindex = i2state->ny + j; */
335			/* mag[newindex] = 0.0; */
336			/* newmag[newindex] = 0.0; */
337			/* } */
338			/* } */
339
340			/* for (i=0; i< 2state->nx; i++) { /
341			/* for(j=0; j< 2state->ny; j++) { /
342			/* newindex = i2state->ny + j; */
343			/* if (present_in_old[newindex] == 0) { */
344			/* // interpolate from surrounding points: */
345
346			/* interpsum = 0.0; */
347			/* ninterp = 0; */
348
349			/* //point1 = bottom; */
350
351			/* px = i - 1; */
352			/* py = j; */
353			/* newp = px2state->ny + py; */
354			/* if (px >= 0) { */
355			/* interpsum += state->r[newp].pos[2]; */
356			/* ninterp++; */
357			/* state->r[newindex].pos[1] = state->r[newp].pos[1]; */
358			/* } */
359
360			/* //point2 = top; */
361
362			/* px = i + 1; */
363			/* py = j; */
364			/* newp = px2state->ny + py; */
365			/* if (px < 2state->nx) { /
366			/* interpsum += state->r[newp].pos[2]; */
367			/* ninterp++; */
368			/* state->r[newindex].pos[1] = state->r[newp].pos[1]; */
369			/* } */
370
371			/* //point3 = left; */
372
373			/* px = i; */
374			/* py = j - 1; */
375			/* newp = px2state->ny + py; */
376			/* if (py >= 0) { */
377			/* interpsum += state->r[newp].pos[2]; */
378			/* ninterp++; */
379			/* state->r[newindex].pos[0] = state->r[newp].pos[0]; */
380			/* } */
381
382			/* //point4 = right; */
383
384			/* px = i; */
385			/* py = j + 1; */
386			/* newp = px2state->ny + py; */
387			/* if (py < 2state->ny) { /
388			/* interpsum += state->r[newp].pos[2]; */
389			/* ninterp++; */
390			/* state->r[newindex].pos[0] = state->r[newp].pos[0]; */
391			/* } */
392
393			/* value = interpsum / (double)ninterp; */
394
395			/* state->r[newindex].pos[2] = value; */
396			/* } */
397			/* } */
398			/* } */
399
400
401
402
403
404			/* for (i=0; i < 2state->nx; i++) { /
405			/* for (j=0; j < 2state->ny; j++) { /
406			/* newindex = i2state->ny + j; */
407
408			/* c_re(in[newindex]) = state->r[newindex].pos[2]; */
409			/* c_im(in[newindex]) = 0.0; */
410			/* } */
411			/* } */
412
413
414			/* fftwnd(p, 1, in, 1, 0, out, 1, 0); */
415
416			/* for (i=0; i< 2state->nx; i++) { /
417			/* for(j=0; j< 2state->ny; j++) { /
418			/* newindex = i2state->ny + j; */
419			/* mag[newindex] += sqrt(pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2)); */
420			/* } */
421			/* } */
422
423
424			temp_state = state->next;
425			state->next = NULL;
426
427			if (temp_state != NULL) {
428			free(temp_state->r);
429			free(temp_state);
430			}
431
432
433			/* fftwnd_destroy_plan(p); */
434			/* fftw_free(out); */
435			/* fftw_free(in); */
436
437			free(present_in_old);
438
439			// Make a new frame
440
441			temp_state = (struct system *) malloc(sizeof(struct system));
442			temp_state->next = state;
443			state = temp_state;
444			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
445			lineCount++;
446
447			}
448
449
450			// This stuff is for printing out the FFT results directly for matlab:
451
452			/* for (i=0; i < 2nx; i++) { /
453			/* for(j=0; j< 2ny; j++) { /
454			/* newindex = i2ny + j; */
455			/* printf("%lf\t", mag[newindex]/(double)nframes); */
456			/* } */
457			/* printf("\n"); */
458			/* } */
459
460
461			// This stuff is for stitching the four quadrants together:
462
463			/* for (i=0; i< (2nx/2); i++) { /
464			/* for(j=0; j< (2ny/2); j++) { /
465			/* index = i2ny + j; */
466			/* new_i = i + (2nx/2); /
467			/* new_j = j + (2ny/2); /
468			/* new_index = new_i2ny + new_j; */
469			/* newmag[new_index] = mag[index]; */
470			/* } */
471			/* } */
472
473			/* for (i=(2nx/2); i< 2nx; i++) { */
474			/* for(j=0; j< (2ny/2); j++) { /
475			/* index = i2ny + j; */
476			/* new_i = i - (2nx/2); /
477			/* new_j = j + (2ny/2); /
478			/* new_index = new_i2ny + new_j; */
479			/* newmag[new_index] = mag[index]; */
480			/* } */
481			/* } */
482
483			/* for (i=0; i< (2nx/2); i++) { /
484			/* for(j=(2ny/2); j< 2ny; j++) { */
485			/* index = i2ny + j; */
486			/* new_i = i + (2nx/2); /
487			/* new_j = j - (2ny/2); /
488			/* new_index = new_i2ny + new_j; */
489			/* newmag[new_index] = mag[index]; */
490			/* } */
491			/* } */
492
493			/* for (i=(2nx/2); i< 2nx; i++) { */
494			/* for(j=(2ny/2); j< 2ny; j++) { */
495			/* index = i2ny + j; */
496			/* new_i = i - (2nx/2); /
497			/* new_j = j - (2ny/2); /
498			/* new_index = new_i2ny + new_j; */
499			/* newmag[new_index] = mag[index]; */
500			/* } */
501			/* } */
502
503			/* maxfreqx = 2.0 / aLat; */
504			/* maxfreqy = 2.0 / bLat; */
505
506			/* // printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy); */
507
508			/* maxfreq = sqrt(maxfreqxmaxfreqx + maxfreqymaxfreqy); */
509			/* dfreq = maxfreq/(double)(nbins-1); */
510
511			/* //printf("%lf\n", dfreq); */
512
513			/* zero_freq_x = nx; */
514			/* zero_freq_y = ny; */
515
516			/* for (i=0; i< 2nx; i++) { /
517			/* for(j=0; j< 2ny; j++) { /
518
519			/* freq_x = (double)(i - zero_freq_x)maxfreqx / nx; /
520			/* freq_y = (double)(j - zero_freq_y)maxfreqy / ny; /
521
522			/* freq = sqrt(freq_xfreq_x + freq_yfreq_y); */
523
524			/* whichbin = (int) (freq / dfreq); */
525			/* newindex = i2ny + j; */
526			/* // printf("%d %d %lf %lf\n", whichbin, newindex, freq, dfreq); */
527			/* bin[whichbin] += mag[newindex]; */
528			/* samples[whichbin]++; */
529			/* } */
530			/* } */
531
532
533
534
535			/* for (i = 0; i < nbins; i++) { */
536			/* if (samples[i] > 0) { */
537			/* printf("%lf\t%lf\n", idfreq, bin[i]/(double)samples[i]); /
538			/* } */
539			/* } */
540
541			// This stuff is for stitching the four quadrants together:
542
543			/* for (i=0; i< (2nx/2); i++) { /
544			/* for(j=0; j< (2ny/2); j++) { /
545			/* index = i2ny + j; */
546			/* new_i = i + (2nx/2); /
547			/* new_j = j + (2ny/2); /
548			/* new_index = new_i2ny + new_j; */
549			/* newmag[new_index] = mag[index]; */
550			/* } */
551			/* } */
552
553			/* for (i=(2nx/2); i< 2nx; i++) { */
554			/* for(j=0; j< (2ny/2); j++) { /
555			/* index = i2ny + j; */
556			/* new_i = i - (2nx/2); /
557			/* new_j = j + (2ny/2); /
558			/* new_index = new_i2ny + new_j; */
559			/* newmag[new_index] = mag[index]; */
560			/* } */
561			/* } */
562
563			/* for (i=0; i< (2nx/2); i++) { /
564			/* for(j=(2ny/2); j< 2ny; j++) { */
565			/* index = i2ny + j; */
566			/* new_i = i + (2nx/2); /
567			/* new_j = j - (2ny/2); /
568			/* new_index = new_i2ny + new_j; */
569			/* newmag[new_index] = mag[index]; */
570			/* } */
571			/* } */
572
573			/* for (i=(2nx/2); i< 2nx; i++) { */
574			/* for(j=(2ny/2); j< 2ny; j++) { */
575			/* index = i2ny + j; */
576			/* new_i = i - (2nx/2); /
577			/* new_j = j - (2ny/2); /
578			/* new_index = new_i2ny + new_j; */
579			/* newmag[new_index] = mag[index]; */
580			/* } */
581			/* } */
582
583			/* for (i=0; i < 2nx; i++) { /
584			/* for(j=0; j< 2ny; j++) { /
585			/* newindex = i2ny + j; */
586			/* printf("%lf\t",newmag[newindex]/(double)nframes); */
587			/* } */
588			/* printf("\n"); */
589			/* } */
590
591			// This stuff is for printing out the FFT results directly:
592
593			/* for (i=0; i < 2nx; i++) { /
594			/* for(j=0; j< 2ny; j++) { /
595			/* newindex = i2ny + j; */
596			/* printf("%lf\t", mag[newindex]/(double)nframes); */
597			/* } */
598			/* printf("\n"); */
599			/* } */
600
601			fclose(in_file);
602			free(newmag);
603			free(mag);
604			free(samples);
605			free(bin);
606			free(state);
607			return(0);
608			}
609
610			/* } */
611			/* for(i=0;i<120;i++) printf("%f\t%f\n", c_re(out[i]), c_im(out[i])); */
612			/* fftw_destroy_plan(p); */
613			/* fftw_free(in); */
614			/* fftw_free(out); */
615			/* averZ = sumZ / (double)nloops; */
616			/* averUx = sumUx / (double)nloops; */
617			/* averUy = sumUy / (double)nloops; */
618			/* averUz = sumUz / (double)nloops; */
619			/* averP = 1.5sumP / (double)nloops-0.5; /
620			/* averTheta = acos(averUz); */
621			/* printf("nloops=%d\n",nloops); */
622			/* printf("sumZ=%f\n",sumZ); */
623			/* printf("average height is : %f\n",averZ); */
624			/* printf("average ux is : %f\n",averUx); */
625			/* printf("average uy is : %f\n",averUy); */
626			/* printf("average uz is : %f\n",averUz); */
627			/* printf("average angle is : %f\n",averTheta); */
628			/* printf("average p is : %f\n",averP); */
629
630
631			/* return 0; */
632			/* }*/