Allink  v0.1
ElPolyRepr.cpp
1 #include "ElPoly.h"
2 int ElPoly::ProjectionF(int NBin,int Coord){
3  if(Coord > 4 || Coord <0) return 1;
4  int NType = 5;
5  double *Plot = (double *)calloc(NBin*NBin*NType,sizeof(double));
6  double InvNBin = 1./(double)NBin;
7  double RefPos[3] = {0.,0.,0.};
8  for(int d=0;d<3;d++){
9  RefPos[d] = Nano->Pos[d]-.5*pEdge(d);
10  }
11  if(Coord == 3){
12  RefPos[0]=pCm(0);RefPos[1]=pCm(1);RefPos[2]=pCm(2);
13  }
14  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
15  for(int f=NFile[0];f<NFile[1];f++){
16  Processing(f);
17  OpenRisk(cFile[f],BF_PART);
18  ShiftSys(SHIFT_CM);
19  int NPlot = 0;
20  //---Projects-against-one-coordinate--
21  if(Coord < 3){
22  int coord1 = (Coord+1)%3;
23  int coord2 = (Coord+2)%3;
24  for(int p=0;p<pNPart();p++){
25  double x = pPos(p,coord1) - RefPos[coord1];
26  x -= floor(x*pInvEdge(coord1))*pEdge(coord1);
27  double y = pPos(p,coord2) - RefPos[coord2];
28  y -= floor(y*pInvEdge(coord2))*pEdge(coord2);
29  int v = (int)(NBin*x*pInvEdge(coord1));
30  if( v < 0 || v >= NBin) continue;
31  int vv = (int)(NBin*y*pInvEdge(coord2));
32  if( vv < 0 || vv >= NBin) continue;
33  int t = pType(p);
34  if( t < 0 || t > 3) continue;
35  if( CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED) )
36  Plot[(v*NBin+vv)*NType+3] += 1.;
37  Plot[(v*NBin+vv)*NType+t] += 1.;
38  if(p<pNPart()-1)
39  if(pType(p+1) == 1 && pType(p) == 0)
40  Plot[(v*NBin+vv)*NType+4] += 1.;
41  }
42  }
43  //---Projects-against-the-radial-coordinate--
44  else if(Coord == 3){
45  SetEdge(.5*MAX((pEdge(CLat1)),(pEdge(CLat2))),3);
46  for(int p=0;p<pNPart();p++){
47  double x = pPos(p,CLat1) - RefPos[CLat1];
48  x -= floor(x*pInvEdge(CLat1))*pEdge(CLat1);
49  double y = pPos(p,CLat2) - RefPos[CLat2];
50  y -= floor(y*pInvEdge(CLat2))*pEdge(CLat2);
51  double z = pPos(p,CNorm) - RefPos[CNorm];
52  z -= floor(z*pInvEdge(CNorm))*pEdge(CNorm);
53  double r = sqrt(SQR(x)+SQR(y));
54  int v = (int)(NBin*r*pInvEdge(3));
55  if( v < 0 || v >= NBin) continue;
56  int vv = (int)(NBin*pPos(p,CNorm)/pEdge(CNorm));
57  if( vv < 0 || vv >= NBin) continue;
58  int t = pType(p);
59  if( t < 0 || t > 3) continue;
60  if( CHAIN_IF_TYPE(Ch[pChain(p)].Type,CHAIN_ADDED) )
61  Plot[(v*NBin+vv)*NType+3] += 1.;
62  Plot[(v*NBin+vv)*NType+t] += 1.;
63  if(p<pNPart()-1)
64  if(pType(p+1) == 1 && pType(p) == 0)
65  Plot[(v*NBin+vv)*NType+4] += 1.;
66  }
67  }
68  }
69  printf("\n");
70  //-----writes-the-output-file-------------------
71  FILE *FileToWrite = NULL;
72  FileToWrite = fopen("Projection.xyd","w");
73  fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NBin,ChooseDraw(EL_QUAD1));
74  int coord1 = (Coord+1)%3;
75  int coord2 = (Coord+2)%3;
76  if(Coord == 3){
77  coord1 = 3;
78  coord2 = CNorm;
79  }
80  double Max[NType];
81  for(int t=0;t<NType;t++){
82  Max[t] = Plot[t];
83  for(int v=0;v<NBin;v++)
84  for(int vv=0;vv<NBin;vv++)
85  if(Max[t] < Plot[(v*NBin+vv)*NType+t]) Max[t] = Plot[(v*NBin+vv)*NType+t];
86  Max[t] = Max[t] <= 0. ? 1. : Max[t];
87  }
88  //for(int t=0;t<NType-1;t++){
89  for(int t=0;t<1;t++){
90  for(int v=0;v<NBin;v++){
91  for(int vv=0;vv<NBin;vv++){
92  int p = (v*NBin+vv)*NType+t;
93  int c = 0;
94  if(Plot[p] < .1) continue;
95  double x = (v)*InvNBin*pEdge(CLat1);
96  double y = (vv)*InvNBin*pEdge(CLat2);
97  double dens = Plot[p]/Max[t]*5.+.5*pEdge(CNorm);
98  double NanoAdded = 0.;//Plot[p]/Max[t]+Plot[((v*NBin+vv)*NType+3]/Max[3];
99  double Phob = t == 0 ? Plot[(v*NBin+vv)*NType+0]/Max[0] : 0.;
100  double Phil = t == 1 ? Plot[(v*NBin+vv)*NType+1]/Max[1] : 0.;
101  fprintf(FileToWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf)}\n",p,c,t,x,y,dens,NanoAdded,Phob,Phil);
102  }
103  }
104  }
105  free(Plot);
106  fclose(FileToWrite);
107  return 0;
108 }
109 // int ElPoly::LateralDistr(int NBin){
110 // }
111 int ElPoly::Surface(int NBin,int Threshold){
112  double *PlotAv = (double *)calloc(NBin*NBin,sizeof(double));
113  int NPoint=0;
114  double AreaAv = 0.,SurfAv = 0.;
115  double Max=0.;
116  printf("Edge %lf Threshold %d\n",pEdge(3),Threshold);
117  for(int f=NFile[0];f<NFile[1];f++){
118  Processing(f);
119  if(OpenRisk(cFile[f],BF_PART))return 1;
120  double Surf=0.,Area =0.;
121  double *Plot = (double *)calloc(NBin*NBin,sizeof(double));
122  for(int p=0;p<pNPart();p++){
123  int v = (int)(pPos(p,CLat1)/pEdge(CLat1)*NBin);
124  if( v < 0 || v > NBin){ printf("Oi 0 < %d < %d\n",v,NBin);return 1;}
125  int vv = (int)(pPos(p,CLat2)/pEdge(CLat2)*NBin);
126  if( vv < 0 || vv > NBin){ printf("Oi 0 < %d < %d\n",vv,NBin);return 1;}
127  Plot[v*NBin+vv] += 1.;
128  }
129  for(int v=0;v<NBin;v++)
130  for(int vv=0,n=0;vv<NBin;vv++)
131  PlotAv[v*NBin+vv] += Plot[v*NBin+vv];
132  for(int v=1;v<NBin-1;v++){
133  for(int vv=1,n=0;vv<NBin-1;vv++){
134  if(Plot[v*NBin+vv] > Threshold){
135  NPoint++;
136  //printf("(%d %d) %d %d %d %d\n",v,vv,Plot[v-1][vv],Plot[v+1][vv],Plot[v][vv-1],Plot[v][vv+1]);
137  if(Plot[(v-1)*NBin+vv] > Threshold && v > 0)
138  n++;
139  if(Plot[(v+1)*NBin+vv] > Threshold && v < NBin -1)
140  n++;
141  if(Plot[v*NBin+vv-1] > Threshold && vv > 0)
142  n++;
143  if(Plot[v*NBin+vv+1] > Threshold && vv < NBin -1)
144  n++;
145  if(n == 4){
146  Area += 1.;
147  }
148  else if(n < 4 && n != 0){
149  if(n == 3){
150  Surf += 1.;
151  Area += 1.;
152  }
153  else if(n == 2){
154  Surf += 2;
155  }
156  else if(n == 1){
157  Surf += sqrt(2);
158  }
159  }
160  n = 0;
161  }
162  }
163  }
164  AreaAv += Area;
165  SurfAv += Surf;
166  free(Plot);
167  }
168  printf("\n");
169  //-----------normalize-------------------
170  for(int v=0;v<NBin;v++)
171  for(int vv=0,n=0;vv<NBin;vv++)
172  if(PlotAv[v*NBin+vv] > Max)
173  Max = PlotAv[v*NBin+vv];
174  //------------write-output-----------------
175  FILE *FileToWrite = NULL;
176  FileToWrite = fopen("Surface.xyz","w");
177  double FNorma = 1./(double)(NFile[1] - NFile[0]);
178  double LengthConv = pEdge(0)*pEdge(1)/(double)(NBin*NBin);
179  fprintf(FileToWrite,"#AreaTot %lf Area %lf Surf %lf Threshold %d NChain %d NBin %d\n",pEdge(CLat1)*pEdge(CLat2),AreaAv*FNorma*LengthConv,SurfAv*FNorma*LengthConv,Threshold,pNChain(),NBin);
180  for(int v=0;v<NBin;v++)
181  for(int vv=0,n=0;vv<NBin;vv++)
182  if(PlotAv[v*NBin+vv] > 0.)
183  fprintf(FileToWrite,"%lf %lf %lf\n",(double)v/NBin*pEdge(CLat1),(double)vv/NBin*pEdge(CLat2),PlotAv[v*NBin+vv]/Max);
184  fclose(FileToWrite);
185  free(PlotAv);
186  return 0;
187 }
188 int ElPoly::From3To2d(int NSample,double Param){
189  char FName[120];
190  for(int f=NFile[0];f<NFile[1];f++){
191  Processing(f);
192  if(OpenRisk(cFile[f],BF_NANO));
193  sprintf(FName,"ProjOnSurf%05d.dat",f);
194  FILE *FLine = fopen(FName,"w");
195  fprintf(FLine,"# l(%lf %lf %lf) d[part]\n",1.,pEdge(CLat2)*pInvEdge(CLat1),.5);
196  for(int p=0;p<pNPart();p++){
197  if(pType(p) != 0) continue;
198  fprintf(FLine,"{x(%lf %lf %lf) v(%lf %lf %lf)}\n",pPos(p,CLat1)*pInvEdge(CLat1),pPos(p,CLat2)*pInvEdge(CLat1),.5,pVel(p,0),pVel(p,1),pVel(p,2));
199  }
200  fclose(FLine);
201  }
202  printf("\n");
203 }
204 int ElPoly::From2To1d(int Coord){
205  if(Coord != 0 && Coord != 1){
206  printf("Coordinate accepted are 0 or 1 not %d\n",Coord);
207  return 1;
208  }
209  int SecCoord = (Coord+1)%2;
210  double *Line = (double *) calloc(3*NEdge,sizeof(double));
211  double *Count = (double *) calloc(3*NEdge,sizeof(double));
212  for(int p=0;p<pNPart();p++){
213  int vz = (int)((pPos(p,Coord)+0.001)*pInvEdge(Coord)*NEdge);
214  if(vz < 0 || vz >= NEdge) continue;
215  for(int d=0;d<3;d++){
216  Line[vz*3+d] += pVel(p,d);
217  Count[vz*3+d] += 1.;
218  }
219  }
220  for(int vz=0;vz<NEdge;vz++){
221  for(int d=0;d<3;d++){
222  double Norm = Count[vz*3+d] > 0. ? Count[vz*3+d] : 1.;
223  //Line[vz*3+d] /= Norm;
224  Line[vz*3+d] *= pEdge(SecCoord)/(double)NEdge;
225  }
226  }
227  FILE *FLine = fopen("ProjOnLine.dat","w");
228  for(int v=0;v<NEdge;v++){
229  fprintf(FLine,"%lf %lf %lf %lf\n",v*pEdge(Coord)/(double)NEdge,Line[v*3+0],Line[v*3+1],Line[v*3+2]);
230  }
231  fclose(FLine);
232  free(Line);
233  free(Count);
234 }
235 int ElPoly::From3To1d(int Coord){
236 
237 
238 }
239 int ElPoly::RadialShell(int NBin){
240  double Volume=0;//Global constant
241  double Area=0.;
242  double Norm=0.;
243  double **Plot;
244  double Ypsilon = 0.;
245  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
246  Plot = (double **)calloc(NBin,sizeof(double));
247  for(int i=0;i<NBin;i++){
248  *(Plot+i) = (double *)calloc(NBin,sizeof(double));
249  }
250  for(int f=NFile[0];f<NFile[1];f++){
251  Processing(f);
252  if(OpenRisk(cFile[f],BF_PART))return 0;
253  SetEdge(MIN((pEdge(CLat1)*.5),(pEdge(CLat2)*.5)),3);
254  for(int p=0;p<pNPart();p++){
255  double Rad = sqrt(SQR((pPos(p,CLat1)-pCm(CLat1)))
256  + SQR((pPos(p,CLat2)-pCm(CLat2))) );
257  int vr = (int)(NBin*Rad*pInvEdge(3));
258  int vz = (int)(NBin*pPos(p,CNorm)*pInvEdge(CNorm));
259  //printf("%lf %lf %d %d \n",Rad,pPos(p,CNorm),v,vv);
260  if( vr < 0 || vr >= NBin) continue;
261  if( vz < 0 || vz >= NBin) continue;
262  Plot[vr][vz] += InvNFile;
263  }
264  Ypsilon += pEdge(2)-pCm(2)-1.;
265  }
266  printf("\n");
267  FILE *FileToWrite = fopen("RadialShell.xye","w");
268  //fprintf(FileToWrite,"# l(%.2f %.2f %.2f) v[60] d[chain]\n",Gen->Edge[0],Gen->Edge[1],Gen->Edge[2]);
269  fprintf(FileToWrite,"%lf %lf %lf\n",0.,0.,0.);
270  fprintf(FileToWrite,"%lf %lf %lf\n",pEdge(0),pEdge(1),1.);
271  double Max=0.;
272  for(int i=0;i<NBin;i++){
273  for(int j=0;j<NBin;j++){
274  if(Plot[i][j] > Max) Max = Plot[i][j];
275  }
276  }
277  int th=0;
278  for(int vz=0;vz<NBin;vz++){
279  //for(int vr=NBin-1;vr>=0;vr--){
280  for(int vr=0;vr<NBin;vr++){
281  if(Plot[vr][vz] > 0.){
282  fprintf(FileToWrite,"%lf %lf %lf\n",(double)vr/NBin*pEdge(3),(double)vz/NBin*pEdge(2),Plot[vr][vz]/Max);
283  Norm += Plot[vr][vz];
284  th++;
285  if(th > 4){
286  th = 0;
287  break;
288  }
289  }
290  }
291  }
292  fclose(FileToWrite);
293  return 0;
294  Mat->Ypsilon = Ypsilon*InvNFile;
295  double Vol = 1.;//pNPart()/(Gen->rho/(double)pNPCh()*CUB(5.));
296  Mat->PreFact = 3./8.*pow((3.*Vol)/DUE_PI,1./3.);
298  int NRadici = 4;
299  printf("Volume %lf Cm %lf Area %lf # to Invert %lf\n",(double)pNPart()/10.,pCm(2),Area,pow(DUE_PI/(2.*pNPart()/10.),.25));
300  double *Radici;
301  Radici = (double *)malloc(NRadici*sizeof(double));
302  int NZeri = Mat->Zeri(0.,DUE_PI/2.,Radici,NRadici);
303  for(int i=0;i<NZeri;i++){
304  Radici[i] *= 360./DUE_PI;
305  printf("Angle %lf\n",Radici[i]);
306  }
307  if(NZeri == 0){
308  printf("The function has no real roots\n");
309  }
310  free(Plot);
311  return 0;
312 }
313 void ElPoly::IsoSurf(int NSample,double *IsoLevel,int NIso){
314  int NPairF = NFile[1]-NFile[0];
315  double OldPos[3] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2)};
316  double DensEl = CUB(NSample)/(pVol()*NPairF);
317  double Dens = 13.3;
318  FILE *FNano = fopen("NanoPos.txt","w");
319  //IsoLevel *= NPairF;
320  for(int ff=NFile[0];ff<NFile[1];ff+=NPairF){
321  double *Plot = (double *)calloc(CUBE(NSample),sizeof(double));
322  double Min = 0.;
323  double Max = 0.;
324  VAR_TRIANGLE *Triang = NULL;
325  double Pos[3];
326  for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
327  Processing(f);
328  if(OpenRisk(cFile[f],BF_PART))return;
329  for(int b=0;b<pNBlock();b++){
330  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
331  if(pType(p) != 0)continue;
332  for(int d=0;d<3;d++){
333  Pos[d] = pPos(p,d) - (pNanoPos(0,d)-.5*pEdge(d));
334  Pos[d] -= floor(Pos[d]*pInvEdge(d))*pEdge(d);
335  }
336  int sx = (int)(Pos[0]*pInvEdge(0)*NSample);
337  int sy = (int)(Pos[1]*pInvEdge(1)*NSample);
338  int sz = (int)(Pos[2]*pInvEdge(2)*NSample);
339  int sTot = (sx*NSample+sy)*NSample+sz;
340  Plot[sTot] += DensEl;
341  if(Max < Plot[sTot]) Max = Plot[sTot];
342  if(Min > Plot[sTot]) Min = Plot[sTot];
343  }
344  }
345  }
346  Matrice Mask(3,3,3);
347  Mask.FillGaussian(.5,3.);
348  Mask.Print();
349  int NDim = 3;
350  int IfMinImConv = 1;
351  Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
352  Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
353  // ConvoluteMatrix(Plot,NSample,&Mask,3);
354  // ConvoluteMatrix(Plot,NSample,&Mask,3);
355  char FString[256];
356  sprintf(FString,"IsoSurf%05d.dat",ff);
357  FILE *F2Write = fopen(FString,"w");
358  fprintf(F2Write,"#l(%lf %lf %lf) v[%d] d[polygon]\n",pEdge(0),pEdge(1),pEdge(2),NSample);
359  HeaderNano(F2Write);
360  int NTri = 0;
361  for(int i=0;i<NIso;i++){
362  printf("Min %lf Max %lf IsoLevel %lf DensEl %lf\n",Min,Max,IsoLevel[i],DensEl);
363  Triang = MarchingCubes(Plot,NSample,IsoLevel[i],&NTri);
364  for(int t=0;t<NTri;t++){
365  for(int i=0;i<3;i++){
366  int l1 = t*3 + (i+1)%3;
367  int l2 = t*3 + (i+2)%3;
368  for(int d=0;d<3;d++) Pos[d] = Triang[t].p[i].x[d];
369  int sx = (int)(Pos[0]*pInvEdge(0)*NSample);
370  int sy = (int)(Pos[1]*pInvEdge(1)*NSample);
371  int sz = (int)(Pos[2]*pInvEdge(2)*NSample);
372  int sTot = (sx*NSample+sy)*NSample+sz;
373  fprintf(F2Write,"{t[%d %d %d] ",sTot,t,0);
374  fprintf(F2Write,"x(%lf %lf %lf) ",Pos[0],Pos[1],Pos[2]);
375  fprintf(F2Write,"v(%lf %lf %lf) ",Triang[t].n[i].x[0],Triang[t].n[i].x[1],Triang[t].n[i].x[2]);
376  fprintf(F2Write,"l[%d] l[%d]}\n",l1,l2);
377  }
378  }
379  }
380  fclose(F2Write);
381  free(Plot);
382  free(Triang);continue;
383  int NVertex = CUBE(2*NSample-1);
384  double *WeightL = (double *) calloc(NVertex,sizeof(double));
385  NormalWeight(Triang,WeightL,NSample,NTri);
386  double CmStalk[3] = {0.,0.,0.};//OldPos[0],OldPos[1],OldPos[2]};
387  double CountStalk = 0.;
388  for(int t=0;t<NTri;t++){
389  for(int v=0;v<3;v++){
390  int vRef = Triang[t].v[v];
391  for(int d=0;d<3;d++){
392  CmStalk[d] = Triang[t].p[v].x[d]*WeightL[vRef];
393  }
394  CountStalk += WeightL[vRef];
395  }
396  }
397  free(WeightL);
398  free(Triang);
399  if(CountStalk <= 0.) CountStalk = 1.;
400  for(int d=0;d<3;d++){
401  CmStalk[d] /= CountStalk;
402  }
403  pPos(CmStalk);
404  SetNNano(1);
405  for(int d=0;d<3;d++){
406  Nano->Pos[d] = CmStalk[d];
407  OldPos[d] = Nano->Pos[d];
408  }
409  SetNanoBkf(0);
410  Nano->Axis[0] = 0.;
411  Nano->Axis[1] = 0.;
412  Nano->Axis[2] = 1.;
413  Nano->Rad = .5;
414  Nano->Height = 4.;
415  Nano->Hamaker = 1.;
416  Nano->OffSet = 1.;
417  Nano->Shape = SHAPE_STALK;
418  for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
419  char String[120];
420  StringNano(String,0);
421  fprintf(FNano,"sed '/Rigid/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
422  fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
423  //command("chmod u+x %s\n","NanoPos.txt");
424  //SubNanoHeader(cFile[f]);
425  }
426  printf("\n");
427  }
428  fclose(FNano);
429 }
430 void ElPoly::IsoLine(int NSample,double *IsoLevel,int NIso,int How){
431  int NPairF = 1;//NFile[1]-NFile[0];
432  double OldPos[3] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2)};
433  double DensEl = SQR(NSample)/(pVol()*NPairF);
434  double Dens = 13.3;
435  //IsoLevel *= NPairF;
436  for(int ff=NFile[0];ff<NFile[1];ff+=NPairF){
437  double *Plot = (double *)calloc(SQR(NSample),sizeof(double));
438  double Min = 0.;
439  double Max = 0.;
440  double Pos[3];
441  for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
442  Processing(f);
443  if(OpenRisk(cFile[f],BF_PART))return;
444  for(int b=0;b<pNBlock();b++){
445  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
446  if(pType(p) != 0)continue;
447  for(int d=0;d<3;d++){
448  Pos[d] = pPos(p,d) - (pNanoPos(0,d)-.5*pEdge(d));
449  Pos[d] -= floor(Pos[d]*pInvEdge(d))*pEdge(d);
450  }
451  int sx = (int)(Pos[CLat1]*pInvEdge(CLat1)*NSample);
452  int sy = (int)(Pos[CLat2]*pInvEdge(CLat2)*NSample);
453  int sTot = sx*NSample+sy;
454  if(How==0)//particle file
455  Plot[sTot] += DensEl;
456  else
457  Plot[sTot] += pPos(p,CNorm);
458  if(Max < Plot[sTot]) Max = Plot[sTot];
459  if(Min > Plot[sTot]) Min = Plot[sTot];
460  }
461  }
462  printf("Min %lf Max %lf DensEl %lf\n",Min,Max,DensEl);
463  }
464  Matrice Mask(5,5);
465  Mask.FillGaussian(.5,3.);
466  Mask.Print();
467  int NDim = 2;
468  int IfMinImConv = 1;
469  Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
470  Mask.ConvoluteMatrix(Plot,NSample,NDim,IfMinImConv);
471  char FString[60];
472  sprintf(FString,"IsoSurf%05d.dat",ff);
473  FILE *F2Write = fopen(FString,"w");
474  fprintf(F2Write,"#l(%lf %lf %lf) v[%d] d[part]\n",pEdge(0),pEdge(1),pEdge(2),NSample);
475  HeaderNano(F2Write);
476  IsoLine(F2Write,Plot,NSample,IsoLevel,NIso);
477  free(Plot);
478  fclose(F2Write);
479  }
480 }
481 void ElPoly::IsoLine(FILE *F2Write,double *Plot,int NSample,double *IsoLevel,int NIso){
482  VAR_LINE *Triang = NULL;
483  int NTri = 0;
484  int p=0;
485  char FName[60];
486  double Pos[3];
487  for(int i=0;i<NIso;i++){
488  Triang = MarchingSquares(Plot,NSample,IsoLevel[i],&NTri);
489  ConnectLineChain(Triang,NSample,NTri);
490  sprintf(FName,"IsoLineChain%d.dat",i);
491  sprintf(cWhat2Draw,"part");
492  Write(FName);
493  // for(int c=0;c<pNChain();c++){
494  // sprintf(FString,"Line%05dChain%02d.dat",ff,c);
495  // FILE *FLine = fopen(FString,"w");
496  // while(
497  // for(int p=0;p<pNChain();p++){
498  // if
499  // fprintf(FString,"%lf %lf\n"
500  // fclose(FLine);
501  // }
502  if(1==1){//Isoline without ordering
503  for(int t=0;t<NTri;t++){
504  for(int v=0;v<2;v++){
505  for(int d=0;d<3;d++) Pos[d] = Triang[t].p[v].x[d];
506  int sTot = Triang[t].v[v];
507  fprintf(F2Write,"{t[%d %d %d] ",sTot,t,0);
508  fprintf(F2Write,"x(%lf %lf %lf) ",Pos[0],Pos[1],Pos[2]);
509  if(v==0)
510  fprintf(F2Write,"l[%d]}\n",p+1);
511  else
512  fprintf(F2Write,"\n");
513  p++;
514  }
515  }
516  }
517  if(1==0){//Isoline separated between inside and outside an elips
518  double Radx = 20.;
519  double Rady = 12.;
520  double Center[3] = {0.,.5*pEdge(CLat2),0.};
521  int NIt = 100;
522  double DistElips = SQR(Radx)+SQR(Rady);
523  for(int i=0;i<NIt;i++){
524  double Ang = i/(double)NIt*M_PI - M_PI*.5;
525  double x = Center[0] + Radx*cos(Ang);
526  double y = Center[1] + Rady*sin(Ang);
527  fprintf(F2Write,"{t[%d %d %d] ",p++,2,2);
528  fprintf(F2Write,"x(%lf %lf %lf) }\n",x,y,0.);
529  }
530  for(int t=0;t<NTri;t++){
531  for(int v=0;v<2;v++){
532  for(int d=0;d<3;d++) Pos[d] = Triang[t].p[v].x[d];
533  int sTot = Triang[t].v[v];
534  int type = 0;
535  if(fabs(Pos[0] - Center[0]) < Radx && fabs(Pos[1] - Center[1]) < Rady) type = 1;
536  fprintf(F2Write,"{t[%d %d %d] ",sTot,t,type);
537  fprintf(F2Write,"x(%lf %lf %lf) ",Pos[0],Pos[1],Pos[2]);
538  if(v==0)
539  fprintf(F2Write,"l[%d]}\n",p+1);
540  else
541  fprintf(F2Write,"\n");
542  p++;
543  }
544  }
545  }
546  }
547  if(1==1){//to write the continuum density values
548  for(int sx=0;sx<NSample;sx++){
549  for(int sy=0;sy<NSample;sy++){
550  double x = pEdge(CLat1)*sx/(double)NSample;
551  double y = pEdge(CLat2)*sy/(double)NSample;
552  int sTot = sx*NSample+sy;
553  fprintf(F2Write,"{t[%d %d %d] ",sTot,NTri+sTot/2,1);
554  fprintf(F2Write,"x(%lf %lf %lf) }\n",x,y,Plot[sx*NSample+sy]);
555  }
556  }
557  }
558  if(1==0){//Write the chains separately
559  for(int c=0;c<pNChain();c++){
560  sprintf(FName,"Chain%d.dat",c);
561  FILE *FChain = fopen(FName,"w");
562  for(int p=0;p<pNPart();p++){
563  if(Pm[p].CId != c) continue;
564  fprintf(FChain,"%lf %lf\n",Pm[p].Pos[0],Pm[p].Pos[1]);
565  }
566  }
567  }
568  free(Triang);
569 }
571  FILE *FNano = fopen("NanoPos.sh","w");
572  double OldPos[5] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2),Nano->Rad,Nano->Height};
573  FILE *StalkArea = fopen("PoreArea.dat","w");
574  fprintf(StalkArea,"#time rad asy\n");
575  for(int f=NFile[0];f<NFile[1];f++){
576  Processing(f);
577  if(OpenRisk(cFile[f],BF_PART)) return;
578  double Asy = PorePos();
579  Nano[0].Shape = ShapeId("pore");
580  Nano->Height = .2;
581  char String[120];
582  StringNano(String,0);
583  fprintf(StalkArea,"%lf %lf %lf\n",pTime(),Nano->Rad,Asy);
584  fprintf(FNano,"sed '/pore/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
585  fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
586  //command("chmod u+x %s\n","NanoPos.txt");
587  //SubNanoHeader(cFile[f]);
588  }
589  fclose(FNano);
590  fclose(StalkArea);
591  printf("\n");
592 }
597  FILE *FNano = fopen("NanoPos.sh","w");
598  double OldPos[5] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2),Nano->Rad,Nano->Height};
599  FILE *StalkArea = fopen("StalkArea.dat","w");
600  char FName[60];
601  for(int f=NFile[0];f<NFile[1];f++){
602  Processing(f);
603  if(OpenRisk(cFile[f],BF_NO)) return;
604  SetNNano(1);
605  if(StalkPos(OldPos)) continue;
606  char String[120];
607  StringNano(String,0);
608  fprintf(StalkArea,"%lf %lf\n",pTime(),Nano->Area);
609  fprintf(FNano,"sed '/Rigid/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
610  fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
611  sprintf(FName,"Centered%05d.dat",f);
612  //Write(FName);
613  //command("chmod u+x %s\n","NanoPos.txt");
614  //SubNanoHeader(cFile[f]);
615  }
616  fclose(FNano);
617  fclose(StalkArea);
618  printf("\n");
619 }
625  FILE *FNano = fopen("NanoPosA.sh","w");
626  FILE *StalkArea = fopen("StalkArea.dat","w");
627  FILE *AreaS = fopen("AreaS.dat","w");
628  double OldPos[5] = {pNanoPos(0,0),pNanoPos(0,1),pNanoPos(0,2),Nano->Rad,Nano->Height};
629  int NBin = 36;
630  double *Count = (double *)calloc(NBin*NBin,sizeof(double));
631  double RadTorus = 1.;//Nano->Rad;
632  //fprintf(AreaS,"#l(%lf %lf %lf) \n",2.*Nano->Height,2.*Nano->Height,10.);
633  fprintf(AreaS,"#l(%lf %lf %lf) \n",pEdge(0),pEdge(1),pEdge(2));
634  HeaderNano(AreaS);
635  char FName[60];
636  for(int f=NFile[0];f<NFile[1];f++){
637  Processing(f);
638  if(OpenRisk(cFile[f],BF_NO)) return;
639  SetNNano(1);
640  if(StalkPos(OldPos)) continue;
641  //Nano->Pos[CNorm] = pCm(CNorm);
642  //Nano->Pos[CNorm] -= floor(Nano->Pos[CNorm]*pInvEdge(CNorm))*pEdge(CNorm);
643  double Cm[3] = {0.,0.,0.};
644  double CountCm = 0;
645  int nPart = 0;
646  double Pos[3];
647  //counts the particles inside the torus
648  for(int p=0;p<pNPart();p++){
649  for(int d=0;d<3;d++){
650  Pos[d] = pPos(p,d) - Nano->Pos[d];
651  Pos[d] -= floor(pPos(p,d)*pInvEdge(d))*pEdge(d);
652  }
653  double Rad = hypot(Pos[CLat1],Pos[CLat2]);
654  if(Rad > Nano->Height) continue;
655  if(fabs(Pos[CNorm]) > RadTorus) continue;
656  // fprintf(AreaS,"{t[%d 0 %d] x(%lf %lf %lf)",nPart++,pType(p),pPos(p,0),pPos(p,1),pPos(p,2));
657  // if(Ln[p].NLink > 0) fprintf(AreaS,"l[%d]",p-Ln[p].Link[0]+nPart+1);
658  // fprintf(AreaS,"}\n");
659  if(pType(p) == 1) continue;
660  for(int d=0;d<3;d++){
661  Cm[d] += pPos(p,d);
662  }
663  CountCm += 1.;
664  int vx = (int)(Pos[CLat1]/(2.*Nano->Height)*NBin);
665  vx += NBin/2;
666  if(vx < 0 || vx >= NBin) continue;
667  int vy = (int)(Pos[CLat2]/(2.*Nano->Height)*NBin);
668  vy += NBin/2;
669  if(vy < 0 || vy >= NBin) continue;
670  Count[vx*NBin+vy] += 1.;
671  }
672  double Area = 0.;
673  double Norm = 0.;
674  for(int vx=0;vx<NBin;vx++){
675  for(int vy=0;vy<NBin;vy++){
676  if(Count[vx*NBin+vy] < 1.) continue;
677  // double x = vx*Nano->Height*2./(double)NBin;
678  // double y = vy*Nano->Height*2./(double)NBin;
679  // fprintf(AreaS,"{t[%d 0 2] x(%lf %lf %lf)}\n",nPart++,x+pNanoPos(0,0)-Nano->Height,y+pNanoPos(0,1)-Nano->Height,pNanoPos(0,2));
680  Area += 1.;
681  }
682  }
683  if(CountCm <= 0.){
684  printf("No particles in the torus %s\n",cFile[f]);
685  return;
686  }
687  Nano->Area = SQR(2.*Nano->Height)*Area/(double)(SQR(NBin));
688  for(int d=0;d<3;d++){
689  Cm[d] /= CountCm;
690  }
691  Cm[CNorm] = pCm(CNorm);
692  //fprintf(AreaS,"%lf %lf %lf\n",Cm[0]-Nano->Pos[0],Cm[1]-Nano->Pos[1],Cm[2]-Nano->Pos[2]);
693  for(int d=0;d<3;d++){
694  Nano->Pos[d] = Cm[d] - floor(Cm[d]*pInvEdge(d))*pEdge(d);
695  fprintf(StalkArea,"%lf %lf\n",pTime(),Nano->Area);
696  }
697  SetNanoBkf(0);
698  Nano->Height = Nano->Rad + sqrt(Nano->Area/DUE_PI);
699  char String[120];
700  StringNano(String,0);
701  fprintf(StalkArea,"%lf %lf\n",pTime(),Nano->Area);
702  fprintf(FNano,"sed '/Rigid/c\\%s' %s > Ciccia.dat\n",String,cFile[f]);
703  fprintf(FNano,"mv Ciccia.dat %s\n",cFile[f]);
704  sprintf(FName,"Centered%05d.dat",f);
705  //Write(FName);
706  //HeaderNano(AreaS);
707  }
708  fclose(AreaS);
709  fclose(StalkArea);
710  fclose(FNano);
711  free(Count);
712  printf("\n");
713 }
715  int NPairF = 10;
716  double InvAv = 1./(double)NPairF;
717  double Norm = 1./(double)(NFile[1]-NFile[0]);
718  char FName[256];
719  // for(int ff=NFile[0];ff<NFile[1];ff+=NPairF){
720  PART *Pn = (PART *)calloc(pNPart(),sizeof(PART));
721  //for(int f=ff;f<ff+NPairF&&f<NFile[1];f++){
722  for(int f=NFile[0];f<NFile[1];f++){
723  Processing(f);
724  if(OpenRisk(cFile[f],BF_NO))return;
725  for(int p=0;p<pNPart();p++){
726  for(int d=0;d<3;d++){
727  Pn[p].Pos[d] += pPos(p,d);
728  }
729  }
730  }
731  for(int p=0;p<pNPart();p++){
732  double Pos[3] = {Pn[p].Pos[0]*Norm,Pn[p].Pos[1]*Norm,Pn[p].Pos[2]*Norm};
733  SetPos(p,Pos);
734  }
735  printf("\n");
736  free(Pn);
737  sprintf(FName,"Av%s",cFile[NFile[0]]);
738  Write(FName);
739 }
740 void ElPoly::Sample(int NSample){
741  MOMENTI m1 = SampleSurfaceMem(NSample);
742  FILE *FWrite = fopen("SolSim.dat","w");
743  double InvNSample = 1./(double)NSample;
744  fprintf(FWrite,"# l(%.2f %.2f %.2f) d[part]\n",pEdge(0),pEdge(1),pEdge(2));
745  for(int sx=0;sx<NSample;sx++){
746  int sy = 0;
747  PlotMem[sx*NSample + sy] = m1.Uno;
748  PlotMem[sy*NSample + sx] = m1.Uno;
749  sy = NSample - 1;
750  PlotMem[sx*NSample + sy] = m1.Uno;
751  PlotMem[sy*NSample + sx] = m1.Uno;
752  }
753  printf("sample average: %lf\n",m1.Uno);
754  double ConvFact = 1.45/m1.Uno;
755  for(int sx=0;sx<NSample;sx++){
756  double x = sx*InvNSample*pEdge(CLat1);
757  int sx1 = sx + 1 == NSample ? 0 : sx + 1;
758  for(int sy=0;sy<NSample;sy++){
759  double y = sy*InvNSample*pEdge(CLat2);
760  int sy1 = sy + 1 == NSample ? 0 : sy + 1;
761  int p = sx*NSample+sy;
762  int l1 = sx1*NSample+sy + SQR(NSample);
763  int l2 = sx*NSample+sy1 + SQR(NSample);
764  //fprintf(FWrite,"{t[%d 0 2] x(%lf %lf %lf) l[%d] l[%d]}\n",p,x,y,PlotMem[p],sx1*NSample+sy,sx*NSample+sy1);
765  fprintf(FWrite,"{t[%d 0 3] x(%lf %lf %lf) l[%d] l[%d]}\n",p,x,y,(PlotMem[p]-m1.Uno)*ConvFact,l1,l2);
766  //fprintf(FWrite,"{t[%d 0 2] x(%lf %lf %lf)}\n",p,x,y,PlotMem[p]-m1.Uno);
767  }
768  }
769 }
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
double pTime()
Total time.
Definition: VarData.h:908
double Ypsilon
External parameter to calculate the contact angle.
Definition: Matematica.h:330
double x[3]
Cartesian coordinates.
Definition: VarData.h:486
void IsoLine(int NSample, double *IsoValue, int NIso, int How)
Calculate the discrete density of the system and the correspondent isolines.
Definition: ElPolyRepr.cpp:430
bool OpenRisk(char *InFile, int BF)
Opens a file without reallocationg.
Definition: VarData.cpp:126
double pVel(int p, int d)
Return the velocity.
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
void ShapeId(int iShape, char *Shape)
Identifier of the shape.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
Define a triangle.
Definition: VarData.h:509
double Hamaker
Strength of the interaction.
Definition: VarData.h:447
int pChain(int p)
Return the chain.
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file".
int Zeri(double a, double b, double *Radici, int NRadici)
Find the.
void StringNano(char *NString, int n)
String for the rigid inclusion in the header file.
char cWhat2Draw[STRSIZE]
What to draw.
Definition: VarData.h:1042
void SetPos(int p, double *Pos)
Set the particle position.
double Rad
Size.
Definition: VarData.h:445
void SetNanoBkf(int n)
Set the back folded array for the nano n.
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
double Height
Height of the cylinder.
Definition: VarData.h:449
Define a triangle.
Definition: VarData.h:498
MOMENTI SampleSurfaceMem(int NSample)
Allocate and fill PlotMem with the particle average position.
Definition: VarDataEl.cpp:49
int pType(int p)
Return the type.
double pCm(int d)
Center of mass of the system.
Definition: VarData.h:924
int CLat2
lateral coordinate
Definition: VarData.h:1078
int v[2]
Reference for the vectors.
Definition: VarData.h:515
VAR_TRIANGLE * MarchingCubes(double *Plot, int NSample, double IsoLevel, int *NTri)
Defines the triangles close to the IsoLevel of the 3d density Plot.
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
void ConnectLineChain(VAR_LINE *Triang, int NGrid, int NTri)
Connect the lines in a chain.
double * PlotMem
Particle position/density on the square lattice.
Definition: VarData.h:1056
XYZ p[3]
The three vertices.
Definition: VarData.h:500
void Sample(int NSample)
Sample the space in NSample lattice points.
Definition: ElPolyRepr.cpp:740
void ConvoluteMatrix(double *Plot, int NGrid, int NDim, int IfMinImConv)
Convolute with a matrix.
VAR_LINE * MarchingSquares(double *Plot, int NSample, double IsoLevel, int *NTri)
Defines the triangles close to the IsoLevel of the 3d density Plot.
double Area
Area of a pore or a stalk.
Definition: VarData.h:471
double pVol()
xyzr edges of the simulation box
Definition: VarData.h:922
FUNC Func
Pointer to a function.
Definition: Matematica.h:111
int v[3]
Reference for the vectors.
Definition: VarData.h:504
double Pos[3]
Position.
Definition: VarData.h:427
int Shape
0 none, 1 spherical, 2 cylindrical 3 wall
Definition: VarData.h:473
int From3To2d(int Coord, double Param)
Project the velocities of a 3d system on a 2d system wrt a coordinate.
Definition: ElPolyRepr.cpp:188
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
double PorePos()
Find the position of the pore.
Moments of a distribution.
void AvSnap()
Average the postion of the lipids over many snapshots.
Definition: ElPolyRepr.cpp:714
void FetchPore()
Write the position of the pore for every snapshot.
Definition: ElPolyRepr.cpp:570
int pNBlock()
Number of blocks.
int StalkPos(double *OldPos)
Find the position of the stalk.
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
void FetchStalk()
Write the position of the stalk for every snapshot.
Definition: ElPolyRepr.cpp:596
double ContactAngle(double x)
Definition of the contact angle.
int HeaderNano(FILE *FileToWrite)
Write the nano section of the header to the file.
bool ShiftSys(int How)
Shift the system accordin to the SHIFT_ definitions.
double PreFact
External parameter in the definition of the contact angle.
Definition: Matematica.h:332
int RadialShell(int NBin)
Outer shell of a projected (rad,normal) graph.
Definition: ElPolyRepr.cpp:239
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
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 OffSet
Reference potential.
Definition: VarData.h:457
void StalkArea()
Area of hydrophobic in the torus.
Definition: ElPolyRepr.cpp:624
double Uno
First moment.
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
double NormalWeight(VAR_TRIANGLE *Triang, double *Weight, int NGrid, int NTri)
Weight of the neighblorung normal on a vertex.
int EndIdx
End particle position.
Definition: VarData.h:267
double pPos(int p, int d)
Return back folded position.
int pNChain()
Number of chain.
int From3To1d(int Coord)
Project a 3d system on one coordinate.
Definition: ElPolyRepr.cpp:235
int NFile[2]
First and last file of the list.
Definition: ElPoly.h:418
void IsoSurf(int NSample, double *IsoValue, int NIso)
Density plot all over different snapshots and calculation of the isolevel surface.
Definition: ElPolyRepr.cpp:313
int CLat1
lateral coordinate
Definition: VarData.h:1076
Information of every particle.
Definition: VarData.h:214
int Surface(int NBin, int Coord)
Area of a surfuce projected on the coordinate Coord surf.
Definition: ElPolyRepr.cpp:111
double Axis[3]
Rotation axis.
Definition: VarData.h:435
int ProjectionF(int NBin, int Coord)
Projection of the system against one direction pro.
Definition: ElPolyRepr.cpp:2
int SetNNano(int Val)
Set NNano.
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
void Print()
Print the entries.
void Processing(int f)
Information on the current file elaborated.
Definition: ElPolyEl.cpp:369
int From2To1d(int Coord)
Projet a 2d system on one of the two coordinates.
Definition: ElPolyRepr.cpp:204
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.