trunk/dp/read_new.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;

  FFTW_COMPLEX *in, *out;
  fftwnd_plan p;
        

  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);

  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];
            } 

            //point2 = top;

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

            //point3 = left;

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

            } 

            //point4 = right;

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

              
            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];
    }
  }
      
  for (i=0; i< xlat; i++) {
    for(j=0; j< ylat; j++) {
      newindex = i*ylat + j;
      printf("%lf\t", newmag[newindex]/(double)nframes);
    }
    printf("\n");
  }
  
  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:	581
Committed:	Wed Jul 9 15:14:05 2003 UTC (20 years, 11 months ago) by xsun
Content type:	text/plain
File size:	6692 byte(s)
Log Message:	This is the dp and read file
#	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
23			FFTW_COMPLEX in, out;
24			fftwnd_plan p;
25
26
27			inFile = fopen("./dp_file","r");
28			if (inFile == NULL) {
29			printf("Error opening file\n");
30			exit(-1);
31			}
32
33			fgets(s, 500, inFile);
34			nsites = atoi(s);
35			N = 2*nsites;
36			fscanf(inFile,"%s",s);
37			xlat=atoi(s);
38			fscanf(inFile,"%s",s);
39			ylat=atoi(s);
40
41			//sanity check:
42			if (2nsites != xlat ylat) {
43			printf("Lattice Mismatch: nsites = %d, xlat = %d, ylat = %d\n", nsites, xlat, ylat);
44			exit(-1);
45			}
46
47			rewind(inFile);
48
49			in = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
50			out = fftw_malloc(sizeof(FFTW_COMPLEX) * N);
51			p = fftw2d_create_plan(xlat, ylat, FFTW_FORWARD, FFTW_ESTIMATE);
52
53			x = (double ) malloc(sizeof(double) N);
54			y = (double ) malloc(sizeof(double) N);
55			z = (double ) malloc(sizeof(double) N);
56			mag = (double ) malloc(sizeof(double) N);
57			newmag = (double ) malloc(sizeof(double) N);
58			present_in_old = (int ) malloc(sizeof(int) N);
59
60			for (i=0; i< xlat; i++) {
61			for(j=0; j< ylat; j++) {
62			newindex = i*ylat + j;
63			mag[newindex] = 0.0;
64			newmag[newindex] = 0.0;
65			}
66			}
67
68			sumZ = 0.0;
69			sumUx = 0.0;
70			sumUy = 0.0;
71			sumUz = 0.0;
72			sumP = 0.0;
73
74			nloops = 0;
75			nframes = 0;
76
77			while(!feof(inFile))
78			{
79			nframes++;
80			fgets(s, 500, inFile);
81			nsites = atoi(s);
82			fscanf(inFile,"%s",s);
83			xlat=atoi(s);
84			fscanf(inFile,"%s",s);
85			ylat=atoi(s);
86
87			//sanity check:
88			if (2nsites != xlat ylat) {
89			printf("Lattice Mismatch: nsites = %d, xlat = %d, ylat = %d\n", nsites, xlat, ylat);
90			exit(-1);
91			}
92
93			for (i = 0; i < N; i++ ) {
94			present_in_old[i] = 0;
95			}
96
97			for (i=0;i<xlat;i=i+2)
98			{
99			for(j=0;j<ylat;j++)
100			{
101
102			newy = j;
103			if ((j % 2) == 0) {
104			newx = i;
105			} else {
106			newx = i + 1;
107			}
108			newindex = newx*ylat + newy;
109
110			fscanf(inFile,"%s",s);
111			fscanf(inFile,"%s",s);
112			x[newindex]=strtod(s,&endptr);
113			fscanf(inFile,"%s",s);
114			y[newindex]=strtod(s,&endptr);
115			fscanf(inFile,"%s",s);
116			z[newindex] = strtod(s,&endptr);
117			present_in_old[newindex] = 1;
118			// printf("z=%f\t",z[newindex]);
119			fscanf(inFile,"%s",s);
120			fscanf(inFile,"%s",s);
121			ux=strtod(s,&endptr);
122			fscanf(inFile,"%s",s);
123			uy=strtod(s,&endptr);
124			fscanf(inFile,"%s",s);
125			uz=strtod(s,&endptr);
126			// p1 = uz * uz;
127			// sumZ = sumZ + z[i];
128			// sumUx = sumUx + ux;
129			// sumUy = sumUy + uy;
130			// sumUz = sumUz + uz;
131			// sumP = sumP + p1;
132			nloops += 1;
133			fgets(s, 500, inFile);
134			}
135			}
136
137			for (i=0; i< xlat; i++) {
138			for(j=0; j< ylat; j++) {
139			newindex = i*ylat + j;
140			if (present_in_old[newindex] == 0) {
141			// interpolate from surrounding points:
142
143			interpsum = 0.0;
144			ninterp = 0;
145
146			//point1 = bottom;
147
148			px = i - 1;
149			py = j;
150			newp = px*ylat + py;
151			if (px >= 0) {
152			interpsum += z[newp];
153			ninterp++;
154			y[newindex] = y[newp];
155			}
156
157			//point2 = top;
158
159			px = i + 1;
160			py = j;
161			newp = px*ylat + py;
162			if (px < xlat) {
163			interpsum += z[newp];
164			ninterp++;
165			y[newindex] = y[newp];
166			}
167
168			//point3 = left;
169
170			px = i;
171			py = j - 1;
172			newp = px*ylat + py;
173			if (py >= 0) {
174			interpsum += z[newp];
175			ninterp++;
176			x[newindex] = x[newp];
177
178			}
179
180			//point4 = right;
181
182			px = i;
183			py = j + 1;
184			newp = px*ylat + py;
185			if (py < ylat) {
186			interpsum += z[newp];
187			ninterp++;
188			x[newindex] = x[newp];
189			}
190
191
192			value = interpsum / (double)ninterp;
193
194			z[newindex] = value;
195
196			}
197			}
198			}
199
200			for (i=0; i < xlat; i++) {
201			for (j=0; j < ylat; j++) {
202			newindex = i*ylat + j;
203			c_re(in[newindex]) = z[newindex];
204			c_im(in[newindex]) = 0.0;
205			}
206			}
207
208			fftwnd(p, 1, in, 1, 0, out, 1, 0);
209			for (i=0; i< xlat; i++) {
210			for(j=0; j< ylat; j++) {
211			newindex = i*ylat + j;
212			mag[newindex] += sqrt(pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2));
213			}
214			}
215
216			}
217
218			for (i=0; i< (xlat/2); i++) {
219			for(j=0; j< (ylat/2); j++) {
220			index = i*ylat + j;
221			new_i = i + (xlat/2);
222			new_j = j + (ylat/2);
223			new_index = new_i*ylat + new_j;
224			newmag[new_index] = mag[index];
225			}
226			}
227
228			for (i=(xlat/2); i< xlat; i++) {
229			for(j=0; j< (ylat/2); j++) {
230			index = i*ylat + j;
231			new_i = i - (xlat/2);
232			new_j = j + (ylat/2);
233			new_index = new_i*ylat + new_j;
234			newmag[new_index] = mag[index];
235			}
236			}
237
238			for (i=0; i< (xlat/2); i++) {
239			for(j=(ylat/2); j< ylat; j++) {
240			index = i*ylat + j;
241			new_i = i + (xlat/2);
242			new_j = j - (ylat/2);
243			new_index = new_i*ylat + new_j;
244			newmag[new_index] = mag[index];
245			}
246			}
247
248			for (i=(xlat/2); i< xlat; i++) {
249			for(j=(ylat/2); j< ylat; j++) {
250			index = i*ylat + j;
251			new_i = i - (xlat/2);
252			new_j = j - (ylat/2);
253			new_index = new_i*ylat + new_j;
254			newmag[new_index] = mag[index];
255			}
256			}
257
258			for (i=0; i< xlat; i++) {
259			for(j=0; j< ylat; j++) {
260			newindex = i*ylat + j;
261			printf("%lf\t", newmag[newindex]/(double)nframes);
262			}
263			printf("\n");
264			}
265
266			free(x);
267			free(y);
268			free(z);
269			free(mag);
270			free(newmag);
271			free(present_in_old);
272			fftwnd_destroy_plan(p);
273			fftw_free(in);
274			fftw_free(out);
275
276			return(0);
277			}
278
279			/* } */
280			/* for(i=0;i<120;i++) printf("%f\t%f\n", c_re(out[i]), c_im(out[i])); */
281			/* fftw_destroy_plan(p); */
282			/* fftw_free(in); */
283			/* fftw_free(out); */
284			/* averZ = sumZ / (double)nloops; */
285			/* averUx = sumUx / (double)nloops; */
286			/* averUy = sumUy / (double)nloops; */
287			/* averUz = sumUz / (double)nloops; */
288			/* averP = 1.5sumP / (double)nloops-0.5; /
289			/* averTheta = acos(averUz); */
290			/* printf("nloops=%d\n",nloops); */
291			/* printf("sumZ=%f\n",sumZ); */
292			/* printf("average height is : %f\n",averZ); */
293			/* printf("average ux is : %f\n",averUx); */
294			/* printf("average uy is : %f\n",averUy); */
295			/* printf("average uz is : %f\n",averUz); */
296			/* printf("average angle is : %f\n",averTheta); */
297			/* printf("average p is : %f\n",averP); */
298
299			/* fclose(inFile); */
300			/* return 0; */
301			/* } */