Allink  v0.1
ElPolyProfile.cpp
1 #include "ElPoly.h"
2 
3 int ElPoly::Temperature(int NBin,int Coord){
4  double *dDensity;
5  dDensity = (double *)malloc(NBin*sizeof(double));
6  for(int v=0;v<NBin;v++){
7  dDensity[v] = 0.;
8  }
9  FILE *FileToWrite = fopen("Temperature.dat","w");
10  for(int f=NFile[0];f<NFile[1];f++){
11  Processing(f);
12  if(OpenRisk(cFile[f],BF_NANO)==0)return 0;
13  //TemperatureProfile(NBin,Coord);
14  for(int v=0;v<NBin;v++){
15  //dDensity[v] += dPoint[v];
16  }
17  }
18  printf("\n");
19  for(int v=0;v<NBin;v++){
20  fprintf(FileToWrite,"%lf \n",dDensity[v]/(double)(NFile[1] - NFile[0]));
21  }
22  //FileToWrite = fopen(OutFile,"w");
23  return 0;
24 }
25 int ElPoly::NanoParticle(int NBin){
26  //SetEdge(.5*sqrt(SQR(pEdge(CLat1))+SQR(pEdge(CLat2))),3);
27  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
28  int NType = 2;
29  double *dDensity = (double *)calloc(NType*NBin,sizeof(double));
30  double *dRad = (double *)calloc(NType*NBin,sizeof(double));
31  double *ThickProf = (double *)calloc(NType*NBin,sizeof(*ThickProf));
32  double *CountCh = (double *)calloc(NType*NBin,sizeof(double));
33  double *CountPart = (double *)calloc(NBin,sizeof(double));
34  double *AngleProf = (double *)calloc(2*NBin,sizeof(double));
35  double *VolEl = (double *)calloc(NBin,sizeof(double));
36  double *VelDistr = (double *)calloc(NBin*2,sizeof(double));
37  double *DisplRad = (double *)calloc(NBin*2,sizeof(double));
38  double *DisplTheta = (double *)calloc(NBin*2,sizeof(double));
39  double *E2E = (double *)calloc(NBin*2,sizeof(double));
40  double *HFluct = (double *)calloc(NBin*2,sizeof(double));
41  CHAIN *Ch2 = (CHAIN *)calloc(pNChain(),sizeof(CHAIN));
42  Vettore Directive(3);
43  Vettore Axis(3);
44  Vettore ZedDir(3);
45  for(int d=0;d<3;d++) ZedDir.Set(0.,d);
46  ZedDir.Set(1.,CNorm);
47  VolumeCircSlab(VolEl,NBin);
48  for(int f=NFile[0];f<NFile[1];f++){
49  Processing(f);
50  if(OpenRisk(cFile[f],BF_CHAIN)) return 0;
51  //ChainDiffusion(DIFF_LATERAL,ChDiff);
52  //Sum radial profile------------------------------------
53  for(int b=0;b<pNBlock();b++){
54  int IfAdded = 0;
55  if(!strncmp(Block[b].Name,"PEP",3) ) continue;
56  for(int p=Block[b].InitIdx,link=0;p<Block[b].EndIdx;p++){
57  if(pPos(p,CNorm) > pNanoPos(0,CNorm) + Nano->Rad*.5) continue;
58  if(pPos(p,CNorm) < pNanoPos(0,CNorm) - Nano->Rad*.5) continue;
59  double RadDist = sqrt( SQR(pPos(p,CLat1)-pNanoPos(0,CLat1)) + SQR(pPos(p,CLat2)-pNanoPos(0,CLat2)) );
60  double Zed = pPos(p,CNorm);
61  int vr = (int)(RadDist/pEdge(3)*NBin);
62  if(vr < 0 || vr >= NBin) continue;
63  int t = 0;
64  if(CHAIN_IF_TYPE(Ch[pType(p)].Type,CHAIN_ADDED) ) t = 1;
65  dRad[vr*NType+t] += 1.;
66  CountPart[vr] += 1.;
67  }
68  }
69  // DensityProfile(3,NBin,NType,dDensity);
70  // for(int v=0;v<NBin;v++){
71  // for(int t=0;t<NType;t++){
72  // dRad[v*NType+t] += dDensity[v*NType+t];
73  // dDensity[v*NType+t] = 0.;
74  // }
75  // }
76  //Sum thickness profile--------------------------------
77  for(int b=0,cOff=0;b<pNBlock();cOff+=pNChain(b++)){
78  if(!strncmp(Block[b].Name,"PEP",3)){
79  continue;
80  }
81  for(int c=cOff;c<cOff+pNChain(b);c++){
82  double RelPos[3];
83  for(int d=0;d<3;d++){
84  RelPos[d] = Ch[c].Pos[d] - pNanoPos(0,d);
85  RelPos[d] -= floor(RelPos[d]*pInvEdge(d) + .5)*pEdge(d);
86  Directive.Set(RelPos[d],d);
87  Axis.Set(Ch[c].Dir[d],d);
88  }
89  double RadDist = sqrt( SQR(RelPos[CLat1]) + SQR(RelPos[CLat2]) );
90  double Zed = (Ch[c].Pos[CNorm] - pCm(CNorm));
91  int vr = (int)(RadDist*pInvEdge(3)*NBin);
92  if(vr < 0 || vr >= NBin) continue;
93  int t = 0;
94  if(VAR_IF_TYPE(Ch[c].Type,CHAIN_DOWN)) t = 1;
95  CountCh[(vr)*2+t] += 1;
96  ThickProf[(vr)*2+t] += Zed;
97  HFluct[vr*2 ] += Zed;
98  HFluct[vr*2+1] += SQR(Zed);
99  double Angle = ZedDir.Angle(&Axis,&Directive);
100  Angle = Angle < .5*M_PI ? Angle : Angle - .5*M_PI;
101  if(!isnan(Angle)){
102  AngleProf[vr*2 ] += Angle;
103  AngleProf[vr*2+1] += SQR(Angle);
104  }
105  double VelRad = sqrt( SQR(Ch[c].Vel[CLat1]) + SQR(Ch[c].Vel[CLat2]) );
106  double VelTheta = fabs(atan2(Ch[c].Vel[CLat2],Ch[c].Vel[CLat1]));
107  double DRad = sqrt(SQR(Ch[c].Pos[CLat1]-Ch2[c].Pos[CLat1]) + SQR(Ch[c].Pos[CLat2]-Ch2[c].Pos[CLat2]));
108  double DTheta = atan2(Ch[c].Pos[CLat1]-Ch2[c].Pos[CLat1],Ch[c].Pos[CLat2]-Ch2[c].Pos[CLat2]);
109  VelDistr[vr*2 ] += VelRad;
110  VelDistr[vr*2+1] += VelTheta < 100. ? VelTheta : 0.;
111  DisplRad[vr*2 ] += DRad;
112  DisplRad[vr*2+1] += SQR(DRad);
113  DisplTheta[vr*2 ] += DTheta;
114  DisplTheta[vr*2+1] += SQR(DTheta);
115  double DistE2E = 0.;
116  int pE1 = c*pNPCh();
117  int pE2 = (c+1)*pNPCh()-1;
118  for(int d=0;d<3;d++){
119  Ch2[c].Pos[d] = Ch[c].Pos[d];
120  DistE2E += SQR(pPos(pE1,d)-pPos(pE2,d));
121  }
122  E2E[vr*2 ] += sqrt(DistE2E);
123  E2E[vr*2+1] += DistE2E;
124  }
125  }
126  }
127 #ifdef OMPI_MPI_H
128  MPI_Allreduce(MPI_IN_PLACE,dDensity,NType*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
129  MPI_Allreduce(MPI_IN_PLACE,ThickProf,NType*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
130  MPI_Allreduce(MPI_IN_PLACE,CountCh,NType*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
131  MPI_Allreduce(MPI_IN_PLACE,CountPart,NType*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
132  MPI_Allreduce(MPI_IN_PLACE,VelDistr,2*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
133  MPI_Allreduce(MPI_IN_PLACE,DisplRad,2*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
134  MPI_Allreduce(MPI_IN_PLACE,DisplTheta,2*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
135  MPI_Allreduce(MPI_IN_PLACE,E2E,2*NBin, MPI_DOUBLE, MPI_SUM, Proc->CommGrid);
136  int Rank=0;
137  MPI_Comm_rank(Proc->CommGrid, &Rank);
138  if(Rank==0){
139 #endif
140  printf("\n");
141  //Normalize-----------------------------------
142  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
143  double VolumeZed = pEdge(0)*pEdge(1)*pEdge(2)/(double)SQR(NBin);
144  double VolumeRad = pEdge(3)*pEdge(3)*Nano->Rad/(double)SQR(NBin);
145  for(int v=0;v<NBin;v++){
146  double r = v*pEdge(3)/(double)NBin;
147  double NormType = 0.;
148  double NormVol = (VolumeRad*VolEl[v])/InvNFile;
149  NormVol = NormVol > 0. ? 1./NormVol : 1.;
150  for(int t=0;t<NType;t++){
151  dRad[v*NType+t] *= NormVol;
152  double Norm = CountCh[v*NType+t] > 0. ? 1./CountCh[v*NType+t] : 1.;
153  ThickProf[v*NType+t] *= Norm;
154  NormType += CountCh[v*NType+t];
155  }
156  NormType = NormType > 0. ? 1./NormType : 1.;
157  if(r < Nano->Rad) NormType = 0.;
158  if(r < Nano->Rad) NormVol = 0.;
159  VelDistr[v*2 ] *= NormType;
160  VelDistr[v*2+1] *= NormType;
161  DisplRad[v*2 ] *= NormType;
162  DisplRad[v*2+1] = sqrt(DisplRad[v*2+1]*NormType - SQR(DisplRad[v*2 ]));
163  DisplTheta[v*2 ] *= NormType;
164  DisplTheta[v*2+1] = sqrt(DisplTheta[v*2+1]*NormType - SQR(DisplTheta[v*2 ]));
165  HFluct[v*2 ] *= NormType;
166  HFluct[v*2+1] = sqrt(HFluct[v*2+1]*NormType - SQR(HFluct[v*2 ]));
167  AngleProf[v*2 ] *= NormType;
168  AngleProf[v*2+1] = sqrt(AngleProf[v*2+1]*NormType - SQR(AngleProf[v*2 ]));
169  E2E[v*2 ] *= NormType;
170  E2E[v*2+1] = sqrt( (E2E[v*2+1]*NormType - SQR(E2E[v*2 ])) );
171  //E2E[v*2+1] *= NormVol;
172  }
173  //--------------radial density-------------------------------
174  char FileName[256];
175  sprintf(FileName,"NanoDensR%.1fS%.1fH%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
176  FILE *FileToWrite = fopen(FileName,"w");
177  char cSystem[STRSIZE];
178  SysDef(cSystem);
179  fprintf(FileToWrite,"# %s",cSystem);
180  fprintf(FileToWrite,"#r DensPhob DensPhil\n");
181  for(int v=0;v<NBin;v++){
182  double r = v*pEdge(3)/(double)NBin;
183  fprintf(FileToWrite,"%lf %lf %lf\n",r,dRad[v*NType],dRad[v*NType+1]);
184  }
185  fclose(FileToWrite);
186  // //---------------thickness profile-----------------------------
187  // sprintf(FileName,"NanoThickR%.1fS%.1fH%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
188  // FileToWrite = fopen(FileName,"w");
189  // fprintf(FileToWrite,"#r ThickUp ThickDown\n");
190  // for(int vr=0;vr<NBin;vr++){
191  // double r = vr*pEdge(3)/(double)NBin;
192  // if(r < Nano->Rad) continue;
193  // if(ThickProf[vr*NType+0]*ThickProf[vr*NType+1] > 0.)
194  // fprintf(FileToWrite,"%lf %lf %lf\n",r,ThickProf[vr*NType+0],ThickProf[vr*NType+1]);
195  // }
196  // fclose(FileToWrite);
197  //----------------velocity profile-----------------------------
198  sprintf(FileName,"NanoRadVelR%.1fS%.1fH%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
199  FileToWrite = fopen(FileName,"w");
200  fprintf(FileToWrite,"#r VelRad VelTheta\n");
201  for(int v=0;v<NBin;v++){
202  double r = v*pEdge(3)/(double)NBin;
203  fprintf(FileToWrite,"%lf %lf %lf\n",r,VelDistr[v*2+0],VelDistr[v*2+1]);
204  }
205  fclose(FileToWrite);
206  //----------------z fluctuation profile-----------------------------
207  sprintf(FileName,"NanoZedFluctR%.1fS%.1fH%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
208  FileToWrite = fopen(FileName,"w");
209  fprintf(FileToWrite,"#r <z> <(z - <z>)^2> \n");
210  for(int v=0;v<NBin;v++){
211  double r = v*pEdge(3)/(double)NBin;
212  fprintf(FileToWrite,"%lf %.5g %.5g\n",r,HFluct[v*2+0],HFluct[v*2+1]);
213  }
214  fclose(FileToWrite);
215  //----------------angle profile-----------------------------
216  sprintf(FileName,"NanoAngleR%.1fS%.1fH%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
217  FileToWrite = fopen(FileName,"w");
218  fprintf(FileToWrite,"#r <a> <(a - <a>)^2> \n");
219  for(int v=0;v<NBin;v++){
220  double r = v*pEdge(3)/(double)NBin;
221  fprintf(FileToWrite,"%lf %.5g %.5g\n",r,AngleProf[v*2+0],AngleProf[v*2+1]);
222  }
223  fclose(FileToWrite);
224  //----------------end to end-----------------------------
225  sprintf(FileName,"NanoE2ER%.1fS%.1fH%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
226  FileToWrite = fopen(FileName,"w");
227  fprintf(FileToWrite,"#r <e2e> <(e2e - <e2e>)^2> \n");
228  for(int v=0;v<NBin;v++){
229  double r = v*pEdge(3)/(double)NBin;
230  fprintf(FileToWrite,"%lf %.5g %.5g\n",r,E2E[v*2+0],E2E[v*2+1]);
231  }
232  fclose(FileToWrite);
233  //----------------radial displacement-----------------------------
234  sprintf(FileName,"NanoRadDisplR%.1fS%.1fH%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
235  FileToWrite = fopen(FileName,"w");
236  fprintf(FileToWrite,"#r JumpRad JumpTheta \n");
237  for(int v=0;v<NBin;v++){
238  double r = v*pEdge(3)/(double)NBin;
239  fprintf(FileToWrite,"%lf %lf %lf %lf %lf\n",r,DisplRad[v*2 ],DisplRad[v*2+1],DisplTheta[v*2 ],DisplTheta[v*2+1]);
240  }
241  fclose(FileToWrite);
242 #ifdef OMPI_MPI_H
243  }
244 #endif
245  //freeing-------------------------------------------------
246  free(CountCh);
247  free(CountPart);
248  free(DisplRad);
249  free(DisplTheta);
250  free(dDensity);
251  free(dRad);
252  free(ThickProf);
253  free(VelDistr);
254  free(AngleProf);
255  return 0;
256 }
257 int ElPoly::WormF(int Partition,int NBin){
258  double *dDensity = (double *)calloc(2*NBin,sizeof(double));
259  double Norma;
260  double Border[2] = {-10.,10.};
261  for(int f=NFile[0];f<NFile[1];f++){
262  Processing(f);
263  if(OpenRisk(cFile[f],BF_PART))return 0;
264  Worm(Partition,NBin,Border,dDensity);
265  }
266  printf("\n");
267  FILE *FileToWrite = fopen("WormDensity.dat","w");
268  for(int v=0;v<NBin;v++){
269  fprintf(FileToWrite,"%lf %lf %lf\n",(double)v/(double)NBin*(Border[1]-Border[0])+Border[0],dDensity[v*2]/(NFile[1]-NFile[0]),dDensity[v*2+1]/(NFile[1]-NFile[0]));
270  }
271  fclose(FileToWrite);
272  free(dDensity);
273  return 1;
274 }
275 //The skin is not closing, better use the marching cubes
276 void ElPoly::StalkF(int NSample){
277  double Threshold = 0.;
278  int NLevel = 4;
279  double **Plot = (double **) calloc(NLevel,sizeof(double));
280  for(int l=0;l<NLevel;l++){
281  Plot[l] = (double *)calloc(NSample*NSample,sizeof(double));
282  }
283  double **PlotA = (double **) calloc(NLevel,sizeof(double));
284  for(int l=0;l<NLevel;l++){
285  PlotA[l] = (double *)calloc(NSample*NSample,sizeof(double));
286  }
287  for(int f=NFile[0];f<NFile[1];f++){
288  Processing(f);
289  if(OpenRisk(cFile[f],BF_PART))return ;
290  Stalk(NSample,NLevel,Plot,Threshold);
291  for(int s=0;s<NSample;s++)
292  for(int ss=0;ss<NSample;ss++)
293  for(int l=0;l<NLevel;l++)
294  PlotA[l][s*NSample+ss] += Plot[l][s*NSample+ss];
295  }
296  printf("\n");
297  Vettore v00(3);
298  Vettore v01(3);
299  Vettore v11(3);
300  Vettore v10(3);
301  Vettore vN(3);
302  char FileName[256];
303  sprintf(FileName,"Stalk.xvl");
304  FILE *FileToWrite = fopen(FileName,"w");
305  fprintf(FileToWrite,"#l(%lf %lf %lf) v[%d] d[%s]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NSample,ChooseDraw(EL_SKIN));
306  double FactC1 = pEdge(CLat1)/(double)NSample;
307  double FactC2 = pEdge(CLat2)/(double)NSample;
308  double FactCN = 1./(double)(NFile[1]-NFile[0]);
309  for(int l=0,p=0,c=0;l<NLevel;l++){
310  for(int ss=0;ss<NSample-1;ss++){
311  for(int s=0;s<NSample-1;s++){
312  v00.x[CLat1]=(s+.5)*FactC1;
313  v00.x[CLat2]=(ss+.5)*FactC2;
314  v00.x[CNorm]=PlotA[l][s*NSample+ss]*FactCN;
315  v01.x[CLat1]=(s+1.5)*FactC1;
316  v01.x[CLat2]=(ss+.5)*FactC2;
317  v01.x[CNorm]=PlotA[l][(s+1)*NSample+ss]*FactCN;
318  v11.x[CLat1]=(s+1.5)*FactC1;
319  v11.x[CLat2]=(ss+1.5)*FactC2;
320  v11.x[CNorm]=PlotA[l][(s+1)*NSample+ss+1]*FactCN;
321  v10.x[CLat1]=(s+.5)*FactC1;
322  v10.x[CLat2]=(ss+1.5)*FactC2;
323  v10.x[CNorm]=PlotA[l][s*NSample+ss+1]*FactCN;
324  //------------defines-the-squares---------------------
325  if(PlotA[l][(s)*NSample+ss] <= Threshold){
326  if(l==0)
327  v00.x[CNorm]=PlotA[3][(s)*NSample+ss+1]*FactCN;
328  if(l==1)
329  v00.x[CNorm]=PlotA[2][(s)*NSample+ss+1]*FactCN;
330  if(l==2)
331  v00.x[CNorm]=PlotA[1][(s)*NSample+ss+1]*FactCN;
332  if(l==3)
333  v00.x[CNorm]=PlotA[0][(s)*NSample+ss+1]*FactCN;
334  if(v00.x[CNorm] <= 0.) continue;
335  }
336  if(PlotA[l][(s+1)*NSample+ss] <= Threshold){
337  if(l==0)
338  v01.x[CNorm]=PlotA[3][(s)*NSample+ss]*FactCN;
339  if(l==1)
340  v01.x[CNorm]=PlotA[2][(s)*NSample+ss]*FactCN;
341  if(l==2)
342  v01.x[CNorm]=PlotA[1][(s)*NSample+ss]*FactCN;
343  if(l==3)
344  v01.x[CNorm]=PlotA[0][(s)*NSample+ss]*FactCN;
345  if(v01.x[CNorm] <= 0.) continue;
346  }
347  if(PlotA[l][s*NSample+ss+1] <= Threshold){
348  if(l==0)
349  v10.x[CNorm]=PlotA[3][(s)*NSample+ss]*FactCN;
350  if(l==1)
351  v10.x[CNorm]=PlotA[2][(s)*NSample+ss]*FactCN;
352  if(l==2)
353  v10.x[CNorm]=PlotA[1][(s)*NSample+ss]*FactCN;
354  if(l==3)
355  v10.x[CNorm]=PlotA[0][(s)*NSample+ss]*FactCN;
356  if(v10.x[CNorm] <= 0.) continue;
357  }
358  if(PlotA[l][(s+1)*NSample+(ss+1)] <= Threshold){
359  if(l==0)
360  v11.x[CNorm]=PlotA[3][(s)*NSample+ss+1]*FactCN;
361  if(l==1)
362  v11.x[CNorm]=PlotA[2][(s)*NSample+ss+1]*FactCN;
363  if(l==2)
364  v11.x[CNorm]=PlotA[1][(s)*NSample+ss+1]*FactCN;
365  if(l==3)
366  v11.x[CNorm]=PlotA[0][(s)*NSample+ss+1]*FactCN;
367  if(v11.x[CNorm] <= 0.) continue;
368  }
369  fprintf(FileToWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf) l[%d] l[%d] l[%d] }\n",p+0,c,l,v00.x[CLat1],v00.x[CLat2],v00.x[CNorm],(double)l*.25,.5,.8,p+1,p+2,p+3);
370  fprintf(FileToWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf) l[%d] l[%d] l[%d] }\n",p+1,c,l,v01.x[CLat1],v01.x[CLat2],v01.x[CNorm],(double)l*.25,.5,.8,p+2,p+3,p+0);
371  fprintf(FileToWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf) l[%d] l[%d] l[%d] }\n",p+2,c,l,v11.x[CLat1],v11.x[CLat2],v11.x[CNorm],(double)l*.25,.5,.8,p+3,p+0,p+1);
372  fprintf(FileToWrite,"{t[%d %d %d] x(%lf %lf %lf) v(%lf %lf %lf) l[%d] l[%d] l[%d] }\n",p+3,c,l,v10.x[CLat1],v10.x[CLat2],v10.x[CNorm],(double)l*.25,.5,.8,p+0,p+1,p+2);
373  p += 4;
374  c++;
375  }
376  }
377  }
378  for(int l=0;l<NLevel;l++){
379  free(Plot[l]);
380  free(PlotA[l]);
381  }
382  free(Plot);
383  free(PlotA);
384 }
388 void ElPoly::StalkLineProfF(int NBin){
389 #ifdef USE_FFTW
390  double *Line = (double *)calloc(NBin,sizeof(double));
391  double *Spe = (double *)calloc(NBin,sizeof(double));
392  char FName[120];
393  int NHalf = (int)(NBin/2);
394  // Excluding the two ends that are at most di
395  double InvNBin = 1./(double)(NBin-2);
396  double InvNFile = 1./(double)(NFile[1]-NFile[0]);
397  fftw_complex *out = (fftw_complex *)fftw_malloc(NBin*sizeof(fftw_complex));
398  fftw_complex *in = (fftw_complex *)fftw_malloc(NBin*sizeof(fftw_complex));
399  fftw_plan plan = fftw_plan_dft_1d(NBin-2,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
400  fftw_plan plan2 = fftw_plan_dft_1d(NBin-2,out,in,FFTW_BACKWARD,FFTW_ESTIMATE);
401  for(int f=NFile[0];f<NFile[1];f++){
402  Processing(f);
403  if(Open(cFile[f],BF_CHAIN)) return;
404  StalkLineProf(Line,NBin);
405  for(int vx=0;vx<NBin-2;vx++) in[vx][0] = Line[vx+1];
406  fftw_execute(plan);
407  for(int vx=0;vx<NBin-2;vx++){
408  //int vvx = vx < NHalf ? vx + NHalf : vx - NHalf;
409  Spe[vx] += (SQR(out[vx][0]) + SQR(out[vx][1]))*SQR(InvNBin);
410  }
411  //-------------------------write-temp-files
412  if(1==1){
413  sprintf(FName,"StalkLine%05d.dat",f);
414  FILE *FWrite = fopen(FName,"w");
415  for(int vx=0;vx<NBin;vx++){
416  fprintf(FWrite,"%lf %lf \n",vx*pEdge(CLat2)/(double)NBin,Line[vx]);
417  }
418  fclose(FWrite);
419  }
420  }
421  printf("\n");
422  double dx = pEdge(CLat1)*InvNBin;
423  FILE *FWrite = fopen("StalkLineSpectrum.dat","w");
424  for(int vx=0;vx<NBin/2-2;vx++){
425  double qx = 2.*M_PI*vx*pInvEdge(CLat1);
426  double qq = SQR(qx/dx - .5/dx);
427  //fprintf(FWrite,"%lf %lf \n",qx,Spe[vx]*InvNFile*pEdge(CLat2)*pEdge(CLat2));
428  fprintf(FWrite,"%lf %lf \n",qx,Spe[vx]*InvNFile);
429  }
430  fclose(FWrite);
431  free(Line);
432  free(Spe);
433  fftw_free(out);
434  fftw_free(in);
435 #else
436  printf("fftw not implemented\n");
437 #endif // USE_FFTW
438 }
440  char FName[60];
441  for(int f=NFile[0];f<NFile[1];f++){
442  Processing(f);
443  if(Open(cFile[f],BF_NO))return;
444  Nano[0].Pos[0] = 12.0;
445  Nano[0].Pos[1] = 12.0;
446  Nano[0].Pos[2] = 9.0;
447  sprintf(FName,"new%07d00.dat",f*2);
448  Write(cFile[f]);
449  }
450  printf("\n");
451 }
452 void ElPoly::RadNormPos(int NBin,int NGrid){
453  //SigErr((NFile[1]-NFile[0]) != SQR(NGrid),"RadNormPos: the number of files is not compatible with the parameter span");
454  char FName[120];
455  int NWeight = 5;
456  double *RadProf = (double *)calloc(NBin,sizeof(double));
457  double *RadProfS = (double *)calloc(NBin,sizeof(double));
458  double *Count = (double *)calloc(NBin,sizeof(double));
459  double *PlotMin = (double *)calloc(NGrid*NGrid,sizeof(double));
460  double *PlotThin = (double *)calloc(NGrid*NGrid,sizeof(double));
461  double *PlotTmp = (double *)calloc(NGrid*NGrid,sizeof(double));
462  double *AngList = (double *)calloc((NFile[1]-NFile[0]),sizeof(double));
463  double *HeiList = (double *)calloc((NFile[1]-NFile[0]),sizeof(double));
464  double *Weight = (double *)calloc(NWeight,sizeof(double));
465  int *WIndex = (int *)calloc(NWeight,sizeof(int));
466  double InvNBin = 1./(NBin);
467  Mat->FillWeightGauss(Weight,WIndex,NWeight,0.0,2.0);
468  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
469  for(int f=NFile[0];f<NFile[1];f++){
470  Processing(f);
471  if(Open(cFile[f],BF_PART)) return;
472  double Angle = Nano->Hamaker;
473  double Hei = Nano->Height;
474  for(int r=0;r<NBin;r++){
475  RadProf[r] = 0.;
476  Count[r] = 0.;
477  }
478  for(int p=0;p<pNPart();p++){
479  double Rad = sqrt(SQR(Pm[p].Pos[CLat1]-Nano->Pos[CLat1])+SQR(Pm[p].Pos[CLat2]-Nano->Pos[CLat2]));
480  //if(Rad < 1.5*Nano->Rad) continue;
481  int r = (int)(Rad*pInvEdge(3)*NBin);
482  if(r < 0 || r >= NBin) continue;
483  RadProf[r] += Pm[p].Pos[CNorm];
484  Count[r] += 1.;
485  }
486  for(int r=0;r<NBin;r++){
487  if(Count[r] <= 0.) continue;
488  RadProf[r] /= Count[r];
489  }
490  //smooth
491  Matrice Mask1(5);
492  Mask1.FillGaussian(.5,3.);
493  int NDim = 1;
494  int IfMinImConv = 0;
495  Mask1.ConvoluteMatrix(RadProf,NBin,NDim,IfMinImConv);
496  // for(int v=0;v<NBin;v++) RadProfS[v] = RadProf[v];
497  // InterBSpline1D(RadProfS,RadProf,NBin,NBin);
498  // for(int v=0;v<NBin;v++) RadProfS[v] = RadProf[v];
499  // InterBSpline1D(RadProfS,RadProf,NBin,NBin);
500  // for(int v=0;v<NBin;v++) RadProfS[v] = RadProf[v];
501  // InterBSpline1D(RadProfS,RadProf,NBin,NBin);
502  // Mat->ConvWeight(RadProf,NBin,Weight,WIndex,NWeight);
503  double Thinning = 100000.;
504  double MinDist = 0.;
505  int rLower = 0;
506  for(int r=3;r<NBin;r++){
507  if(Count[r] <= 0.) continue;
508  if(Thinning > RadProf[r]){
509  Thinning = RadProf[r];
510  MinDist = r*pEdge(3)/(double)NBin;
511  rLower = r;
512  }
513  }
514  //printf("%lf %lf %lf %lf\n",Angle,Hei,MinDist,Thinning);
515  AngList[f] = Angle;
516  HeiList[f] = Hei;
517  PlotMin[f] = MinDist;
518  PlotThin[f] = Thinning;
519  sprintf(FName,"ThinProfHei%.1fAng%02d.dat",Hei,(int)Angle);
520  FILE *FProf = fopen(FName,"w");
521  for(int r=0;r<NBin;r++){
522  if(Count[r] <= 0.) continue;
523  double Rad = r*pEdge(3)*InvNBin;
524  fprintf(FProf,"%lf %lf\n",Rad,RadProf[r]);
525  }
526  fclose(FProf);
527  }
528  //smooth and write dens
529  Matrice Mask(5,5);
530  Mask.FillGaussian(.5,3.);
531  Mask.Print();
532  int NDim = 2;
533  int IfMinImConv = 1;
534  for(int t=0;t<0;t++){
535  Mask.ConvoluteMatrix(PlotMin,NGrid,NDim,IfMinImConv);
536  Mask.ConvoluteMatrix(PlotThin,NGrid,NDim,IfMinImConv);
537  }
538  for(int s=0;s<NGrid*NGrid;s++) PlotTmp[s] = PlotMin[s];
539  InterBSpline2D(PlotTmp,PlotMin,NGrid,NGrid);
540  for(int s=0;s<NGrid*NGrid;s++) PlotTmp[s] = PlotThin[s];
541  InterBSpline2D(PlotTmp,PlotThin,NGrid,NGrid);
542  FILE *FMinDist = fopen("HeiAngMinDist.dat","w");
543  FILE *FThin = fopen("HeiAngThin.dat","w");
544  for(int f=NFile[0];f<NFile[1];f++){
545  double Angle = AngList[f];
546  double Hei = HeiList[f];
547  fprintf(FMinDist,"%lf %lf %lf\n",Hei,Angle,PlotMin[f]);
548  fprintf(FThin,"%lf %lf %lf\n",Hei,Angle,PlotThin[f]);
549  }
550  free(Count);
551  free(RadProf);
552  free(RadProfS);
553  free(PlotMin);
554  free(PlotThin);
555  fclose(FThin);
556  fclose(FMinDist);
557  printf("\n");
558 }
559 #ifdef USE_CGAL
560 void ElPoly::AreaDistrF(int NBin){
561  SetEdge(.5*MIN(pEdge(CLat1),pEdge(CLat2)),3);
562  double *Distr = (double *)calloc(NBin,sizeof(double));
563  double *DistrRadAv = (double *)calloc(NBin,sizeof(double));
564  double *DistrRad = (double *)calloc(NBin,sizeof(double));
565  double *VolContr = (double *)calloc(NBin,sizeof(double));
566  VolumeCircSlab(VolContr,NBin);
567  double FNormaInv = 1./(double)(NFile[1]-NFile[0]);
568  double BoxRadInv = .5*sqrt( SQR(pEdge(CLat1)) + SQR(pEdge(CLat2)) )/(double)NBin;
569  double AreaMean = 3.*.5*pNChain()/(pEdge(CLat1)*pEdge(CLat2));
570  for(int f=NFile[0];f<NFile[1];f++){
571  Processing(f);
572  if(OpenRisk(cFile[f],BF_CHAIN))return ;
573  BfDefChain();
574  AreaDistr(Distr,DistrRad,NBin);
575  for(int v=0;v<NBin;v++)
576  DistrRadAv[v] += DistrRad[v];
577  }
578  printf("\n");
579  //---------normalize------------------
580  double NormF = NFile[1] - NFile[0];
581  for(int v=0;v<NBin;v++){
582  DistrRadAv[v] /= VolContr[v]*NormF;
583  Distr[v] /= VolContr[v]*NormF;
584  }
585  //---------write----------------------
586  char FileName[120];
587  sprintf(FileName,"NanoAreaRadDistrR%.1fH%.1fh%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
588  FILE *File2Write = fopen(FileName,"w");
589  for(int v=0;v<NBin;v++){
590  fprintf(File2Write,"%lf %lf\n",v*BoxRadInv,DistrRadAv[v]*FNormaInv);
591  }
592  fclose(File2Write);
593  sprintf(FileName,"NanoAreaDistrR%.1fH%.1fh%.1f.dat",Nano->Rad,Nano->Hamaker,Nano->Height);
594  File2Write = fopen(FileName,"w");
595  for(int v=0;v<NBin;v++){
596  fprintf(File2Write,"%lf %lf\n",v*AreaMean/(double)NBin,Distr[v]*FNormaInv);
597  }
598  fclose(File2Write);
599  free(Distr);
600  free(DistrRadAv);
601  free(DistrRad);
602 }
603 #else
604 void ElPoly::AreaDistrF(int NBin){
605  printf("CGAL libraries not implemented\n");
606 }
607 #endif
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
void SysDef(char *cSystem)
Print a string with the system definitions.
Definition: VarData.cpp:167
bool OpenRisk(char *InFile, int BF)
Opens a file without reallocationg.
Definition: VarData.cpp:126
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
void AreaDistr(double *Distr, double *RadDistr, int NSample)
Calculate the (temporal/radial) area distribution.
Definition: VarDataCGAL.cpp:62
void FillWeightGauss(double *st, int *WIndex, int NWeight, double CutOff, double Sigma)
Fill the weight array with a gaussian fuction.
NANO * Nano
Extra particle.
Definition: VarData.h:1044
double Hamaker
Strength of the interaction.
Definition: VarData.h:447
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file".
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
int BfDefChain()
Definition of the chain.
double Rad
Size.
Definition: VarData.h:445
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
void AreaDistrF(int NBin)
Distribution of the areas around the protein.
double Height
Height of the cylinder.
Definition: VarData.h:449
void StalkLineProfF(int NBin)
Write the line describing a linear stalk.
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
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
double * x
Where the data are stored.
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.
Information of every chain.
Definition: VarData.h:236
int WormF(int Partition, int NBin)
Density profile of a non straight worm shape membrane worm.
void RadNormPos(int NBin, int NGrid)
Radial profile of the normal position of the particles.
int Angle(int NBin)
Calculation of the contact angle angle Boh?
Definition: ElPolyEl.cpp:78
double Pos[3]
Position.
Definition: VarData.h:427
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
int pNPCh()
Number of particle per chain.
int pNBlock()
Number of blocks.
int Worm(int Partition, int NSample, double *Border, double *dPoint)
Density profile along a worm like micelle.
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
void Stalk(int NSample, int NLevel, double **Plot, double Threshold)
Following the contour of a stalk.
double Pos[4]
xyzr Postion of the chain
Definition: VarData.h:238
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
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
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 NFile[2]
First and last file of the list.
Definition: ElPoly.h:418
double Angle(Vettore *u, Vettore *v)
Computes the angle between two Vetttore.
int CLat1
lateral coordinate
Definition: VarData.h:1076
void Prova()
Dummy function to operate on a sequence of files.
int Temperature(int NBin, int Coord)
Spatial distribution of the temperature.
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int InterBSpline2D(double **PlIn, double **PmOut, int NIn, int NOut)
2-d BSpline
void StalkF(int NSample)
Separate the leaflets.
void Print()
Print the entries.
int NanoParticle(int NBin)
Diffusion and density profile around the nanoparticle nano.
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
void StalkLineProf(double *Line, int NBin)
Describe the line for a linear stalk.
int pNPart()
Number of particle.