Allink  v0.1
ElPolyProfPre.cpp
1 #include "ElPoly.h"
2 
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 }
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 }
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 }
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 }
int PressTrace()
Trace of the pressure profile.
double pVel(int p, int d)
Return the velocity.
int Tens2dCartRad()
Change the pressure profile from cartesian to radial.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
int PressRadial()
Contour plot around the inclusion.
double pCm(int d)
Center of mass of the system.
Definition: VarData.h:924
int CLat2
lateral coordinate
Definition: VarData.h:1078
int Diff2Files(int NBin, int How)
Density and thickness profile arond the nanoparticle.
void SumTens()
Sum more tension profile files.
double PerpTo3(Vettore *Pos, Vettore *Axis)
The vector perpendicolar in three dimension (faster)
bool Open(char *InFile, int BF)
Open the.
Definition: VarData.cpp:105
void ConvoluteMatrix(double *Plot, int NGrid, int NDim, int IfMinImConv)
Convolute with a matrix.
double Pos[3]
Position.
Definition: VarData.h:427
int SurfTens(int NBin)
Radial summation of the 3d tension profile.
double pNanoPos(int n, int d)
Return back folded nano position.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
void RestPress(int NBin)
Pressure difference between the virial and the ideal gas term.
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
void PrintDens(FILE *FileToWrite, double **Plot, double *LatDim, int NBin)
Print the density profile in the surfaces representation.
Matrice computes the algebric operations on matrices.
char * ChooseDraw(int ExtWhat2Draw)
Convert the internal definition for the menu of ElPoly in string.
Definition: ElPolyEl.cpp:279
double pPos(int p, int d)
Return back folded position.
int NFile[2]
First and last file of the list.
Definition: ElPoly.h:418
int CLat1
lateral coordinate
Definition: VarData.h:1076
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int Typ
Type.
Definition: VarData.h:226
void Print()
Print the entries.
void Set(double Val, int Col)
Set the N column.
void Processing(int f)
Information on the current file elaborated.
Definition: ElPolyEl.cpp:369
void VolumeCircSlab(double *VolContr, int NSample)
Fill an array of.
void SetEdge(double Val, int d)
Set Edge.
Definition: VarData.h:976
int CNorm
Normal coordinate.
Definition: VarData.h:1074
int pNPart()
Number of particle.