Allink  v0.1
ElPolyProfPre.cpp
1 #include "ElPoly.h"
2 
3 void ElPoly::SumTens(){
4  int NType = 3;
5  double NumDiff = 0.001;
6  double **Plot = (double **)calloc(NType,sizeof(double));
7  int v[3];
8  double FNorma = 1./(double)(NFile[1]-NFile[0]);
9  char FileName[256];
10  for(int t=0;t<3;t++) Plot[t] = (double *)calloc(CUB(NEdge),sizeof(double));
11  double LatLim[3] = {pEdge(0),pEdge(1),pEdge(2)};
12  double InvLatLim[3] = {1./pEdge(0),1./pEdge(1),1./pEdge(2)};
13  for(int f=NFile[0];f<NFile[1];f++){
14  Processing(f);
15  if(Open(cFile[f],BF_NO))return;
16  //ShiftSys(SHIFT_CM);
17  for(int p=0;p<pNPart();p++){
18  for(int d=0;d<3;d++){
19  v[d] = (int)((pPos(p,d)+NumDiff)*InvLatLim[d]*NEdge);
20  if(v[d] < 0 || v[d] >=NEdge) continue;
21  }
22  int vTot = v[0] + NEdge*(v[1] + NEdge*v[2]);
23  for(int t=0;t<NType;t++)
24  Plot[t][vTot] += pVel(p,t);
25  }
26  }
27  sprintf(FileName,"Av%s",cFile[NFile[0]]);
28  FILE *FileToWrite = fopen(FileName,"w");
29  if(FileToWrite == NULL){printf("Can't open %s\n",FileName);return;}
30  fprintf(FileToWrite,"# l(%.1f %.1f %.1f) r(%.2f %.2f %.2f) v[%d] d[color]\n",LatLim[0],LatLim[1],LatLim[2],Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],NEdge);
31  double NEdgeInv = 1./(double)NEdge;
32  Matrice Mask(5,5);
33  Mask.FillGaussian(.5,3.);
34  Mask.Print();
35  int NDim = 2;
36  int IfMinImConv = 0;
37  for(int t=0;t<NType;t++){
38  Mask.ConvoluteMatrix(Plot[t],NEdge,NDim,IfMinImConv);
39  Mask.ConvoluteMatrix(Plot[t],NEdge,NDim,IfMinImConv);
40  }
41  for(int vx=0;vx<NEdge;vx++)
42  for(int vy=0;vy<NEdge;vy++)
43  for(int vz=0;vz<NEdge;vz++){
44  int vTot = vx + NEdge*(vy + NEdge*vz);
45  double x = vx*LatLim[0]*NEdgeInv;
46  double y = vy*LatLim[1]*NEdgeInv;
47  double z = vz*LatLim[2]*NEdgeInv;
48  if(fabs(Plot[0][vTot] + Plot[1][vTot] + Plot[2][vTot]) > 0.){
49  fprintf(FileToWrite,"{x(%.3f %.3f %.3f) v(%lf %.2f %.2f)}\n",
50  x,y,z,
51  Plot[0][vTot]*FNorma,Plot[1][vTot]*FNorma,Plot[2][vTot]*FNorma);
52  }
53  }
54  for(int t=0;t<NType;t++)
55  free(Plot[t]);
56  free(Plot);
57  fclose(FileToWrite);
58 }
59 int ElPoly::PressRadial(){
60  int NZed = NEdge;
61  int NRad = NEdge;
62  int NType = 3;
63  double *Plot = (double *)calloc(NZed*NRad*NType,sizeof(double));
64  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
65  //Vettore NanoAxis(0,0,1);
66  //Vettore PosRel(3);
67  //Vettore Dist(3);
68  for(int p=0;p<pNPart();p++){
69  double Rad = 0.;
70  for(int d=0;d<2;d++){
71  Rad += SQR(remainder(Nano->Pos[d] - pPos(p,d) - 0.001,pEdge(d)));
72  //PosRel.Set(remainder(Nano->PosBf[d] - pPos(p,d),pEdge(d)),d);
73  }
74  //double RadDist = ABS(Dist.PerpTo3(&PosRel,&NanoAxis));
75  double RadDist = sqrt(Rad);
76  int vr = (int)(RadDist/pEdge(3)*NRad);
77  int vz = (int)((pPos(p,CNorm)+.01)/pEdge(CNorm)*NZed);
78  if(vr < 0 || vr >= NRad) continue;
79  if(vz < 0 || vz >= NZed) continue;
80  for(int t=0;t<3;t++){
81  Plot[(vr*NZed+vz)*NType+t] += pVel(p,t);
82  }
83  }
84  //------------Normalize----------------------
85  double *VolContr = (double *)calloc(NRad,sizeof(double));
86  VolumeCircSlab(VolContr,NRad);
87  double Bound[3];
88  double InvNZed = 1./(double)NZed;
89  for(int t=0,n=0;t<NType;t++){
90  Bound[t] = 0.;
91  for(int vz=0;vz<NZed;vz++){
92  for(int vr=0;vr<NRad;vr++){
93  Plot[(vr*NZed+vz)*NType+t] /= VolContr[vr]*3.;
94  if(Bound[t] < Plot[(vr*NZed+vz)*NType+t])
95  Bound[t] = Plot[(vr*NZed+vz)*NType+t];
96  }
97  }
98  if(Bound[t] < 0.) Bound[t] = 1.;
99  }
100  free(VolContr);
101  //-------------------Write---------------------------
102  FILE *PRadial = fopen("ContourPress.xvl","w");
103  fprintf(PRadial,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(3),pEdge(CNorm),1.,NZed,ChooseDraw(EL_QUAD1));
104  FILE *PNormal = fopen("PressNormal.dat","w");
105  for(int vr=0;vr<NRad;vr++){
106  double PRad[3] = {0.,0.,0.};
107  double r = vr*pEdge(3)/(double)NRad;
108  for(int vz=0;vz<NZed;vz++){
109  //if(Plot[(v*NBin+vv)*NType+0] + Plot[(v*NBin+vv)*NType+1] + Plot[(v*NBin+vv)*NType+2] + Plot[(v*NBin+vv)*NType+3] < .1) continue;
110  double z = vz*pEdge(CNorm)/(double)(NZed);
111  double Press = Plot[(vr*NZed+vz)*NType+0];
112  double Phob = Plot[(vr*NZed+vz)*NType+1];
113  double Phil = Plot[(vr*NZed+vz)*NType+2];
114  if(ABS( Press + Phob + Phil) > 0.)
115  fprintf(PRadial,"{x(%lf %lf 0.) v(%lf %lf %lf)}\n",r,z,Press,Phob,Phil);
116  PRad[0] += Press;PRad[1] += Phob;PRad[2] += Phil;
117  }
118  fprintf(PNormal,"%lf %lf %lf %lf\n",r,PRad[0]*InvNZed,PRad[1]*InvNZed,PRad[2]*InvNZed);
119  }
120  fclose(PRadial);
121  fclose(PNormal);
122  free(Plot);
123  free(VolContr);
124 }
125 int ElPoly::Tens2dCartRad(){
126  int NType = 3;
127  // Press, phob, phil
128  double *TensRad = (double *)calloc(NEdge*NEdge*NType,sizeof(double));
129  double *TensNorm = (double *)calloc(NEdge*NEdge*NType,sizeof(double));
130  double *TensAng = (double *)calloc(NEdge*NEdge*NType,sizeof(double));
131  //x y z
132  double *TensRadTemp = (double *)calloc(NEdge*NEdge*3,sizeof(double));
133  if(NFile[1]-NFile[0] != 3) return 1;
134  for(int f=NFile[0];f<NFile[1];f++){
135  if(Open(cFile[f],BF_PART))return 1;
136  for(int p=0;p<pNPart();p++){
137  int vr = (int)((pPos(p,0)+.0001)*NEdge*pInvEdge(0));
138  int vz = (int)((pPos(p,1)+.0001)*NEdge*pInvEdge(1));
139  int vTot = (vr*NEdge+vz)*NType;
140  //pressure
141  TensRadTemp[vTot+f] += pVel(p,0);
142  //phob density
143  TensRad[vTot+1] += pVel(p,1)/3.;
144  //phil density
145  TensRad[vTot+2] += pVel(p,2)/3.;
146  }
147  }
148  for(int vr=0;vr<NEdge;vr++){
149  for(int vz=0;vz<NEdge;vz++){
150  int vTot = (vr*NEdge+vz)*NType;
151  double Rad = sqrt(SQR(TensRadTemp[vTot+0])+SQR(TensRadTemp[vTot+1]));
152  double Trace = TensRadTemp[vTot+0] + TensRadTemp[vTot+1] + TensRadTemp[vTot+2];
153  double Ang = .5*(Trace - Rad);
154  //double Ang = atan2(TensRadTemp[vTot+1],TensRadTemp[vTot+0]);
155  TensRad[vTot+0] = Rad - .5*(Ang+TensRadTemp[vTot+2]);
156  //TensNorm[vTot] = - TensRadTemp[vTot+2] + .5*(Ang+Rad);
157  TensNorm[vTot] = TensRadTemp[vTot+2] - .5*(TensRadTemp[vTot+0]+TensRadTemp[vTot+1]);
158  TensAng[vTot+0] = Ang - .5*(Rad+TensRadTemp[vTot+2]);
159  }
160  }
161  FILE *RadSurfTens = fopen("SurfTensRad.xvl","w");
162  FILE *NormSurfTens = fopen("SurfTensNorm.xvl","w");
163  FILE *AngSurfTens = fopen("SurfTensAng.xvl","w");
164  fprintf(RadSurfTens,"# l(%.2f %.2f %.2f) r(%.2f %.2f %.2f) v[%d] d[color]\n",pEdge(0),pEdge(1),pEdge(2),Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],NEdge);
165  fprintf(NormSurfTens,"# l(%.2f %.2f %.2f) r(%.2f %.2f %.2f) v[%d] d[color]\n",pEdge(0),pEdge(1),pEdge(2),Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],NEdge);
166  fprintf(AngSurfTens,"# l(%.2f %.2f %.2f) r(%.2f %.2f %.2f) v[%d] d[color]\n",pEdge(0),pEdge(1),pEdge(2),Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],NEdge);
167  for(int vr=0;vr<NEdge;vr++){
168  for(int vz=0;vz<NEdge;vz++){
169  int vTot = (vr*NEdge+vz)*NType;
170  double r = vr*pEdge(0)/NEdge;
171  double z = vz*pEdge(1)/NEdge;
172  if(fabs(TensRad[vTot+0]+TensRad[vTot+1]+TensRad[vTot+2])>0.){
173  fprintf(RadSurfTens,"{x(%.4f %.4f %.4f) v(%lf %lf %lf)}\n",
174  r,z,0.,TensRad[vTot+0],TensRad[vTot+1],TensRad[vTot+2]);
175  fprintf(NormSurfTens,"{x(%.4f %.4f %.4f) v(%lf %lf %lf)}\n",
176  r,z,0.,TensNorm[vTot],TensRad[vTot+1],TensRad[vTot+2]);
177  fprintf(AngSurfTens,"{x(%.4f %.4f %.4f) v(%lf %lf %lf)}\n",
178  r,z,0.,TensAng[vTot],TensRad[vTot+1],TensRad[vTot+2]);
179  }
180  }
181  }
182  free(TensRad);
183  free(TensNorm);
184  free(TensAng);
185  free(TensRadTemp);
186  fclose(RadSurfTens);
187  fclose(NormSurfTens);
188  fclose(AngSurfTens);
189 }
190 int ElPoly::SurfTens(int NBin){
191  int NType = 3;
192  double **TensRad = (double **)calloc(NBin,sizeof(double));
193  double **NormRad = (double **)calloc(NBin,sizeof(double));
194  double **TensAng = (double **)calloc(NBin,sizeof(double));
195  double **NormAng = (double **)calloc(NBin,sizeof(double));
196  double **TensNorm= (double **)calloc(NEdge,sizeof(double));
197  double **NormNorm= (double **)calloc(NEdge,sizeof(double));
198  double **TensCart = (double **)calloc(NBin,sizeof(double));
199  double **NormCart = (double **)calloc(NBin,sizeof(double));
200  for(int v=0;v<NBin;v++){
201  TensRad[v] = (double *)calloc(NType,sizeof(double));
202  NormRad[v] = (double *)calloc(NType,sizeof(double));
203  TensNorm[v] = (double *)calloc(NType,sizeof(double));
204  NormNorm[v] = (double *)calloc(NType,sizeof(double));
205  TensAng[v] = (double *)calloc(NType,sizeof(double));
206  NormAng[v] = (double *)calloc(NType,sizeof(double));
207  TensCart[v] = (double *)calloc(NType,sizeof(double));
208  NormCart[v] = (double *)calloc(NType,sizeof(double));
209  }
210  //double *Plot = (double *)calloc(NBin*NBin,sizeof(double));
211  SetEdge(MIN((.5*pEdge(CLat1)),(.5*pEdge(CLat2))),3);
212  double *VolContr = (double *)calloc(NBin,sizeof(double));
213  VolumeCircSlab(VolContr,NBin);
214  double Prex = 0.;
215  double Prey = 0.;
216  int IfTens = 0;
217  if(NFile[1]-NFile[0] == 3) IfTens = 1;
218  for(int f=NFile[0];f<NFile[1];f++){
219  Processing(f);
220  if(Open(cFile[f],BF_PART))return 0;
221  //BackFold(BF_PART);
222  //Vettore NanoAxis(Nano->Axis[0],Nano->Axis[1],Nano->Axis[2]);
223  Vettore NanoAxis(0,0,1);
224  Vettore PosRel(3);
225  Vettore Dist(3);
226  for(int p=0;p<pNPart();p++){
227  for(int d=0;d<3;d++){
228  PosRel.Set(remainder(pNanoPos(0,d) - pPos(p,d),pEdge(d)),d);
229  }
230  double RadDist = ABS(Dist.PerpTo3(&PosRel,&NanoAxis));
231  int vr = (int)(RadDist/pEdge(3)*NBin);
232  int vz = (int)((pPos(p,CNorm)+.01)/pEdge(CNorm)*NEdge);
233  if(vr < 0 || vr >= NBin){continue;}
234  if(vz < 0 || vz >= NEdge){printf("%d out of %d\n",vz,NBin);continue;}
235  for(int t=1;t<3;t++){
236  TensRad[vr][t] += pVel(p,t);
237  NormRad[vr][t] += 1.;
238  TensNorm[vz][t] += pVel(p,t);
239  NormNorm[vz][t] += 1.;
240  }
241  if(IfTens == 0){
242  TensRad[vr][0] += pVel(p,0);
243  NormRad[vr][0] += 1.;
244  TensNorm[vz][0] += pVel(p,0);
245  NormNorm[vz][0] += 1.;
246  }
247  else if(IfTens){
248  TensCart[vr][f] += pVel(p,0);
249  NormCart[vr][f] += 1.;
250  double vPre = (pVel(p,1)+pVel(p,2)-pVel(p,0));
251  if(f==CLat1 || f == CLat2){
252  TensNorm[vz][0] -= .5*vPre;
253  if(NType>3){
254  TensNorm[vz][3] -= .5*vPre*pPos(p,CNorm);
255  NormNorm[vz][3] += 1.;
256  }
257  if(NType>4){
258  TensNorm[vz][4] -= .5*vPre*pPos(p,CNorm)*pPos(p,CNorm);
259  NormNorm[vz][4] += 1.;
260  }
261  if(f==CLat1)
262  Prex += vPre;
263  else
264  Prey += vPre;
265  }
266  if(f==CNorm){
267  TensNorm[vz][0] += vPre;
268  if(NType>3){
269  TensNorm[vz][3] += vPre*pPos(p,CNorm);
270  NormNorm[vz][3] += 1.;
271  }
272  if(NType>4){
273  TensNorm[vz][4] += vPre*pPos(p,CNorm)*pPos(p,CNorm);
274  NormNorm[vz][4] += 1.;
275  }
276  }
277  NormNorm[vz][0] += 1.;
278  }
279  }
280  }
281  printf("\n");
282  double VolElm = 1.;//pEdge(CNorm)/(double)NBin;
283  FILE *File2Write = fopen("TensRadial.dat","w");
284  for(int v=0;v<NBin;v++){
285  for(int d=0;d<3;d++){
286  TensCart[v][d] /= NormCart[v][d] > 0. ? VolContr[v]*NormCart[v][d] : 1.;
287  }
288  TensRad[v][0] = sqrt( SQR(TensCart[v][CLat1])+SQR(TensCart[v][CLat2]) );
289  TensRad[v][0] -= 5.*.5*(TensCart[v][0]+TensCart[v][1]+TensCart[v][2]-sqrt( SQR(TensCart[v][CLat1])+SQR(TensCart[v][CLat2])) );
290  //TensRad[v][0] += atan2(TensCart[v][CLat2],TensCart[v][CLat1]);
291  TensRad[v][0] -= .5*TensCart[v][CNorm];
292  }
293  for(int v=0;v<NBin;v++){
294  fprintf(File2Write,"%lf ",v/(double)NBin*pEdge(3));
295  for(int t=0;t<NType;t++){
296  TensRad[v][t] /= NormRad[v][t] > 0. ? NormRad[v][t]: 1.;
297  fprintf(File2Write,"%lf ",TensRad[v][t]);
298  }
299  fprintf(File2Write,"\n");
300  }
301  fclose(File2Write);
302  File2Write = fopen("TensNormal.dat","w");
303  fprintf(File2Write,"#Press DensPhob DensPhil SponCurv SaddleSplay \n");
304  for(int v=0;v<NEdge;v++){
305  fprintf(File2Write,"%lf ",v/(double)NEdge*pEdge(CNorm));
306  for(int t=0;t<NType;t++){
307  TensNorm[v][t] /= NormNorm[v][t] > 0. ? VolElm*NormNorm[v][t] : 1.;
308  fprintf(File2Write,"%lf ",TensNorm[v][t]);
309  }
310  fprintf(File2Write,"\n");
311  }
312  fclose(File2Write);
313  File2Write = fopen("TensAngle.dat","w");
314  fprintf(File2Write,"#Press DensPhob DensPhil SponCurv SaddleSplay \n");
315  for(int v=0;v<NEdge;v++){
316  fprintf(File2Write,"%lf ",v/(double)NEdge*pEdge(CNorm));
317  for(int t=0;t<NType;t++){
318  TensAng[v][t] /= NormNorm[v][t] > 0. ? VolElm*NormNorm[v][t] : 1.;
319  fprintf(File2Write,"%lf ",TensNorm[v][t]);
320  }
321  fprintf(File2Write,"\n");
322  }
323  fclose(File2Write);
324  for(int v=0;v<NBin;v++){
325  free(TensRad[v]);
326  free(NormRad[v]);
327  }
328  for(int v=0;v<NEdge;v++){
329  free(TensNorm[v]);
330  free(NormNorm[v]);
331  }
332  free(TensRad);
333  free(NormRad);
334  free(TensNorm);
335  free(NormNorm);
336  free(VolContr);
337  //free(Plot);
338 }
339 int ElPoly::PressTrace(){
340  int NType = 3;
341  double *Sum = (double *)calloc(NType*CUBE(NEdge),sizeof(double));
342  if(NFile[1]-NFile[0] != 3){
343  printf("Just three files\n");
344  return 1;
345  }
346  for(int f=NFile[0];f<NFile[1];f++){
347  Processing(f);
348  if(Open(cFile[f],BF_NO))return 1;
349  for(int p=0;p<pNPart();p++){
350  int vx = (int)((pPos(p,CLat1)+.1)/pEdge(CLat1)*NEdge);
351  int vy = (int)((pPos(p,CLat2)+.1)/pEdge(CLat2)*NEdge);
352  int vz = (int)((pPos(p,CNorm)+.1)/pEdge(CNorm)*NEdge);
353  if(vx < 0 || vx >= NEdge) continue;
354  if(vy < 0 || vy >= NEdge) continue;
355  if(vz < 0 || vz >= NEdge) continue;
356  int vTot = (vx*NEdge+vy)*NEdge+vz;
357  Sum[vTot*NType+0] += pVel(p,0)/3.;
358  Sum[vTot*NType+1] += pVel(p,1)/3.;
359  Sum[vTot*NType+2] += pVel(p,2)/3.;
360  }
361  }
362  printf("\n");
363  FILE *FileToWrite = fopen("PressTrace.xvl","w");
364  fprintf(FileToWrite,"# l(%.1f %.1f %.1f) r(%.2f %.2f %.2f) v[%d] d[color]\n",pEdge(0),pEdge(1),pEdge(2),Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],NEdge);
365  for(int vx=0;vx<NEdge;vx++){
366  double x = vx*pEdge(CLat1)/(double)NEdge;
367  for(int vy=0;vy<NEdge;vy++){
368  double y = vy*pEdge(CLat2)/(double)NEdge;
369  for(int vz=0;vz<NEdge;vz++){
370  double z = vz*pEdge(CNorm)/(double)NEdge;
371  int v = (vx*NEdge+vy)*NEdge+vz;
372  if(Sum[v*NType+0] + Sum[v*NType+1] + Sum[v*NType+2] < 0.1) continue;
373  fprintf(FileToWrite,"{x(%.2f %.2f %.2f)",x,y,z);
374  fprintf(FileToWrite," v( %lf %.2f %.2f)}\n",Sum[v*NType+0],Sum[v*NType+1],Sum[v*NType+2]);
375  }
376  }
377  }
378  fclose(FileToWrite);
379  free(Sum);
380  //free(Plot);
381 }
382 void ElPoly::RestPress(int NBin){
383  int NType = 3;
384  double Round = 0.00001;
385  double InvNBin = 1./(double)NBin;
386  double **Plot1 = (double **)calloc(NType,sizeof(double));
387  for(int t=0;t<NType;t++){
388  Plot1[t] = (double *)calloc(NBin*NBin,sizeof(double));
389  }
390  double *Count1 = (double *)calloc(NType*NBin*NBin,sizeof(double));
391  double *Count2 = (double *)calloc(NType*NBin*NBin,sizeof(double));
392  double LatDim[3] = {pEdge(CLat1),pEdge(CLat2),pEdge(CNorm)};
393  double InvLatDim[2] = {1./LatDim[0],1./LatDim[1]};
394  double FirstCm[3] = {pCm(CLat1),pCm(CLat2),pCm(CNorm)};
395  for(int p=0;p<pNPart();p++){
396  int vr = (int)((pPos(p,CLat1)+Round)*InvLatDim[0]*NBin);
397  if (vr < 0 || vr >= NBin) continue;
398  int vz = (int)((pPos(p,CLat2)+Round)*InvLatDim[1]*NBin);
399  if (vz < 0 || vz >= NBin) continue;
400  int t = Pm[p].Typ;
401  for(int t=0;t<3;t++){
402  Plot1[t][(vr*NBin+vz)] += pVel(p,t);//pPos(p,2);
403  Count1[(vr*NBin+vz)*NType+t] += 1.;
404  }
405  }
406  FILE *FTecPlot = fopen("TecPlotPressDiff.dat","w");
407  fprintf(FTecPlot,"VARIABLES = \"R\", \"Z\", \"diff\", \"d1\",\"d2\"\n");
408  fprintf(FTecPlot,"ZONE J=%d, K=%d, F=POINT\n",NBin,NBin);
409  for(int vx=0;vx<NBin;vx++){
410  for(int vy=0;vy<NBin;vy++){
411  double r = vx*pEdge(0)/(double)NBin;
412  double z = vy*pEdge(1)/(double)NBin - pEdge(1)*.5;
413  //fprintf(FTecPlot,"%lf %lf %lf %lf %lf\n",r,z,diff0,diff1,diff2);
414  double Diff = Plot1[0][vx*NBin+vy] - Plot1[1][vx*NBin+vy] - Plot1[2][vx*NBin+vy];
415  fprintf(FTecPlot,"%lf %lf %lf %lf %lf\n",r,z,Diff,Plot1[1][vx*NBin+vy],Plot1[2][vx*NBin+vy]);
416  }
417  }
418  fclose(FTecPlot);
419  //difference
420  FILE *Difference = fopen("PressDifference.xvl","w");
421  fprintf(Difference,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NBin,ChooseDraw(EL_COLOR));
422  int link[4] = {0,0,0,0};
423  for(int t=0,p=0,c=0;t<NType;t++){
424  for(int vx=0;vx<NBin-1;vx++){
425  for(int vy=0;vy<NBin-1;vy++){
426  double Diff = Plot1[0][vx*NBin+vy] - Plot1[1][vx*NBin+vy] - Plot1[2][vx*NBin+vy];
427  fprintf(Difference,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf)}\n",p,c,t,vx*InvNBin*pEdge(0),vy*InvNBin*pEdge(1),0.,Diff,Plot1[1][vx*NBin+vy],Plot1[2][vx*NBin+vy]);
428  p++;
429  }
430  }
431  }
432  fclose(Difference);
433  for(int t=0;t<NType;t++){
434  free(Plot1[t]);
435  }
436  free(Plot1);
437  free(Count1);
438 }
439 int ElPoly::Diff2Files(int NBin,int How){
440  if(NFile[1] - NFile[0] != 2){
441  printf("The number of files must be two\n");
442  return 1;
443  }
444  int NType = 3;
445  double Round = 0.00001;
446  NBin = MIN(NBin,NEdge);
447  double InvNBin = 1./(double)NBin;
448  double **Plot1 = (double **)calloc(NType,sizeof(double));
449  double **Plot2 = (double **)calloc(NType,sizeof(double));
450  double **PlotDiff = (double **)calloc(NType,sizeof(double));
451  for(int t=0;t<NType;t++){
452  Plot1[t] = (double *)calloc(NBin*NBin,sizeof(double));
453  Plot2[t] = (double *)calloc(NBin*NBin,sizeof(double));
454  PlotDiff[t] = (double *)calloc(NBin*NBin,sizeof(double));
455  }
456  double *Count1 = (double *)calloc(NType*NBin*NBin,sizeof(double));
457  double *Count2 = (double *)calloc(NType*NBin*NBin,sizeof(double));
458  double LatDim[3] = {pEdge(CLat1),pEdge(CLat2),pEdge(CNorm)};
459  double InvLatDim[2] = {1./LatDim[0],1./LatDim[1]};
460  double FirstCm[3] = {pCm(CLat1),pCm(CLat2),pCm(CNorm)};
461  for(int p=0;p<pNPart();p++){
462  int vr = (int)((pPos(p,CLat1)+Round)*InvLatDim[0]*NBin);
463  if(vr < 0 || vr >= NBin) continue;
464  int vz = (int)((pPos(p,CLat2)+Round)*InvLatDim[1]*NBin);
465  if(vz < 0 || vz >= NBin) continue;
466  if(How==0){//Press
467  for(int t=0;t<3;t++){
468  Plot1[t][vr*NBin+vz] += pVel(p,t);
469  Count1[(vr*NBin+vz)*NType+t] += 1.;
470  }
471  }
472  else{//Dens
473  int t = Pm[p].Typ;
474  Plot1[t][vr*NBin+vz] += pPos(p,2);
475  Count1[(vr*NBin+vz)*NType+t] += 1.;
476  }
477  }
478  if(Open(cFile[NFile[1]-1],BF_NO))return 1;
479  for(int p=0;p<pNPart();p++){
480  int vr = (int)((pPos(p,CLat1)+Round+1.)*InvLatDim[0]*NBin);
481  if (vr < 0 || vr >= NBin) continue;
482  int vz = (int)((pPos(p,CLat2)+Round)*InvLatDim[1]*NBin);
483  if (vz < 0 || vz >= NBin) continue;
484  if(How==0){//Press
485  for(int t=0;t<3;t++){
486  Plot2[t][vr*NBin+vz] += pVel(p,t);
487  Count2[(vr*NBin+vz)*NType+t] += 1.;
488  }
489  }
490  else{//Dens
491  int t = Pm[p].Typ;
492  Plot2[t][vr*NBin+vz] += pPos(p,2);
493  Count2[(vr*NBin+vz)*NType+t] += 1.;
494  }
495  }
496  Matrice Mask(5,5);
497  Mask.FillGaussian(.5,3.);
498  Mask.Print();
499  for(int t=0;t<NType;t++){
500  //ConvoluteMatrix(Plot1[t],NBin,&Mask,2);
501  // ConvoluteMatrix(Plot1[t],NBin,&Mask,2);
502  //ConvoluteMatrix(Plot2[t],NBin,&Mask,2);
503  // ConvoluteMatrix(Plot2[t],NBin,&Mask,2);
504  }
505  for(int vr=0;vr<NBin;vr++){
506  for(int vz=0;vz<NBin;vz++){
507  for(int t=0;t<NType;t++){
508  Plot1[t][(vr*NBin+vz)] /= Count1[(vr*NBin+vz)*NType+t] > 0. ? Count1[(vr*NBin+vz)*NType+t] : 1.;
509  Plot2[t][(vr*NBin+vz)] /= Count2[(vr*NBin+vz)*NType+t] > 0. ? Count2[(vr*NBin+vz)*NType+t] : 1.;
510  PlotDiff[t][(vr*NBin+vz)] = Plot1[t][(vr*NBin+vz)] - Plot2[t][(vr*NBin+vz)];
511  }
512  int t = 0;
513  }
514  }
515  int NDim = 2;
516  int IfMinImConv = 0;
517  for(int t=0;t<NType;t++){
518  Mask.ConvoluteMatrix(PlotDiff[t],NBin,NDim,IfMinImConv);
519  // Mask.ConvoluteMatrix(PlotDiff[t],NBin,2);
520  }
521  //tecplot
522  FILE *FTecPlot = fopen("TecPlotDiff.dat","w");
523  fprintf(FTecPlot,"VARIABLES = \"R\", \"Z\", \"diff\", \"d1\",\"d2\"\n");
524  fprintf(FTecPlot,"ZONE J=%d, K=%d, F=POINT\n",NBin,NBin);
525  //difference
526  FILE *Difference;
527  if(How) Difference = fopen("DensDifference.rzd","w");
528  else Difference = fopen("PressDifference.xvl","w");
529  if(How==1){//Dens
530  PrintDens(Difference,PlotDiff,LatDim,NBin);
531  for(int vx=0;vx<NBin;vx++){
532  for(int vy=0;vy<NBin;vy++){
533  double r = vx*LatDim[0]/(double)NBin;
534  double z = vy*LatDim[1]/(double)NBin - LatDim[1]*.5;
535  fprintf(FTecPlot,"%lf %lf %lf %lf %lf\n",r,z,PlotDiff[0][vx*NBin+vy],Plot1[0][vx*NBin+vy],Plot2[0][vx*NBin+vy]);
536  }
537  }
538  }
539  else{//Pre
540  fprintf(Difference,"# l(%.1f %.1f %.1f) v[%d] d[color]\n",LatDim[0],LatDim[1],LatDim[2],NBin);
541  for(int vr=0;vr<NBin;vr++){
542  for(int vz=0;vz<NBin;vz++){
543  double NanoAdded = Plot1[2][(vr*NBin+vz)];
544  double Phob = PlotDiff[0][(vr*NBin+vz)];
545  double Phil = PlotDiff[1][(vr*NBin+vz)];
546  double r = (vr)*InvNBin*LatDim[0];
547  double z = (vz)*InvNBin*LatDim[1];
548  double dens = (PlotDiff[2][(vr*NBin+vz)]);
549  fprintf(Difference,"{x(%lf %lf %lf) v(%lf %lf %lf)}\n",r,z,dens,NanoAdded,Phob,Phil);
550  double z1 = vz*LatDim[1]/(double)NBin - LatDim[1]*.5;
551  fprintf(FTecPlot,"%lf %lf %lf %lf %lf\n",r,z1,PlotDiff[0][vr*NBin+vz],Plot1[1][vr*NBin+vz],Plot2[1][vr*NBin+vz]);
552  }
553  }
554  }
555  fclose(FTecPlot);
556  fclose(Difference);
557  for(int t=0;t<NType;t++){
558  free(Plot1[t]);
559  free(Plot2[t]);
560  }
561  free(Plot1);
562  free(Plot2);
563  free(Count1);
564  free(Count2);
565  return 0;
566 }
ElPoly::PressTrace
int PressTrace()
Trace of the pressure profile.
Definition: ElPolyProfPre.cpp:339
VarData::pVel
double pVel(int p, int d)
Return the velocity.
Definition: VarDataComm.cpp:202
ElPoly::Tens2dCartRad
int Tens2dCartRad()
Change the pressure profile from cartesian to radial.
Definition: ElPolyProfPre.cpp:125
VarData::Nano
NANO * Nano
Extra particle.
Definition: VarData.h:1044
Vettore
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
VarData::pInvEdge
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
ElPoly::PressRadial
int PressRadial()
Contour plot around the inclusion.
Definition: ElPolyProfPre.cpp:59
VarData::pCm
double pCm(int d)
Center of mass of the system.
Definition: VarData.h:924
VarData::CLat2
int CLat2
lateral coordinate
Definition: VarData.h:1078
ElPoly::Diff2Files
int Diff2Files(int NBin, int How)
Density and thickness profile arond the nanoparticle.
Definition: ElPolyProfPre.cpp:439
ElPoly::SumTens
void SumTens()
Sum more tension profile files.
Definition: ElPolyProfPre.cpp:3
Vettore::PerpTo3
double PerpTo3(Vettore *Pos, Vettore *Axis)
The vector perpendicolar in three dimension (faster)
Definition: MatematicaVect.cpp:415
VarData::Open
bool Open(char *InFile, int BF)
Open the.
Definition: VarData.cpp:105
Matrice::ConvoluteMatrix
void ConvoluteMatrix(double *Plot, int NGrid, int NDim, int IfMinImConv)
Convolute with a matrix.
Definition: MatematicaAlgebra.cpp:937
NANO::Pos
double Pos[3]
Position.
Definition: VarData.h:427
ElPoly::SurfTens
int SurfTens(int NBin)
Radial summation of the 3d tension profile.
Definition: ElPolyProfPre.cpp:190
VarData::pNanoPos
double pNanoPos(int n, int d)
Return back folded nano position.
Definition: VarDataComm.cpp:225
VarData::pEdge
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
ElPoly::RestPress
void RestPress(int NBin)
Pressure difference between the virial and the ideal gas term.
Definition: ElPolyProfPre.cpp:382
VarData::NEdge
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
Matrice::FillGaussian
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
Definition: MatematicaAlgebra.cpp:375
ElPoly::PrintDens
void PrintDens(FILE *FileToWrite, double **Plot, double *LatDim, int NBin)
Print the density profile in the surfaces representation.
Definition: ElPolyProfDens.cpp:414
Matrice
Matrice computes the algebric operations on matrices.
Definition: MatematicaMatrix.h:10
ElPoly::ChooseDraw
char * ChooseDraw(int ExtWhat2Draw)
Convert the internal definition for the menu of ElPoly in string.
Definition: ElPolyEl.cpp:279
VarData::pPos
double pPos(int p, int d)
Return back folded position.
Definition: VarDataComm.cpp:192
ElPoly::NFile
int NFile[2]
First and last file of the list.
Definition: ElPoly.h:418
VarData::CLat1
int CLat1
lateral coordinate
Definition: VarData.h:1076
VarData::Pm
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
PART::Typ
int Typ
Type.
Definition: VarData.h:226
Matrice::Print
void Print()
Print the entries.
Definition: MatematicaAlgebra.cpp:492
Vettore::Set
void Set(double Val, int Col)
Set the N column.
Definition: MatematicaVect.cpp:218
ElPoly::Processing
void Processing(int f)
Information on the current file elaborated.
Definition: ElPolyEl.cpp:369
VarData::VolumeCircSlab
void VolumeCircSlab(double *VolContr, int NSample)
Fill an array of.
Definition: VarDataProfile.cpp:75
VarData::SetEdge
void SetEdge(double Val, int d)
Set Edge.
Definition: VarData.h:976
VarData::CNorm
int CNorm
Normal coordinate.
Definition: VarData.h:1074
VarData::pNPart
int pNPart()
Number of particle.
Definition: VarDataComm.cpp:246
Generated by   doxygen 1.8.11