trunk/dp/read.c

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

int main()
{
  FILE *inFile;
  char *endptr, s[100000];
  double *x, *y, *z;
  double *mag, *newmag;
  int *present_in_old;
  double ux, uy, uz, p1;
  double sumZ, sumUx, sumUy, sumUz, sumP;
  double interpsum, value;
  int ninterp, px, py, newp;
  int i, j, nsites, nloops, xlat, ylat;
  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, pt1x, pt1y, pt2x, pt2y;
  int whichbin;
  int nbins;
  int *samples;
  double *bin;

  FFTW_COMPLEX *in, *out;
  fftwnd_plan p;
        
  nbins = 200;

  inFile = fopen("./dp_file","r");
  if (inFile == NULL) {
    printf("Error opening file\n");
    exit(-1);
  }

  fgets(s, 500, inFile);
  nsites = atoi(s);
  N = 2*nsites;
  fscanf(inFile,"%s",s);
  xlat=atoi(s);
  fscanf(inFile,"%s",s);
  ylat=atoi(s);
 
  
  //sanity check:
  if (2*nsites != xlat * ylat) {
    printf("Lattice Mismatch: nsites = %d, xlat = %d, ylat = %d\n", nsites, xlat, ylat);
    exit(-1);
  }

  rewind(inFile);

  in = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
  out = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
  p =  fftw2d_create_plan(xlat, ylat, FFTW_FORWARD, FFTW_ESTIMATE);

  x = (double *) malloc(sizeof(double) * N);
  y = (double *) malloc(sizeof(double) * N);
  z = (double *) malloc(sizeof(double) * N);
  mag = (double *) malloc(sizeof(double) * N);
  newmag = (double *) malloc(sizeof(double) * N);
  present_in_old = (int *) malloc(sizeof(int) * N);

  bin = (double *) malloc(sizeof(double) * nbins);
  samples = (int *) malloc(sizeof(int) * nbins);


  for (i=0; i< xlat; i++) {
    for(j=0; j< ylat; j++) {
      newindex = i*ylat + j;
      mag[newindex] = 0.0;
      newmag[newindex] = 0.0;
    }
  }

  sumZ = 0.0;
  sumUx = 0.0;
  sumUy = 0.0;
  sumUz = 0.0;
  sumP = 0.0;

  nloops = 0;
  nframes = 0;
        
  while(!feof(inFile))
    {
      nframes++;
      fgets(s, 500, inFile);
      nsites = atoi(s);
      fscanf(inFile,"%s",s);
      xlat=atoi(s);
      fscanf(inFile,"%s",s);
      ylat=atoi(s);

      //sanity check:
      if (2*nsites != xlat * ylat) {
        printf("Lattice Mismatch: nsites = %d, xlat = %d, ylat = %d\n", nsites, xlat, ylat);
        exit(-1);
      }
      
      for (i = 0; i < N; i++ ) {
        present_in_old[i] = 0;
      }

      for (i=0;i<xlat;i=i+2)
        {
          for(j=0;j<ylat;j++)
            {

              newy = j;
              if ((j % 2) == 0) {
                newx = i;
              } else {
                newx = i + 1;
              }
              newindex = newx*ylat + newy;

              fscanf(inFile,"%s",s);
              fscanf(inFile,"%s",s);
              x[newindex]=strtod(s,&endptr);
              fscanf(inFile,"%s",s);
              y[newindex]=strtod(s,&endptr);
              fscanf(inFile,"%s",s);
              z[newindex] = strtod(s,&endptr);
              present_in_old[newindex] = 1;
//              printf("z=%f\t",z[newindex]);
              fscanf(inFile,"%s",s);
              fscanf(inFile,"%s",s);
              ux=strtod(s,&endptr);
              fscanf(inFile,"%s",s);
              uy=strtod(s,&endptr);
              fscanf(inFile,"%s",s);
              uz=strtod(s,&endptr);
//            p1 = uz * uz;
//            sumZ = sumZ + z[i];
//            sumUx = sumUx + ux;
//            sumUy = sumUy + uy;
//            sumUz = sumUz + uz;
//            sumP = sumP + p1;
              nloops += 1;
              fgets(s, 500, inFile);
            }
        }
      
      for (i=0; i< xlat; i++) {
        for(j=0; j< ylat; j++) {
          newindex = i*ylat + 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*ylat + py;
            if (px >= 0) {
              interpsum += z[newp];
              ninterp++;
              y[newindex] = y[newp];
            }
            if (px < 0) {
                interpsum += z[newp+xlat*ylat];
                ninterp++;
                y[newindex] = y[newp+xlat*ylat];
            }

            //point2 = top;

            px = i + 1;
            py = j;
            newp = px*ylat + py;
            if (px < xlat) {
              interpsum += z[newp];
              ninterp++;
              y[newindex] = y[newp];
            } 
            if (px == xlat) {
              interpsum += z[newp-xlat*ylat];
              ninterp++;
              y[newindex] = y[newp-xlat*ylat];
            } 

            //point3 = left;

            px = i;
            py = j - 1;
            newp = px*ylat + py;
            if (py >= 0) {
              interpsum += z[newp];
              ninterp++;
              x[newindex] = x[newp];
            } 
            if (py < 0) {
              interpsum += z[newp+ylat];
              ninterp++;
              x[newindex] = x[newp+ylat];
            } 

            //point4 = right;

            px = i;
            py = j + 1;
            newp = px*ylat + py;
            if (py < ylat) {
              interpsum += z[newp];
              ninterp++;
              x[newindex] = x[newp];
            } 
            if (py == ylat) {
              interpsum += z[newp-ylat];
              ninterp++;
              x[newindex] = x[newp-ylat];
            } 

            value = interpsum / (double)ninterp;

            z[newindex] = value;

          }
        }
      }

      for (i=0; i < xlat; i++) {
        for (j=0; j < ylat; j++) {
                newindex = i*ylat + j;
                c_re(in[newindex]) = z[newindex];
                c_im(in[newindex]) = 0.0;
        } 
      }

      fftwnd(p, 1, in, 1, 0, out, 1, 0); 
      for (i=0; i< xlat; i++) {
        for(j=0; j< ylat; j++) {
          newindex = i*ylat + j;
          mag[newindex] += sqrt(pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2));
        }
      }
      
    }

  for (i=0; i< (xlat/2); i++) {
    for(j=0; j< (ylat/2); j++) {
      index = i*ylat + j;
      new_i = i + (xlat/2);
      new_j = j + (ylat/2);
      new_index = new_i*ylat + new_j;
      newmag[new_index] = mag[index];
    }
  }
  
  for (i=(xlat/2); i< xlat; i++) {
    for(j=0; j< (ylat/2); j++) {
      index = i*ylat + j;
      new_i = i - (xlat/2);
      new_j = j + (ylat/2);
      new_index = new_i*ylat + new_j;
      newmag[new_index] = mag[index];
    }
  }
    
  for (i=0; i< (xlat/2); i++) {
    for(j=(ylat/2); j< ylat; j++) {
      index = i*ylat + j;
      new_i = i + (xlat/2);
      new_j = j - (ylat/2);
      new_index = new_i*ylat + new_j;
      newmag[new_index] = mag[index];
    }
  }
      
  for (i=(xlat/2); i< xlat; i++) {
    for(j=(ylat/2); j< ylat; j++) {
      index = i*ylat + j;
      new_i = i - (xlat/2);
      new_j = j - (ylat/2);
      new_index = new_i*ylat + new_j;
      newmag[new_index] = mag[index];
    }
  }
      
  pt1x = x[0];
  pt1y = y[0];
  pt2x = x[ylat + 1];
  pt2y = y[ylat + 1];

  dx = pt2x - pt1x;
  dy = pt2y - pt1y;

  maxfreqx = 1.0/dx;
  maxfreqy = 1.0/dy;


  maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
  dfreq = maxfreq/(double)nbins;

  for (i=0; i < nbins; i++) {
    bin[i] = 0.0;
    samples[i] = 0;
  }

  zero_freq_x = xlat/2;
  zero_freq_y = ylat/2;

  for (i=0; i< xlat; i++) {
    for(j=0; j< ylat; j++) {

      freq_x = (double)(i - zero_freq_x)*maxfreqx/(double)xlat;
      freq_y = (double)(j - zero_freq_y)*maxfreqy/(double)ylat;

      freq = sqrt(freq_x*freq_x + freq_y*freq_y);

      whichbin = (int) (freq / dfreq);
      newindex = i*ylat + j;
      bin[whichbin] += newmag[newindex];
      samples[whichbin]++;

      //newindex = i*ylat + j;
      //printf("%lf\t", newmag[newindex]/(double)nframes);
    }
    //printf("\n");
  }
  

  for (i = 0; i < nbins; i++) {

    if (samples[i] > 0) {
      printf("%lf\t%lf\n", i*dfreq, bin[i]/(double)samples[i]);
    }
  }


  free(samples);
  free(bin);
  free(x);
  free(y);
  free(z);
  free(mag);
  free(newmag);
  free(present_in_old);
  fftwnd_destroy_plan(p);
  fftw_free(in);
  fftw_free(out);
  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); */

