trunk/ripple/read.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;

  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;


  while(eof_test != NULL){
    
    nframes++;
    (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] );
        
        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;
      }
    }
      
    
    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, 3 months ago) by xsun
Content type:	text/plain
File size:	15928 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
76			int done;
77			char current_flag;
78
79			program_name = argv[0];
80
81			for(i = 1; i < argc; i++){
82			if(argv[i][0] == '-'){
83			done = 0;
84			j = 1;
85			current_flag = argv[i][j];
86			while( (current_flag != '\0') && ( !done ) ){
87			switch(current_flag){
88			case 'i':
89			i++;
90			strcpy( in_name, argv[i] );
91			done = 1;
92			break;
93			}
94			}
95			}
96			}
97
98			struct system* state;
99			struct system* temp_state;
100			struct coords* r;
101
102			FFTW_COMPLEX in, out;
103			fftwnd_plan p;
104
105			nbins = 100;
106			bin = (double ) malloc(sizeof(double) nbins);
107			samples = (int ) malloc(sizeof(int) nbins);
108			for (i=0; i < nbins; i++) {
109			bin[i] = 0.0;
110			samples[i] = 0;
111			}
112
113			in_file = fopen(in_name, "r");
114			if(in_file == NULL){
115			printf("Cannot open file \"%s\" for reading.\n", in_name);
116			exit(8);
117			}
118
119			nframes = 0;
120
121			// start reading the first frame
122
123			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
124			lineCount++;
125
126			state = (struct system *) malloc(sizeof(struct system));
127			state->next = NULL;
128			state->strength = 10.0;
129
130
131			while(eof_test != NULL){
132
133			nframes++;
134			(void)sscanf(read_buffer, "%d", &state->nAtoms);
135			N = 2 * state->nAtoms;
136
137			state->r = (struct coords *)calloc(N, sizeof(struct coords));
138			mag = (double ) malloc(sizeof(double) N);
139			newmag = (double ) malloc(sizeof(double) N);
140			present_in_old = (int ) malloc(sizeof(int) N);
141
142			// read and the comment line and grab the time and box dimensions
143
144			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
145			lineCount++;
146			if(eof_test == NULL){
147			printf("error in reading file at line: %d\n", lineCount);
148			exit(8);
149			}
150
151			foo = strtok( read_buffer, " ,;\t\n" );
152			(void)sscanf( read_buffer, "%d", &state->iCycle );
153
154			foo = strtok(NULL, " ,;\t\0");
155			if(foo == NULL){
156			printf("error in reading file at line: %d\n", lineCount);
157			exit(8);
158			}
159			(void)sscanf(foo, "%d", &state->nx);
160
161			nx = state->nx;
162
163			foo = strtok(NULL, " ,;\t\0");
164			if(foo == NULL){
165			printf("error in reading file at line: %d\n", lineCount);
166			exit(8);
167			}
168			(void)sscanf(foo, "%d", &state->ny);
169
170			ny = state->ny;
171
172			foo = strtok(NULL, " ,;\t\0");
173			if(foo == NULL){
174			printf("error in reading file at line: %d\n", lineCount);
175			exit(8);
176			}
177			(void)sscanf(foo, "%lf",&state->Hmat[0][0]);
178
179			foo = strtok(NULL, " ,;\t\0");
180			if(foo == NULL){
181			printf("error in reading file at line: %d\n", lineCount);
182			exit(8);
183			}
184			(void)sscanf(foo, "%lf",&state->Hmat[1][0]);
185
186			foo = strtok(NULL, " ,;\t\0");
187			if(foo == NULL){
188			printf("error in reading file at line: %d\n", lineCount);
189			exit(8);
190			}
191			(void)sscanf(foo, "%lf",&state->Hmat[0][1]);
192
193			foo = strtok(NULL, " ,;\t\0");
194			if(foo == NULL){
195			printf("error in reading file at line: %d\n", lineCount);
196			exit(8);
197			}
198			(void)sscanf(foo, "%lf",&state->Hmat[1][1]);
199
200			// Length of the two box vectors:
201
202			dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
203			dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));
204
205			aLat = dx / (double)(state->nx);
206			bLat = dy / (double)(state->ny);
207
208
209			// FFT stuff depends on nx and ny, so delay allocation until we have
210			// that information
211
212			in = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
213			out = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
214			p = fftw2d_create_plan(2*state->nx,
215			2*state->ny,
216			FFTW_FORWARD,
217			FFTW_ESTIMATE);
218
219
220			for (i = 0; i < N; i++) {
221			present_in_old[i] = 0;
222			}
223
224			for (i=0;i<2*state->nx;i=i+2){
225			for(j=0;j<2*state->ny;j++){
226			newy = j;
227			if ((j % 2) == 0) {
228			newx = i;
229			} else {
230			newx = i + 1;
231			}
232			newindex = newx2state->ny + newy;
233
234			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
235			lineCount++;
236			if(eof_test == NULL){
237			printf("error in reading file at line: %d\n", lineCount);
238			exit(8);
239			}
240
241			foo = strtok(read_buffer, " ,;\t\0");
242			(void)strcpy(state->r[newindex].name, foo); //copy the atom name
243
244			// next we grab the positions
245
246			foo = strtok(NULL, " ,;\t\0");
247			if(foo == NULL){
248			printf("error in reading postition x from %s\n"
249			"natoms = %d, line = %d\n",
250			in_name, state->nAtoms, lineCount );
251			exit(8);
252			}
253			(void)sscanf( foo, "%lf", &state->r[newindex].pos[0] );
254
255			foo = strtok(NULL, " ,;\t\0");
256			if(foo == NULL){
257			printf("error in reading postition y from %s\n"
258			"natoms = %d, line = %d\n",
259			in_name, state->nAtoms, lineCount );
260			exit(8);
261			}
262			(void)sscanf( foo, "%lf", &state->r[newindex].pos[1] );
263
264			foo = strtok(NULL, " ,;\t\0");
265			if(foo == NULL){
266			printf("error in reading postition z from %s\n"
267			"natoms = %d, line = %d\n",
268			in_name, state->nAtoms, lineCount );
269			exit(8);
270			}
271			(void)sscanf( foo, "%lf", &state->r[newindex].pos[2] );
272
273			foo = strtok(NULL, " ,;\t\0");
274			if(foo == NULL){
275			printf("error in reading angle phi from %s\n"
276			"natoms = %d, line = %d\n",
277			in_name, state->nAtoms, lineCount );
278			exit(8);
279			}
280			(void)sscanf( foo, "%lf", &state->r[newindex].phi );
281
282			foo = strtok(NULL, " ,;\t\0");
283			if(foo == NULL){
284			printf("error in reading unit vector x from %s\n"
285			"natoms = %d, line = %d\n",
286			in_name, state->nAtoms, lineCount );
287			exit(8);
288			}
289			(void)sscanf( foo, "%lf", &uxi );
290
291			foo = strtok(NULL, " ,;\t\0");
292			if(foo == NULL){
293			printf("error in reading unit vector y from %s\n"
294			"natoms = %d, line = %d\n",
295			in_name, state->nAtoms, lineCount );
296			exit(8);
297			}
298			(void)sscanf( foo, "%lf", &uyi );
299
300			foo = strtok(NULL, " ,;\t\0");
301			if(foo == NULL){
302			printf("error in reading unit vector z from %s\n"
303			"natoms = %d, line = %d\n",
304			in_name, state->nAtoms, lineCount );
305			exit(8);
306			}
307			(void)sscanf( foo, "%lf", &uzi );
308
309			state->r[newindex].theta = acos(uzi);
310
311			// The one parameter not stored in the dump file is the dipole strength
312			state->r[newindex].mu = state->strength;
313
314			present_in_old[newindex] = 1;
315			}
316			}
317
318
319			for (i=0; i< 2*state->nx; i++) {
320			for(j=0; j< 2*state->ny; j++) {
321			newindex = i2state->ny + j;
322			mag[newindex] = 0.0;
323			newmag[newindex] = 0.0;
324			}
325			}
326
327			for (i=0; i< 2*state->nx; i++) {
328			for(j=0; j< 2*state->ny; j++) {
329			newindex = i2state->ny + j;
330			if (present_in_old[newindex] == 0) {
331			// interpolate from surrounding points:
332
333			interpsum = 0.0;
334			ninterp = 0;
335
336			//point1 = bottom;
337
338			px = i - 1;
339			py = j;
340			newp = px2state->ny + py;
341			if (px >= 0) {
342			interpsum += state->r[newp].pos[2];
343			ninterp++;
344			state->r[newindex].pos[1] = state->r[newp].pos[1];
345			}
346
347			//point2 = top;
348
349			px = i + 1;
350			py = j;
351			newp = px2state->ny + py;
352			if (px < 2*state->nx) {
353			interpsum += state->r[newp].pos[2];
354			ninterp++;
355			state->r[newindex].pos[1] = state->r[newp].pos[1];
356			}
357
358			//point3 = left;
359
360			px = i;
361			py = j - 1;
362			newp = px2state->ny + py;
363			if (py >= 0) {
364			interpsum += state->r[newp].pos[2];
365			ninterp++;
366			state->r[newindex].pos[0] = state->r[newp].pos[0];
367			}
368
369			//point4 = right;
370
371			px = i;
372			py = j + 1;
373			newp = px2state->ny + py;
374			if (py < 2*state->ny) {
375			interpsum += state->r[newp].pos[2];
376			ninterp++;
377			state->r[newindex].pos[0] = state->r[newp].pos[0];
378			}
379
380			value = interpsum / (double)ninterp;
381
382			state->r[newindex].pos[2] = value;
383			}
384			}
385			}
386
387
388
389
390
391			for (i=0; i < 2*state->nx; i++) {
392			for (j=0; j < 2*state->ny; j++) {
393			newindex = i2state->ny + j;
394
395			c_re(in[newindex]) = state->r[newindex].pos[2];
396			c_im(in[newindex]) = 0.0;
397			}
398			}
399
400
401			fftwnd(p, 1, in, 1, 0, out, 1, 0);
402
403			for (i=0; i< 2*state->nx; i++) {
404			for(j=0; j< 2*state->ny; j++) {
405			newindex = i2state->ny + j;
406			mag[newindex] += sqrt(pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2));
407			}
408			}
409
410
411			temp_state = state->next;
412			state->next = NULL;
413
414			if (temp_state != NULL) {
415			free(temp_state->r);
416			free(temp_state);
417			}
418
419
420			fftwnd_destroy_plan(p);
421			fftw_free(out);
422			fftw_free(in);
423
424			free(present_in_old);
425
426			// Make a new frame
427
428			temp_state = (struct system *) malloc(sizeof(struct system));
429			temp_state->next = state;
430			state = temp_state;
431			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
432			lineCount++;
433
434			}
435
436
437			// This stuff is for printing out the FFT results directly for matlab:
438
439			/* for (i=0; i < 2nx; i++) { /
440			/* for(j=0; j< 2ny; j++) { /
441			/* newindex = i2ny + j; */
442			/* printf("%lf\t", mag[newindex]/(double)nframes); */
443			/* } */
444			/* printf("\n"); */
445			/* } */
446
447
448			// This stuff is for stitching the four quadrants together:
449
450			for (i=0; i< (2*nx/2); i++) {
451			for(j=0; j< (2*ny/2); j++) {
452			index = i2ny + j;
453			new_i = i + (2*nx/2);
454			new_j = j + (2*ny/2);
455			new_index = new_i2ny + new_j;
456			newmag[new_index] = mag[index];
457			}
458			}
459
460			for (i=(2nx/2); i< 2nx; i++) {
461			for(j=0; j< (2*ny/2); j++) {
462			index = i2ny + j;
463			new_i = i - (2*nx/2);
464			new_j = j + (2*ny/2);
465			new_index = new_i2ny + new_j;
466			newmag[new_index] = mag[index];
467			}
468			}
469
470			for (i=0; i< (2*nx/2); i++) {
471			for(j=(2ny/2); j< 2ny; j++) {
472			index = i2ny + j;
473			new_i = i + (2*nx/2);
474			new_j = j - (2*ny/2);
475			new_index = new_i2ny + new_j;
476			newmag[new_index] = mag[index];
477			}
478			}
479
480			for (i=(2nx/2); i< 2nx; i++) {
481			for(j=(2ny/2); j< 2ny; j++) {
482			index = i2ny + j;
483			new_i = i - (2*nx/2);
484			new_j = j - (2*ny/2);
485			new_index = new_i2ny + new_j;
486			newmag[new_index] = mag[index];
487			}
488			}
489
490			maxfreqx = 2.0 / aLat;
491			maxfreqy = 2.0 / bLat;
492
493			// printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
494
495			maxfreq = sqrt(maxfreqxmaxfreqx + maxfreqymaxfreqy);
496			dfreq = maxfreq/(double)(nbins-1);
497
498			//printf("%lf\n", dfreq);
499
500			zero_freq_x = nx;
501			zero_freq_y = ny;
502
503			for (i=0; i< 2*nx; i++) {
504			for(j=0; j< 2*ny; j++) {
505
506			freq_x = (double)(i - zero_freq_x)*maxfreqx / nx;
507			freq_y = (double)(j - zero_freq_y)*maxfreqy / ny;
508
509			freq = sqrt(freq_xfreq_x + freq_yfreq_y);
510
511			whichbin = (int) (freq / dfreq);
512			newindex = i2ny + j;
513			// printf("%d %d %lf %lf\n", whichbin, newindex, freq, dfreq);
514			bin[whichbin] += mag[newindex];
515			samples[whichbin]++;
516			}
517			}
518
519
520
521
522			for (i = 0; i < nbins; i++) {
523			if (samples[i] > 0) {
524			printf("%lf\t%lf\n", i*dfreq, bin[i]/(double)samples[i]);
525			}
526			}
527
528			// This stuff is for stitching the four quadrants together:
529
530			/* for (i=0; i< (2nx/2); i++) { /
531			/* for(j=0; j< (2ny/2); j++) { /
532			/* index = i2ny + j; */
533			/* new_i = i + (2nx/2); /
534			/* new_j = j + (2ny/2); /
535			/* new_index = new_i2ny + new_j; */
536			/* newmag[new_index] = mag[index]; */
537			/* } */
538			/* } */
539
540			/* for (i=(2nx/2); i< 2nx; i++) { */
541			/* for(j=0; j< (2ny/2); j++) { /
542			/* index = i2ny + j; */
543			/* new_i = i - (2nx/2); /
544			/* new_j = j + (2ny/2); /
545			/* new_index = new_i2ny + new_j; */
546			/* newmag[new_index] = mag[index]; */
547			/* } */
548			/* } */
549
550			/* for (i=0; i< (2nx/2); i++) { /
551			/* for(j=(2ny/2); j< 2ny; j++) { */
552			/* index = i2ny + j; */
553			/* new_i = i + (2nx/2); /
554			/* new_j = j - (2ny/2); /
555			/* new_index = new_i2ny + new_j; */
556			/* newmag[new_index] = mag[index]; */
557			/* } */
558			/* } */
559
560			/* for (i=(2nx/2); i< 2nx; i++) { */
561			/* for(j=(2ny/2); j< 2ny; j++) { */
562			/* index = i2ny + j; */
563			/* new_i = i - (2nx/2); /
564			/* new_j = j - (2ny/2); /
565			/* new_index = new_i2ny + new_j; */
566			/* newmag[new_index] = mag[index]; */
567			/* } */
568			/* } */
569
570			/* for (i=0; i < 2nx; i++) { /
571			/* for(j=0; j< 2ny; j++) { /
572			/* newindex = i2ny + j; */
573			/* printf("%lf\t",newmag[newindex]/(double)nframes); */
574			/* } */
575			/* printf("\n"); */
576			/* } */
577
578			// This stuff is for printing out the FFT results directly:
579
580			/* for (i=0; i < 2nx; i++) { /
581			/* for(j=0; j< 2ny; j++) { /
582			/* newindex = i2ny + j; */
583			/* printf("%lf\t", mag[newindex]/(double)nframes); */
584			/* } */
585			/* printf("\n"); */
586			/* } */
587
588			fclose(in_file);
589			free(newmag);
590			free(mag);
591			free(samples);
592			free(bin);
593			free(state);
594			return(0);
595			}
596
597			/* } */
598			/* for(i=0;i<120;i++) printf("%f\t%f\n", c_re(out[i]), c_im(out[i])); */
599			/* fftw_destroy_plan(p); */
600			/* fftw_free(in); */
601			/* fftw_free(out); */
602			/* averZ = sumZ / (double)nloops; */
603			/* averUx = sumUx / (double)nloops; */
604			/* averUy = sumUy / (double)nloops; */
605			/* averUz = sumUz / (double)nloops; */
606			/* averP = 1.5sumP / (double)nloops-0.5; /
607			/* averTheta = acos(averUz); */
608			/* printf("nloops=%d\n",nloops); */
609			/* printf("sumZ=%f\n",sumZ); */
610			/* printf("average height is : %f\n",averZ); */
611			/* printf("average ux is : %f\n",averUx); */
612			/* printf("average uy is : %f\n",averUy); */
613			/* printf("average uz is : %f\n",averUz); */
614			/* printf("average angle is : %f\n",averTheta); */
615			/* printf("average p is : %f\n",averP); */
616
617
618			/* return 0; */
619			/* }*/