/*       fclose(inFile); */
/*       return 0; */
/*     } */
Revision:	585
Committed:	Wed Jul 9 15:52:32 2003 UTC (21 years ago) by xsun
Content type:	text/plain
File size:	8330 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	xsun	581	#include<stdio.h>
2			#include<string.h>
3			#include<stdlib.h>
4			#include<math.h>
5			#include<fftw.h>
6
7			int main()
8			{
9			FILE *inFile;
10			char *endptr, s[100000];
11			double x, y, *z;
12			double mag, newmag;
13			int *present_in_old;
14			double ux, uy, uz, p1;
15			double sumZ, sumUx, sumUy, sumUz, sumP;
16			double interpsum, value;
17			int ninterp, px, py, newp;
18			int i, j, nsites, nloops, xlat, ylat;
19			int newx, newy, newindex, index;
20			int new_i, new_j, new_index;
21			int N, nframes;
22			double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
23			double maxfreqx, maxfreqy, maxfreq, dfreq;
24			double dx, dy, pt1x, pt1y, pt2x, pt2y;
25			int whichbin;
26			int nbins;
27			int *samples;
28			double *bin;
29
30			FFTW_COMPLEX in, out;
31			fftwnd_plan p;
32
33			nbins = 200;
34
35			inFile = fopen("./dp_file","r");
36			if (inFile == NULL) {
37			printf("Error opening file\n");
38			exit(-1);
39			}
40
41			fgets(s, 500, inFile);
42			nsites = atoi(s);
43			N = 2*nsites;
44			fscanf(inFile,"%s",s);
45			xlat=atoi(s);
46			fscanf(inFile,"%s",s);
47			ylat=atoi(s);
48
49
50			//sanity check:
51			if (2nsites != xlat ylat) {
52			printf("Lattice Mismatch: nsites = %d, xlat = %d, ylat = %d\n", nsites, xlat, ylat);
53			exit(-1);
54			}
55
56			rewind(inFile);
57
58			in = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
59			out = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
60			p = fftw2d_create_plan(xlat, ylat, FFTW_FORWARD, FFTW_ESTIMATE);
61
62			x = (double ) malloc(sizeof(double) N);
63			y = (double ) malloc(sizeof(double) N);
64			z = (double ) malloc(sizeof(double) N);
65			mag = (double ) malloc(sizeof(double) N);
66			newmag = (double ) malloc(sizeof(double) N);
67			present_in_old = (int ) malloc(sizeof(int) N);
68
69			bin = (double ) malloc(sizeof(double) nbins);
70			samples = (int ) malloc(sizeof(int) nbins);
71
72
73
74			for (i=0; i< xlat; i++) {
75			for(j=0; j< ylat; j++) {
76			newindex = i*ylat + j;
77			mag[newindex] = 0.0;
78			newmag[newindex] = 0.0;
79			}
80			}
81
82			sumZ = 0.0;
83			sumUx = 0.0;
84			sumUy = 0.0;
85			sumUz = 0.0;
86			sumP = 0.0;
87
88			nloops = 0;
89			nframes = 0;
90
91			while(!feof(inFile))
92			{
93			nframes++;
94			fgets(s, 500, inFile);
95			nsites = atoi(s);
96			fscanf(inFile,"%s",s);
97			xlat=atoi(s);
98			fscanf(inFile,"%s",s);
99			ylat=atoi(s);
100
101			//sanity check:
102			if (2nsites != xlat ylat) {
103			printf("Lattice Mismatch: nsites = %d, xlat = %d, ylat = %d\n", nsites, xlat, ylat);
104			exit(-1);
105			}
106
107			for (i = 0; i < N; i++ ) {
108			present_in_old[i] = 0;
109			}
110
111			for (i=0;i<xlat;i=i+2)
112			{
113			for(j=0;j<ylat;j++)
114			{
115
116			newy = j;
117			if ((j % 2) == 0) {
118			newx = i;
119			} else {
120			newx = i + 1;
121			}
122			newindex = newx*ylat + newy;
123
124			fscanf(inFile,"%s",s);
125			fscanf(inFile,"%s",s);
126			x[newindex]=strtod(s,&endptr);
127			fscanf(inFile,"%s",s);
128			y[newindex]=strtod(s,&endptr);
129			fscanf(inFile,"%s",s);
130			z[newindex] = strtod(s,&endptr);
131			present_in_old[newindex] = 1;
132			// printf("z=%f\t",z[newindex]);
133			fscanf(inFile,"%s",s);
134			fscanf(inFile,"%s",s);
135			ux=strtod(s,&endptr);
136			fscanf(inFile,"%s",s);
137			uy=strtod(s,&endptr);
138			fscanf(inFile,"%s",s);
139			uz=strtod(s,&endptr);
140			// p1 = uz * uz;
141			// sumZ = sumZ + z[i];
142			// sumUx = sumUx + ux;
143			// sumUy = sumUy + uy;
144			// sumUz = sumUz + uz;
145			// sumP = sumP + p1;
146			nloops += 1;
147			fgets(s, 500, inFile);
148			}
149			}
150
151			for (i=0; i< xlat; i++) {
152			for(j=0; j< ylat; j++) {
153			newindex = i*ylat + j;
154			if (present_in_old[newindex] == 0) {
155			// interpolate from surrounding points:
156
157			interpsum = 0.0;
158			ninterp = 0;
159
160			//point1 = bottom;
161
162			px = i - 1;
163			py = j;
164			newp = px*ylat + py;
165			if (px >= 0) {
166			interpsum += z[newp];
167			ninterp++;
168			y[newindex] = y[newp];
169			}
170			if (px < 0) {
171			interpsum += z[newp+xlat*ylat];
172			ninterp++;
173			y[newindex] = y[newp+xlat*ylat];
174			}
175
176			//point2 = top;
177
178			px = i + 1;
179			py = j;
180			newp = px*ylat + py;
181			if (px < xlat) {
182			interpsum += z[newp];
183			ninterp++;
184			y[newindex] = y[newp];
185			}
186			if (px == xlat) {
187			interpsum += z[newp-xlat*ylat];
188			ninterp++;
189			y[newindex] = y[newp-xlat*ylat];
190			}
191
192			//point3 = left;
193
194			px = i;
195			py = j - 1;
196			newp = px*ylat + py;
197			if (py >= 0) {
198			interpsum += z[newp];
199			ninterp++;
200			x[newindex] = x[newp];
201			}
202			if (py < 0) {
203			interpsum += z[newp+ylat];
204			ninterp++;
205			x[newindex] = x[newp+ylat];
206			}
207
208			//point4 = right;
209
210			px = i;
211			py = j + 1;
212			newp = px*ylat + py;
213			if (py < ylat) {
214			interpsum += z[newp];
215			ninterp++;
216			x[newindex] = x[newp];
217			}
218			if (py == ylat) {
219			interpsum += z[newp-ylat];
220			ninterp++;
221			x[newindex] = x[newp-ylat];
222			}
223
224			value = interpsum / (double)ninterp;
225
226			z[newindex] = value;
227
228			}
229			}
230			}
231
232			for (i=0; i < xlat; i++) {
233			for (j=0; j < ylat; j++) {
234			newindex = i*ylat + j;
235			c_re(in[newindex]) = z[newindex];
236			c_im(in[newindex]) = 0.0;
237			}
238			}
239
240			fftwnd(p, 1, in, 1, 0, out, 1, 0);
241			for (i=0; i< xlat; i++) {
242			for(j=0; j< ylat; j++) {
243			newindex = i*ylat + j;
244			mag[newindex] += sqrt(pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2));
245			}
246			}
247
248			}
249
250			for (i=0; i< (xlat/2); i++) {
251			for(j=0; j< (ylat/2); j++) {
252			index = i*ylat + j;
253			new_i = i + (xlat/2);
254			new_j = j + (ylat/2);
255			new_index = new_i*ylat + new_j;
256			newmag[new_index] = mag[index];
257			}
258			}
259
260			for (i=(xlat/2); i< xlat; i++) {
261			for(j=0; j< (ylat/2); j++) {
262			index = i*ylat + j;
263			new_i = i - (xlat/2);
264			new_j = j + (ylat/2);
265			new_index = new_i*ylat + new_j;
266			newmag[new_index] = mag[index];
267			}
268			}
269
270			for (i=0; i< (xlat/2); i++) {
271			for(j=(ylat/2); j< ylat; j++) {
272			index = i*ylat + j;
273			new_i = i + (xlat/2);
274			new_j = j - (ylat/2);
275			new_index = new_i*ylat + new_j;
276			newmag[new_index] = mag[index];
277			}
278			}
279
280			for (i=(xlat/2); i< xlat; i++) {
281			for(j=(ylat/2); j< ylat; j++) {
282			index = i*ylat + j;
283			new_i = i - (xlat/2);
284			new_j = j - (ylat/2);
285			new_index = new_i*ylat + new_j;
286			newmag[new_index] = mag[index];
287			}
288			}
289
290			pt1x = x[0];
291			pt1y = y[0];
292			pt2x = x[ylat + 1];
293			pt2y = y[ylat + 1];
294
295			dx = pt2x - pt1x;
296			dy = pt2y - pt1y;
297
298			maxfreqx = 1.0/dx;
299			maxfreqy = 1.0/dy;
300
301
302			maxfreq = sqrt(maxfreqxmaxfreqx + maxfreqymaxfreqy);
303			dfreq = maxfreq/(double)nbins;
304
305			for (i=0; i < nbins; i++) {
306			bin[i] = 0.0;
307			samples[i] = 0;
308			}
309
310			zero_freq_x = xlat/2;
311			zero_freq_y = ylat/2;
312
313			for (i=0; i< xlat; i++) {
314			for(j=0; j< ylat; j++) {
315
316			freq_x = (double)(i - zero_freq_x)*maxfreqx/(double)xlat;
317			freq_y = (double)(j - zero_freq_y)*maxfreqy/(double)ylat;
318
319			freq = sqrt(freq_xfreq_x + freq_yfreq_y);
320
321			whichbin = (int) (freq / dfreq);
322			newindex = i*ylat + j;
323			bin[whichbin] += newmag[newindex];
324			samples[whichbin]++;
325
326			//newindex = i*ylat + j;
327			//printf("%lf\t", newmag[newindex]/(double)nframes);
328			}
329			//printf("\n");
330			}
331
332
333
334			for (i = 0; i < nbins; i++) {
335
336			if (samples[i] > 0) {
337			printf("%lf\t%lf\n", i*dfreq, bin[i]/(double)samples[i]);
338			}
339			}
340
341
342			free(samples);
343			free(bin);
344			free(x);
345			free(y);
346			free(z);
347			free(mag);
348			free(newmag);
349			free(present_in_old);
350			fftwnd_destroy_plan(p);
351			fftw_free(in);
352			fftw_free(out);
353			return(0);
354			}
355
356			/* } */
357			/* for(i=0;i<120;i++) printf("%f\t%f\n", c_re(out[i]), c_im(out[i])); */
358			/* fftw_destroy_plan(p); */
359			/* fftw_free(in); */
360			/* fftw_free(out); */
361			/* averZ = sumZ / (double)nloops; */
362			/* averUx = sumUx / (double)nloops; */
363			/* averUy = sumUy / (double)nloops; */
364			/* averUz = sumUz / (double)nloops; */
365			/* averP = 1.5sumP / (double)nloops-0.5; /
366			/* averTheta = acos(averUz); */
367			/* printf("nloops=%d\n",nloops); */
368			/* printf("sumZ=%f\n",sumZ); */
369			/* printf("average height is : %f\n",averZ); */
370			/* printf("average ux is : %f\n",averUx); */
371			/* printf("average uy is : %f\n",averUy); */
372			/* printf("average uz is : %f\n",averUz); */
373			/* printf("average angle is : %f\n",averTheta); */
374			/* printf("average p is : %f\n",averP); */
375
376			/* fclose(inFile); */
377			/* return 0; */
378			/* } */