Allink  v0.1
VarDataBackFold.cpp
1 /***********************************************************************
2 VarData: This Program reads and writes a specific file format
3 storing all the information relative to a set of equal structure
4 polymers to the CHAIN, PART and GENERAL structures. It provides
5 two different ways to backfold the coordinates, a fanction that
6 creates an initial system with different option and some function
7 for the data analisys. The first calculate the distribution of the
8 monomer in the box, the second the distribution of the bonds.
9 Copyright (C) 2008 by Giovanni Marelli <sabeiro@virgilio.it>
10 
11 
12 This program is free software; you can redistribute it and/or modify
13 it under the terms of the GNU General Public License as published by
14 the Free Software Foundation; either version 2 of the License, or
15 (at your option) any later version.
16 
17 This program is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 GNU General Public License for more details.
21 
22 You should have received a copy of the GNU General Public License
23 along with this program; if not, write to the Free Software
24 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 ***********************************************************************/
26 #include "../include/VarData.h"
27 ;
29  if(VAR_IF_TYPE(SysType,VAR_CHAIN_DEF)){
30  for(int c=0,b=0;c<Gen->NChain;c++){
31  for(int d=0;d<3;d++){
32  Ch[c].Pos[d] -= floor(Ch[c].Pos[d]*pInvEdge(d))*Gen->Edge[d];
33  }
34  }
35  return 0;
36  }
37  VarMessage("BackFold.DefChains");
38  Vettore Ax0(1.,0.,0.);
39  Vettore Ax2(0.,0.,1.);
40  for(int b=0,NCh=0;b<pNBlock();NCh+=Block[b++].NChain){
41  double *xc = (double *)calloc(pNPCh(b),sizeof(double));
42  double *yc = (double *)calloc(pNPCh(b),sizeof(double));
43  double *zc = (double *)calloc(pNPCh(b),sizeof(double));
44  //if(strcasestr(Block[b].Name, "PEP") == Block[b].Name) continue;
45  int ChType = CHAIN_POLY;
46  if(Block[b].Asym == 0) ChType = CHAIN_ADDED;
47  //printf("\n%s type %d asym %d ppc %d #chain %d offset %d %d\n",Block[b].Name,ChType,Block[b].Asym,pNPCh(b),Block[b].NChain,NCh,Block[b].NChain+NCh);
48  for(int c=NCh;c<NCh+Block[b].NChain;c++){
49  Ch[c].Type = 0;
50  if(Block[b].Asym == 0){
51  VAR_ADD_TYPE(Ch[c].Type,CHAIN_ADDED);
52  }
53  for(int d=0;d<3;d++){
54  Ch[c].Pos[d] = 0.;
55  Ch[c].Vel[d] = 0.;
56  }
57  int p1 = Block[b].InitIdx + (c-NCh)*pNPCh(b);
58  for(int p=p1;p<p1+pNPCh(b);p++){
59  for(int d=0;d<3;d++){
60  Ch[c].Pos[d] += Pm[p].Pos[d];
61  Ch[c].Vel[d] += Pm[p].Vel[d];
62  }
63  }
64  for(int d=0;d<3;d++){
65  Ch[c].Pos[d] /= (double)pNPCh(b);
66  Ch[c].Vel[d] /= (double)pNPCh(b);
67  }
68  for(int ppc=0;ppc<pNPCh(b);ppc++){
69  xc[ppc] = Pm[p1+ppc].Pos[0];
70  yc[ppc] = Pm[p1+ppc].Pos[1];
71  zc[ppc] = Pm[p1+ppc].Pos[2];
72  }
73  RETTA rzx = Mat->InterRett(zc,xc,pNPCh(b));
74  RETTA rzy = Mat->InterRett(zc,yc,pNPCh(b));
75  double x1 = Pm[p1].Pos[2]*rzx.m + rzx.q;
76  double x2 = Pm[p1+pNPCh(b)-1].Pos[2]*rzx.m + rzx.q;
77  double y1 = Pm[p1].Pos[2]*rzy.m + rzy.q;
78  double y2 = Pm[p1+pNPCh(b)-1].Pos[2]*rzy.m + rzy.q;
79  Ch[c].Dir[2] = Pm[p1].Pos[2] - Pm[p1+pNPCh(b)-1].Pos[2];
80  Ch[c].Dir[1] = y1 - y2;
81  Ch[c].Dir[0] = x1 - x2;
82  Vettore ChDir(Ch[c].Dir[0],Ch[c].Dir[1],Ch[c].Dir[2]);
83  Ch[c].Angle = Ax0.Angle(&Ax2,&ChDir);
84  if(Ch[c].Angle > .5*M_PI)
85  VAR_ADD_TYPE(Ch[c].Type,CHAIN_UP);
86  else
87  VAR_ADD_TYPE(Ch[c].Type,CHAIN_DOWN);
88  Vettore Origin(pCm(0)-Ch[c].Pos[0],pCm(1)-Ch[c].Pos[1],pCm(2)-Ch[c].Pos[2]);
89  double Angle = Origin.Angle(&ChDir);
90  if(Angle < .5*M_PI)
91  VAR_ADD_TYPE(Ch[c].Type,CHAIN_OUTER);
92  else
93  VAR_ADD_TYPE(Ch[c].Type,CHAIN_INNER);
94  }
95  free(xc);
96  free(yc);
97  free(zc);
98  }
99  VAR_ADD_TYPE(SysType,VAR_CHAIN_DEF);
100  return 0;
101 }
103  int nBlock = 0;
104  for(int n=0;n<pNNano();n++){
105  if(Nano[n].Shape != SHAPE_CLUSTER) continue;
106  int nChain = 0;
107  for(int b=0;b<Gen->NBlock;b++){
108  if(!strncmp(Block[b].Name,"PEP",3)){
109  if(b > nBlock){
110  nBlock = b;
111  break;
112  }
113  }
114  nChain += Block[b].NChain;
115  }
116  int p1 = Block[nBlock].InitIdx;
117  double Cm[3] = {0.,0.,0.};
118  double NCm = 0.;
119  double *xc = (double *)calloc(pNPCh(nBlock),sizeof(double));
120  double *yc = (double *)calloc(pNPCh(nBlock),sizeof(double));
121  double *zc = (double *)calloc(pNPCh(nBlock),sizeof(double));
122  for(int ppc=0;ppc<pNPCh(nBlock);ppc++){
123  xc[ppc] = Pm[p1+ppc].Pos[0];
124  yc[ppc] = Pm[p1+ppc].Pos[1];
125  zc[ppc] = Pm[p1+ppc].Pos[2];
126  }
127  RETTA rzx = Mat->InterRett(zc,xc,pNPCh(nBlock));
128  RETTA rzy = Mat->InterRett(zc,yc,pNPCh(nBlock));
129  double x1 = Pm[p1].Pos[2]*rzx.m + rzx.q;
130  double x2 = Pm[p1+pNPCh(nBlock)-1].Pos[2]*rzx.m + rzx.q;
131  double y1 = Pm[p1].Pos[2]*rzy.m + rzy.q;
132  double y2 = Pm[p1+pNPCh(nBlock)-1].Pos[2]*rzy.m + rzy.q;
133  Nano[n].Axis[2] = Pm[p1].Pos[2] - Pm[p1+pNPCh(nBlock)-1].Pos[2];
134  Nano[n].Axis[1] = y1 - y2;
135  Nano[n].Axis[0] = x1 - x2;
136  for(int p=Block[nBlock].InitIdx;p<Block[nBlock].EndIdx;p++){
137  for(int d=0;d<3;d++) Cm[d] += Pm[p].Pos[d];
138  NCm += 1.;
139  }
140  for(int d=0;d<3;d++) Nano[n].Pos[d] = Cm[d]/NCm;
141  SetNanoBkf(n);
142  double Norm = 0.;
143  for(int d=0;d<3;d++){
144  //SigErr(isnan(Nano[n].Axis[d]),"Wrong nano axis %d %lf ",d,Nano[n].Axis[d]);
145  Norm += SQR(Nano[n].Axis[d]);
146  }
147  Norm = sqrt(Norm);
148  for(int d=0;d<3;d++){
149  Nano[n].Axis[d] /= -Norm;
150  }
151  free(xc);
152  free(yc);
153  free(zc);
154  }
155 }
157  double Edge[3][2];
158  if( VAR_IF_TYPE(SysType,VAR_EDGE) ) return 0;
159  Edge[0][0] = Pm[0].Pos[0];Edge[0][1] = Pm[0].Pos[0];
160  Edge[1][0] = Pm[0].Pos[1];Edge[1][1] = Pm[0].Pos[1];
161  Edge[2][0] = Pm[0].Pos[2];Edge[2][1] = Pm[0].Pos[2];
162  for(int p=0;p<Gen->NPart;p++){
163  for(int d=0;d<3;d++){
164  if(Edge[d][1] < Pm[p].Pos[d] )
165  Edge[d][1] = Pm[p].Pos[d];
166  if(Edge[d][0] > Pm[p].Pos[d])
167  Edge[d][0] = Pm[p].Pos[d];
168  }
169  }
170  for(int d=0;d<3;d++)
171  SetEdge(Edge[d][1] - Edge[d][0],d);
172  // SetEdge(Edge[d][1],d);
173  VAR_ADD_TYPE(SysType,VAR_EDGE);
174  return 0;
175 }
177  VarMessage("BackFold");
178  if(BackFold == BF_SKIP) return ;
179  BfEdge();
180  BfPep();
181  double Shift[3];
182  for(int d=0;d<3;d++){
183  Shift[d] = -ShiftPos[d]*Gen->Edge[d];
184  }
185  for(int c=0;c<Gen->NChain;c++){
186  for(int d=0;d<3;d++){
187  Ch[c].Pos[d] -= Shift[d];
188  //Ch[c].Pos[d] -= floor(Ch[c].Pos[d]/Gen->Edge[d])*Gen->Edge[d];
189  }
190  }
191  if(BackFold == BF_CHAIN){
192  BfDefChain();
193  for(int p=0;p<Gen->NPart;p++){
194  for(int d=0;d<3;d++){
195  Pm[p].Pos[d] -= Shift[d];
196  Pm[p].Bkf[d] = -floor((Ch[Pm[p].CId].Pos[d]-Shift[d])*pInvEdge(d))*Gen->Edge[d];
197  }
198  }
199  for(int c=0;c<Gen->NChain;c++){
200  for(int d=0;d<3;d++){
201  Ch[c].Pos[d] -= floor(Ch[c].Pos[d]*pInvEdge(d))*pEdge(d);
202  }
203  }
204  }
205  else if(BackFold == BF_PART){
206  for(int p=0;p<Gen->NPart;p++){
207  for(int d=0;d<3;d++){
208  Pm[p].Pos[d] -= Shift[d];
209  Pm[p].Bkf[d] = -floor(Pm[p].Pos[d]*pInvEdge(d))*Gen->Edge[d];
210  }
211  }
212  }
213  else if(BackFold == BF_NANO){
214  for(int d=0;d<3;d++){
215  ShiftPos[d] = pNanoPos(0,d)*pInvEdge(d);
216  Shift[d] = (ShiftPos[d]-.5)*Gen->Edge[d];
217  }
218  for(int p=0;p<Gen->NPart;p++){
219  for(int d=0;d<3;d++){
220  Pm[p].Pos[d] -= Shift[d];
221  Pm[p].Bkf[d] = -floor(Pm[p].Pos[d]*pInvEdge(d))*Gen->Edge[d];
222  }
223  }
224  }
225  else if(BackFold == BF_TILT){
226  int n=0;
227  for(int d=0;d<3;d++){
228  ShiftPos[d] = pNanoPos(n,d);
229  Shift[d] = ShiftPos[d]-.5*Gen->Edge[d];
230  }
231  Vettore NanoAx(Nano[n].Axis[0],Nano[n].Axis[1],Nano[n].Axis[2]);
232  NanoAx.Set(0.,CNorm);
233  NanoAx.Normalize();
234  Vettore Ex(0.,0.,0.);
235  Vettore Zed(0.,0.,0.);
236  Zed.Set(1.,CNorm);
237  Ex.Set(1.,CLat1);
238  //double Angle = atan(Nano[0].Axis[CLat2]/Nano[0].Axis[CLat1]);//Zed.Angle(&NanoAx,&Ex);
239  double Angle = Ex.Angle(&NanoAx);
240  if(isnan(Angle)) Angle = 0.;
241  int NRow = 4;
242  Matrice M(Zed.x,.5*Angle,NRow);
243  // printf("angle %lf\n",Angle*360./DUE_PI);
244  // M.Print();
245  Vettore PosOld(NRow);
246  Vettore PosNew(NRow);
247  for(int p=0;p<pNPart();p++){
248  for(int d=0;d<3;d++){
249  PosOld.Set(Pm[p].Pos[d]-pNanoPos(0,d),d);
250  }
251  for(int r=0;r<NRow;r++){
252  double Temp=0.;
253  for(int c=0;c<NRow;c++)
254  Temp += M.Val(r,c)*PosOld.Val(c);
255  PosNew.Set(Temp,r);
256  }
257  //PosNew = M.Mult(PosOld);
258  for(int d=0;d<3;d++){
259  Pm[p].Pos[d] = PosNew.Val(d) +.5*pEdge(d);
260  Pm[p].Bkf[d] = -floor(Pm[p].Pos[d]*pInvEdge(d))*Gen->Edge[d];
261  }
262  }
263  Nano[0].Axis[CLat2] = -sqrt(1.-SQR(Nano[0].Axis[CNorm]));
264  Nano[0].Axis[CLat1] = 0.;
265  }
266  for(int n=0;n<MAX(1,Gen->NNano);n++){
267  for(int d=0;d<3;d++){
268  Nano[n].Pos[d] -= Shift[d];
269  }
270  SetNanoBkf(n);
271  }
272 }
273 //obsolete?
274 void VarData::DefBlock(int *NChStep,int How){
275  if(How == VAR_OPPOSED){
276  Gen->NBlock = 4;
277  Block = (BLOCK *)realloc(Block,Gen->NBlock*sizeof(*Block));
278  sprintf(Block[0].Name,"Lower1");
279  sprintf(Block[1].Name,"Upper1");
280  sprintf(Block[2].Name,"Lower2");
281  sprintf(Block[3].Name,"Upper2");
282  }
283  else if(How == VAR_VESICLE || How == VAR_TUBE){
284  Gen->NBlock = 2;
285  Block = (BLOCK *)realloc(Block,Gen->NBlock*sizeof(*Block));
286  sprintf(Block[0].Name,"Inner");
287  sprintf(Block[1].Name,"Outer");
288  }
289  Block[0].InitIdx= 0;
290  Block[0].EndIdx = NChStep[0]*Gen->NPCh;
291  Block[0].NChain = NChStep[0];
292  Block[0].NPCh = Gen->NPCh;
293  for(int b=1;b<Gen->NBlock;b++){
294  Block[b].InitIdx= Block[b-1].EndIdx;
295  Block[b].EndIdx = NChStep[b]*Gen->NPCh+Block[b-1].EndIdx;
296  Block[b].NChain = NChStep[b];
297  Block[b].NPCh = Gen->NPCh;
298  }
299  for(int b=0;b<Gen->NBlock;b++)
300  printf("%d %d %d %d\n",Block[b].NChain,Block[b].InitIdx,Block[b].EndIdx,Block[b].NChain);
301 }
302 //obsolete?
303 bool VarData::BackFold(int How){
304  VarMessage("BackFold");
305  if(How == BF_SKIP) return 0;
306  BfEdge();
307  double ShiftNano[3];
308  SetNanoBkf(0);
309  if(How == BF_PART)
310  for(int p=0;p<Gen->NPart;p++)
311  for(int d=0;d<3;d++)
312  Pm[p].Bkf[d] = -floor(Pm[p].Pos[d]/Gen->Edge[d])*Gen->Edge[d];
313  else if(How == BF_CHAIN)
314  for(int p=0;p<Gen->NPart;p++){
315  int c = Pm[p].CId;
316  if(c < 0 || c >= Gen->NChain) continue;
317  for(int d=0;d<3;d++)
318  Pm[p].Bkf[d] = -floor(Ch[c].Pos[d]/Gen->Edge[d])*Gen->Edge[d];
319  }
320  else if(How == BF_NANO)
321  for(int p=0;p<Gen->NPart;p++)
322  for(int d=0;d<3;d++)
323  Pm[p].Pos[d] += ShiftNano[d];
324  BfDefChain();
325  return 0;
326 }
327 //Obsolete
328 bool VarData::ShiftSys(int How){
329  if(How == SHIFT_NO) return 0;
330  double Shift[3] = {0.,0.,0.};
331  SetNanoBkf(0);
332  if(How == SHIFT_CM){//Cm
333  for(int d=0;d<3;d++)
334  Gen->Cm[d] = 0.;
335  for(int p=0;p<Gen->NPart;p++){
336  Gen->Cm[0] += Pm[p].Pos[0];
337  Gen->Cm[1] += Pm[p].Pos[1];
338  Gen->Cm[2] += Pm[p].Pos[2];
339  }
340  for(int d=0;d<3;d++){
341  Gen->Cm[d] /= (double)Gen->NPart;
342  Shift[d] = (Gen->Edge[d]*.5-Gen->Cm[d]);
343  Nano->Pos[d] += Shift[d];
344  }
345  SetNanoBkf(0);
346  }
347  else if(How == SHIFT_NANO){//Nano Particle
348  for(int d=0;d<3;d++){
349  Shift[d] = (Gen->Edge[d]*.5 - pNanoPos(0,d));
350  Nano->Pos[d] += Shift[d];
351  }
352  SetNanoBkf(0);
353  }
354  else if(How == SHIFT_CM_NANO){//x,y Nano; z Cm
355  for(int d=0;d<3;d++)
356  Gen->Cm[d] = 0.;
357  for(int p=0;p<Gen->NPart;p++){
358  Gen->Cm[2] += Pm[p].Pos[2];
359  }
360  Gen->Cm[0] = pNanoPos(0,0);
361  Gen->Cm[1] = pNanoPos(0,1);
362  Gen->Cm[2] /= (double)Gen->NPart;
363  for(int d=0;d<3;d++){
364  Shift[d] = (Gen->Edge[d]*.5 - Gen->Cm[d]);
365  if(d < 2) Nano->Pos[d] += Shift[d];
366  }
367  SetNanoBkf(0);
368  }
369  else return 1;
370  for(int p=0;p<Gen->NPart;p++){
371  // if(Pm[p].Typ == 2){
372  // Pm[p].Pos[0] = Nano->PosBf[0];
373  // Pm[p].Pos[1] = Nano->PosBf[1];
374  // Pm[p].Pos[2] = Nano->PosBf[2];
375  // continue;
376  // }
377  Pm[p].Pos[0] += Shift[0];
378  Pm[p].Pos[1] += Shift[1];
379  Pm[p].Pos[2] += Shift[2];
380  }
381  //printf("NanoPos %lf %lf %lf %lf %lf %lf\n",Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],Nano->PosBf[0],Nano->PosBf[1],Nano->PosBf[2]);
382  return 0;
383 }
384 void VarData::BackBone(double *Line,int NBin){
385  double *Count = new double[NBin];
386  double *LineS = new double[NBin];
387  for(int b=0;b<NBin;b++){
388  Line[b] = 0.;
389  Count[b] = 0.;
390  }
391  CLat1 = 1;
392  CLat2 = 0;
393  //average
394  for(int p=0;p<pNPart();p++){
395  if(Pm[p].Pos[2] < .02*pEdge(2))continue;
396  int b = (int)(Pm[p].Pos[CLat1]*pInvEdge(CLat1)*NBin);
397  if(b<0 || b >= NBin)continue;
398  double Weight = Pm[p].Pos[2];
399  Line[b] += Pm[p].Pos[CLat2]*Weight;
400  Count[b] += Weight;
401  }
402  for(int b=0;b<NBin;b++){
403  Line[b] /= Count[b] > 0. ? Count[b] : 1.;
404  }
405  //smooth
406  for(int v=0;v<NBin;v++) LineS[v] = Line[v];
407  InterBSpline1D(LineS,Line,NBin,NBin);
408  for(int v=0;v<NBin;v++) LineS[v] = Line[v];
409  InterBSpline1D(LineS,Line,NBin,NBin);
410  // //reweighting
411  // for(int b=0;b<NBin;b++){
412  // Line[b] = 0.;
413  // Count[b] = 0.;
414  // }
415  // double NormDist = 1./(.5*pEdge(CLat2));
416  // for(int p=0;p<pNPart();p++){
417  // if(Pm[p].Pos[2] < .02*pEdge(2))continue;
418  // int b = (Pm[p].Pos[CLat1]*pInvEdge(CLat1)*NBin);
419  // if(b < 0 || b >= NBin)continue;
420  // double Dist = SQR(Pm[p].Pos[CLat2] - LineS[b]);
421  // if(Dist > SQR(.1)) continue;
422  // int bm1 = b-3;
423  // if(bm1 < 0) bm1 = b;
424  // double Distm1 = SQR(Pm[p].Pos[CLat2] - LineS[bm1]);
425  // if(Distm1 > SQR(.3)) continue;
426  // int bp1 = b+3;
427  // if(bp1 > NBin-1) bp1 = b;
428  // double Distp1 = SQR(Pm[p].Pos[CLat2] - LineS[bp1]);
429  // if(Distp1 > SQR(.3)) continue;
430  // double Weight = 100.*(1. - Dist*Distm1*Distp1*CUBE(NormDist));
431  // Line[b] += Pm[p].Pos[CLat2]*Weight;
432  // Count[b] += Weight;
433  // }
434  // for(int b=0;b<NBin;b++){
435  // Line[b] /= Count[b] > 0. ? Count[b] : 1.;
436  // }
437  // //smooth
438  // for(int v=0;v<NBin;v++) LineS[v] = Line[v];
439  // InterBSpline1D(LineS,Line,NBin,NBin);
440  // for(int v=0;v<NBin;v++) LineS[v] = Line[v];
441  // InterBSpline1D(LineS,Line,NBin,NBin);
442  //fourier
443 #ifdef USE_FFTW
444  fftw_complex *FouOut = (fftw_complex *)fftw_malloc(NBin*sizeof(fftw_complex));
445  fftw_complex *FouIn = (fftw_complex *)fftw_malloc(NBin*sizeof(fftw_complex));
446  fftw_plan direct = fftw_plan_dft_1d(NBin,
447  FouIn,FouOut,FFTW_FORWARD,FFTW_PATIENT);
448  fftw_plan reverse = fftw_plan_dft_1d(NBin,
449  FouOut,FouIn,FFTW_BACKWARD,FFTW_PATIENT);
450  for(int b=0;b<NBin;b++) FouIn[b][0] = Line[b];
451  fftw_execute(direct);
452  int NComp = NBin;// - 2;//NBin/2;
453  for(int b=NBin-1;b>=NComp;b--){
454  FouOut[b][0] = 0.;
455  FouOut[b][1] = 0.;
456  }
457  fftw_execute(reverse);
458  for(int b=0;b<NBin;b++) Line[b] = FouIn[b][0]/(double)NBin;
459 #endif
460  //write
461  if(1==0){
462  FILE *FOut = fopen("BackBone.dat","w");
463  fprintf(FOut,"#l(%lf %lf %lf) v[%d] d[part]\n",pEdge(CLat1),pEdge(CLat2),pEdge(CNorm),NBin);
464  int NPart = 0;
465  for(int p=0;p<pNPart();p++){
466  if(Pm[p].Pos[2] < .02*pEdge(2)) continue;
467  fprintf(FOut,"{t[0 0 0] x(%lf %lf %lf)}\n",Pm[p].Pos[CLat1],Pm[p].Pos[CLat2],Pm[p].Pos[CNorm]);
468  NPart++;
469  }
470  for(int b=0;b<NBin;b++){
471  double x = (b+.5) * pEdge(CLat1)/(double)NBin;
472  double y = b * pEdge(CLat2)/(double)NBin;
473  fprintf(FOut,"{t[0 0 1] x(%lf %lf 0.) l[%d]} \n",x,Line[b],NPart+1);
474  NPart++;
475  }
476  fclose(FOut);
477  }
478  FILE *FOut1 = fopen("LineProf.dat","w");
479  for(int b=0;b<NBin;b++){
480  double x = (b+.5) * pEdge(CLat1)/(double)NBin;
481  double y = b * pEdge(CLat2)/(double)NBin;
482  fprintf(FOut1,"%lf %lf\n",x,Line[b]);
483  }
484  fclose(FOut1);
485 }
486 void VarData::StalkLineProf(double *Line,int NBin){
487  //-------------------alloc
488  Vettore Ax0(1.,0.,0.);
489  Vettore Ax2(0.,0.,1.);
490  double Dist[4];
491  BfDefChain();
492  //double CmStalk[3] = {0.,0.,0.};
493  double *CountStalk = (double *)calloc(NBin,sizeof(double));
494  int NSample = 8;
495  int NIter = 6;
496  Matrice Mask(5);
497  Mask.FillGaussian(.5,3.);
498  int NDim = 1;
499  int IfMinImConv = 1;
500  int IfSpline = 1;
501  int IfGaussian = 1;
502  double InvNSample = 1./(double)NSample;
503  double **PlotUp = (double **)calloc(NSample,sizeof(double));
504  double **PlotUpS = (double **)calloc(NSample,sizeof(double));
505  double **CountUp = (double **)calloc(NSample,sizeof(double));
506  double **PlotDown = (double **)calloc(NSample,sizeof(double));
507  double **PlotDownS = (double **)calloc(NSample,sizeof(double));
508  double **CountDown = (double **)calloc(NSample,sizeof(double));
509  double *LineS = (double *)calloc(NBin,sizeof(double));
510  for(int s=0;s<NSample;s++){
511  PlotUp[s] = (double *)calloc(NSample,sizeof(double));
512  PlotUpS[s] = (double *)calloc(NSample,sizeof(double));
513  CountUp[s] = (double *)calloc(NSample,sizeof(double));
514  PlotDown[s] = (double *)calloc(NSample,sizeof(double));
515  PlotDownS[s] = (double *)calloc(NSample,sizeof(double));
516  CountDown[s] = (double *)calloc(NSample,sizeof(double));
517  }
518  for(int vx=0;vx<NBin;vx++){
519  Line[vx] = 0.;
520  }
521  //---------------defines the two midplanes
522  //center of mass of the upper layer of the lower membrane
523  for(int b=0,cOff=0;b<pNBlock();b++,cOff+=Block[b].NChain){
524  if(strcasestr(Block[b].Name, "TT0") != Block[b].Name) continue;
525  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
526  if(!VAR_IF_TYPE(Ch[Pm[p].CId].Type,CHAIN_UP))continue;
527  if(Pm[p].Typ == 0) continue;
528  double Posx = Pm[p].Pos[CLat1] - floor(Pm[p].Pos[CLat1]*pInvEdge(CLat1))*pEdge(CLat1);
529  double Posy = Pm[p].Pos[CLat2] - floor(Pm[p].Pos[CLat2]*pInvEdge(CLat2))*pEdge(CLat2);
530  int sx = (int)(Posx*pInvEdge(CLat1)*NSample);
531  int sy = (int)(Posy*pInvEdge(CLat2)*NSample);
532  PlotDown[sx][sy] += Pm[p].Pos[CNorm];
533  CountDown[sx][sy] += 1.;
534  }
535  }
536  //center of mass of the lower layer of the upper membrane
537  for(int b=0,cOff=0;b<pNBlock();b++,cOff+=Block[b].NChain){
538  if(strcasestr(Block[b].Name, "TT1") != Block[b].Name) continue;
539  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
540  if(!VAR_IF_TYPE(Ch[Pm[p].CId].Type,CHAIN_DOWN))continue;
541  if(Pm[p].Typ == 0) continue;
542  //if(Pm[p].Pos[2] < .5*pEdge(2))continue;
543  double Posx = Pm[p].Pos[CLat1] - floor(Pm[p].Pos[CLat1]*pInvEdge(CLat1))*pEdge(CLat1);
544  double Posy = Pm[p].Pos[CLat2] - floor(Pm[p].Pos[CLat2]*pInvEdge(CLat2))*pEdge(CLat2);
545  int sx = (int)(Posx*pInvEdge(CLat1)*NSample);
546  int sy = (int)(Posy*pInvEdge(CLat2)*NSample);
547  PlotUp[sx][sy] += Pm[p].Pos[CNorm];
548  CountUp[sx][sy] += 1.;
549  }
550  }
551  for(int sx=0;sx<NSample;sx++){
552  for(int sy=0;sy<NSample;sy++){
553  PlotUp[sx][sy] /= CountUp[sx][sy] > 0. ? CountUp[sx][sy] : 1.;
554  PlotDown[sx][sy] /= CountDown[sx][sy] > 0. ? CountDown[sx][sy] : 1.;
555  }
556  }
557  //InterBSpline2D(PlotUp,PlotUpS,NSample,NSample);
558  //InterBSpline2D(PlotDown,PlotDownS,NSample,NSample);
559  for(int sx=0;sx<NSample;sx++){
560  for(int sy=0;sy<NSample;sy++){
561  PlotUpS[sx][sy] = PlotUp[sx][sy];
562  PlotDownS[sx][sy] = PlotDown[sx][sy];
563  }
564  }
565  //----------------count the chains between the two midplanes
566  double AngleMax = 1./(.5*M_PI);
567  for(int p=0;p<pNPart();p++){
568  if(Pm[p].Typ == 1) continue;
569  int sx = (int)(Pm[p].Pos[CLat1]*pInvEdge(CLat1)*NSample);
570  if(sx < 0 || sx >= NSample) continue;
571  int sy = (int)(Pm[p].Pos[CLat2]*pInvEdge(CLat2)*NSample);
572  if(sy < 0 || sy >= NSample) continue;
573  if(Pm[p].Pos[CNorm] > PlotUpS[sx][sy] || Pm[p].Pos[CNorm] < PlotDownS[sx][sy]) continue;
574  int vx = (int)(Pm[p].Pos[CLat2]*pInvEdge(CLat2)*NBin);
575  if(vx < 0 || vx >= NBin) continue;
576  int c = Pm[p].CId;
577  Vettore ChDir(Ch[c].Dir[0],Ch[c].Dir[1],Ch[c].Dir[2]);
578  double Angle = Ax0.Angle(&Ax2,&ChDir);
579  if(Angle < .7 || Angle > 2.3) continue;
580  double Weight = 1.;//1. - SQR((Angle-.5*M_PI)*AngleMax);
581  Line[vx] += Pm[p].Pos[CLat1]*Weight;
582  CountStalk[vx] += Weight;
583  }
584  for(int vx=0;vx<NBin;vx++){
585  Line[vx] /= CountStalk[vx] > 0. ? CountStalk[vx] : 1.;
586  }
587  //-------------------------smoothing
588  if(IfSpline){
589  for(int v=0;v<NBin;v++) LineS[v] = Line[v];
590  InterBSpline1D(LineS,Line,NBin,NBin);
591  for(int v=0;v<NBin;v++) LineS[v] = Line[v];
592  InterBSpline1D(LineS,Line,NBin,NBin);
593  }
594  if(IfGaussian){
595  Mask.ConvoluteMatrix(Line,NBin,NDim,IfMinImConv);
596  Mask.ConvoluteMatrix(Line,NBin,NDim,IfMinImConv);
597  }
598  // // //------------------------Squared-average------------
599  for(int i=0;i<NIter;i++){
600  for(int v=0;v<NBin;v++){
601  LineS[v] = Line[v];
602  Line[v] = 0.;
603  CountStalk[v] = 0.;
604  }
605  double NormDist = 1./(.5*pEdge(CLat2));
606  for(int p=0;p<pNPart();p++){
607  if(Pm[p].Typ == 1) continue;
608  int sx = (int)(Pm[p].Pos[CLat1]*pInvEdge(CLat1)*NSample);
609  if(sx < 0 || sx >= NSample) continue;
610  int sy = (int)(Pm[p].Pos[CLat2]*pInvEdge(CLat2)*NSample);
611  if(sy < 0 || sy >= NSample) continue;
612  if(Pm[p].Pos[CNorm] > PlotUpS[sx][sy] || Pm[p].Pos[CNorm] < PlotDownS[sx][sy]) continue;
613  int vx = (int)(Pm[p].Pos[CLat2]*pInvEdge(CLat2)*NBin);
614  if(vx < 0 || vx >= NBin) continue;
615  double Dist = SQR(Pm[p].Pos[CLat1] - LineS[vx]);
616  if(Dist > SQR(3.)) continue;
617  double Weight = 100.*(1. - SQR((Pm[p].Pos[CLat1]-LineS[vx])*NormDist));
618  Line[vx] += Pm[p].Pos[CLat1]*Weight;
619  CountStalk[vx] += Weight;
620  }
621  double Average = 0.;
622  for(int vx=0;vx<NBin;vx++){
623  Line[vx] /= CountStalk[vx] > 0. ? CountStalk[vx] : 1.;
624  Average += Line[vx];
625  }
626  // //-------------------------smoothing
627  if(IfSpline){
628  for(int v=0;v<NBin;v++) LineS[v] = Line[v];
629  InterBSpline1D(LineS,Line,NBin,NBin);
630  for(int v=0;v<NBin;v++) LineS[v] = Line[v];
631  InterBSpline1D(LineS,Line,NBin,NBin);
632  }
633  if(IfGaussian){
634  Mask.ConvoluteMatrix(Line,NBin,NDim,IfMinImConv);
635  Mask.ConvoluteMatrix(Line,NBin,NDim,IfMinImConv);
636  }
637  }
638  //-------------------------print-temp-files
639  if(1==0){
640  FILE *F2Write = fopen("TwoMid.dat","w");
641  FILE *FWrite = fopen("StalkLine.dat","w");
642  FILE *FWrite1 = fopen("StalkLine1.dat","w");
643  fprintf(F2Write,"#l(%lf %lf %lf) v[%d] d[part]\n",pEdge(0),pEdge(1),pEdge(2),NSample);
644  WriteSurf(F2Write,PlotUpS,NSample,0);
645  WriteSurf(F2Write,PlotDownS,NSample,SQR(NSample));
646  for(int p=0;p<pNPart();p++){
647  if(Pm[p].Typ == 1) continue;
648  int sx = (int)(Pm[p].Pos[CLat1]*pInvEdge(CLat1)*NSample);
649  if(sx < 0 || sx >= NSample) continue;
650  int sy = (int)(Pm[p].Pos[CLat2]*pInvEdge(CLat2)*NSample);
651  if(sy < 0 || sy >= NSample) continue;
652  if(Pm[p].Pos[CNorm] > PlotUpS[sx][sy] || Pm[p].Pos[CNorm] < PlotDownS[sx][sy]) continue;
653  int c = Pm[p].CId;
654  Vettore ChDir(Ch[c].Dir[0],Ch[c].Dir[1],Ch[c].Dir[2]);
655  double Angle = Ax0.Angle(&Ax2,&ChDir);
656  if(Angle < .7 || Angle > 2.3) continue;
657  fprintf(F2Write,"{t[0 0 1] x(%lf %lf %lf)}\n",Pm[p].Pos[0],Pm[p].Pos[1],Pm[p].Pos[2]);
658  fprintf(FWrite,"%lf %lf \n",Pm[p].Pos[1],Pm[p].Pos[0]);
659  }
660  for(int vx=0;vx<NBin;vx++){
661  fprintf(FWrite1,"%lf %lf \n",vx*pEdge(CLat2)/(double)NBin,Line[vx]);
662  }
663  fclose(FWrite);
664  fclose(F2Write);
665  //exit(0);
666  }
667  for(int s=0;s<NSample;s++){
668  free(PlotUp[s]);
669  free(PlotUpS[s]);
670  free(CountUp[s]);
671  free(PlotDown[s]);
672  free(PlotDownS[s]);
673  free(CountDown[s]);
674  }
675  free(PlotUp);
676  free(PlotUpS);
677  free(CountUp);
678  free(CountDown);
679  free(PlotDown);
680  free(PlotDownS);
681  free(CountStalk);
682 }
683 void VarData::StalkPos2(double *OldPos,double *CmStalk){
684  Vettore Ax0(1.,0.,0.);
685  Vettore Ax2(0.,0.,1.);
686  double Dist[4];
687  BfDefChain();
688  double CountStalk = 1.;
689  for(int b=0,cOff=0;b<pNBlock();b++,cOff+=Block[b].NChain){
690  if(strncmp(Block[b].Name,"TT",2)) continue;
691  for(int c=cOff;c<cOff+Block[b].NChain;c++){
692  Vettore ChDir(Ch[c].Dir[0],Ch[c].Dir[1],Ch[c].Dir[2]);
693  double Angle = Ax0.Angle(&Ax2,&ChDir);
694  if(Ch[c].Pos[2] < OldPos[2] - 1.5 || Ch[c].Pos[2] > OldPos[2] + 1.5) continue;
695  if(Angle < 1. || Angle > 2.) continue;
696  for(int p=0;p<Block[b].NPCh;p++){
697  int pCurr = Block[b].InitIdx + (c-cOff)*Block[b].NPCh + p;
698  if(Pm[p].Typ == 1) continue;
699  for(int d=0;d<3;d++){
700  Dist[d] = Pm[p].Pos[d] - OldPos[d];
701  Dist[d] -= floor(Dist[d]*pInvEdge(d))*pEdge(d);
702  }
703  Dist[3] = (SQR(Dist[0])+SQR(Dist[1])+SQR(Dist[2]));
704  double Weight = .00001/(Dist[3]*SQR(Angle-.5*M_PI));
705  //double Weight = .1/(Dist[3]);
706  if(Weight > 1.) continue;
707  //printf("%lf %lf %lf\n",Angle,Dist[3],Weight);
708  for(int d=0;d<3;d++){
709  CmStalk[d] += Pm[p].Pos[d]*Weight;//(Pm[p].Pos[d]-OldPos[d])*Weight;
710  }
711  CountStalk += Weight;
712  }
713  }
714  }
715  if(CountStalk <= 0.) CountStalk = 1.;
716  for(int d=0;d<3;d++){
717  CmStalk[d] /= CountStalk;
718  }
719 }
720 void VarData::StalkPos3(double *OldPos,double *CmStalk){
721  int NSample = 32;
722  int NGrid = NSample-1;
723  double VolEl = pVol()/(double)CUB(NSample);
724  double *Plot3d = (double *)calloc(CUBE(NSample),sizeof(double));
725  double *Plot2d = (double *)calloc(SQR(NGrid),sizeof(double));
726  double Min = 0.;
727  double Max = 0.;
728  double CountStalk = 0.;
729  VAR_TRIANGLE *Triang = NULL;
730  double Dist[4];
731  for(int d=0;d<3;d++) CmStalk[d] = 0.;
732  for(int p=0;p<pNPart();p++){
733  //if(Pm[p].Typ != 0) continue;
734  double Posx = Pm[p].Pos[0] - floor(Pm[p].Pos[0]*pInvEdge(0))*pEdge(0);
735  double Posy = Pm[p].Pos[1] - floor(Pm[p].Pos[1]*pInvEdge(1))*pEdge(1);
736  double Posz = Pm[p].Pos[2] - floor(Pm[p].Pos[2]*pInvEdge(2))*pEdge(2);
737  int sx = (int)(Posx*pInvEdge(0)*NSample);
738  int sy = (int)(Posy*pInvEdge(1)*NSample);
739  int sz = (int)(Posz*pInvEdge(2)*NSample);
740  int sTot = (sx*NSample+sy)*NSample+sz;
741  Plot3d[sTot] += VolEl;
742  CmStalk[CNorm] += Pm[p].Pos[CNorm];
743  if(Max < Plot3d[sTot]) Max = Plot3d[sTot];
744  if(Min > Plot3d[sTot]) Min = Plot3d[sTot];
745  }
746  double IsoLevel = .1*Max;
747  int NTri = 0;
748  Triang = MarchingCubes(Plot3d,NSample,IsoLevel,&NTri);
749  free(Plot3d);
750  for(int gx=0;gx<NGrid;gx++){
751  for(int gy=0;gy<NGrid;gy++){
752  Plot2d[gx*NGrid+gy] = 1.;
753  }
754  }
755  for(int t=0;t<NTri;t++){
756  for(int v=0;v<3;v++){
757  for(int d=0;d<3;d++){
758  Dist[d] = OldPos[d] - Triang[t].p[v].x[d];
759  }
760  Dist[3] = SQR(Dist[0])+SQR(Dist[1])+SQR(Dist[2]);
761  int gx = (int)(Triang[t].p[v].x[CLat1]*pInvEdge(CLat1)*NGrid);
762  int gy = (int)(Triang[t].p[v].x[CLat2]*pInvEdge(CLat2)*NGrid);
763  if(gx < 0 || gx >= NGrid) continue;
764  if(gy < 0 || gy >= NGrid) continue;
765  Plot2d[gx*NGrid+gy] += sqrt(Dist[3]);
766  }
767  }
768  double Dx = .5*pEdge(CLat1)/(double)NGrid;
769  double Dy = .5*pEdge(CLat2)/(double)NGrid;
770  double Count = 0.;
771  if(1==0){
772  FILE *Ciccia = fopen("Ciccia.dat","w");
773  for(int gx=0;gx<NGrid;gx++){
774  double x = gx/(double)NGrid*pEdge(CLat1) + Dx;
775  for(int gy=0;gy<NGrid;gy++){
776  double y = gy/(double)NGrid*pEdge(CLat1) + Dy;
777  double Weight = 1./(Plot2d[gx*NGrid+gy]);
778  fprintf(Ciccia,"%lf %lf %lf\n",x,y,Weight*10000.);
779  }
780  }
781  fclose(Ciccia);
782  exit(0);
783  }
784  for(int gx=0;gx<NGrid;gx++){
785  double x = gx/(double)NGrid*pEdge(CLat1) + Dx;
786  for(int gy=0;gy<NGrid;gy++){
787  double y = gy/(double)NGrid*pEdge(CLat1) + Dy;
788  double Weight = 1./(Plot2d[gx*NGrid+gy]*SQR(x-OldPos[0])*SQR(y-OldPos[1]));
789  CmStalk[CLat1] += x*Weight;
790  CmStalk[CLat2] += y*Weight;
791  CountStalk += Weight;
792  //printf(" %lf %lf %lf\n",CmStalk[CLat1],CmStalk[CLat2],Weight);
793  }
794  }
795  CmStalk[0] /= CountStalk;
796  CmStalk[1] /= CountStalk;
797  CmStalk[2] /= pNPart();
798  free(Plot2d);
799 }
800 int VarData::StalkPos4(double *OldPos,double *CmStalk){
801  Nano->Shape = SHAPE_TORUS;
803  int NInt = 30;
804  int NBin = 12;
805  double kElPhob = 40000.;
806  double kElPhil = -1.;//-2.;
807  double MinRad = .2;//.4;
808  double MinHei = .1;//2.;
809  double RadMax = .3;
810  double RadMin = .3;//.4;
811  double HeiMax = 6.;
812  double HeiMin = 1.0;
813  double GainRad = 100.;
814  double GainHei = 2000.;
815  double MoveStep = .05;
816  double RadStep = .01;
817  double HeiStep = .1;
818  //Yuliya
819  double *Count = (double *)calloc(NBin*NBin,sizeof(double));
820  // first configuration
821  for(int d=0;d<3;d++) Nano->Pos[d] = OldPos[d];
822  Nano->Rad = MAX(MIN(RadMax,OldPos[3] + 1.0 ),RadMin);
823  Nano->Height = MAX(MIN(HeiMax,OldPos[4] + 0.),HeiMin);
824  double OldNrg = 0.;
825  for(int p=0;p<pNPart();p++){
826  double Dist2 = NanoDist2(Pm[p].Pos,0);
827  if(Dist2 > SQR(2.*Nano->Rad)) continue;
828  if(Pm[p].Typ == 0) OldNrg += kElPhob*Dist2;
829  if(Pm[p].Typ == 1) OldNrg += kElPhil*Dist2;
830  }
831  OldNrg += (GainRad*SQR(Nano->Rad-MinRad) + GainHei*SQR(Nano->Height-MinHei));
832  // fake Monte Carlo
833  for(int i=0;i<NInt;i++){
834  //change position
835  double OldNPos[5] = {Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],Nano->Rad,Nano->Height};
836  for(int m=0;m<3;m++){
837  //change the torus
838  if(m==0){
839  for(int d=0;d<3;d++){
840  Nano->Pos[d] += MoveStep*(2.*Mat->Casuale()-1.);
841  }
842  SetNanoBkf(0);
843  }
844  if(m==1){
845  Nano->Rad += RadStep*(2.*Mat->Casuale()-1.);
846  if(Nano->Rad < RadMin || Nano->Rad > RadMax){
847  Nano->Rad = OldNPos[3];
848  continue;
849  }
850  }
851  if(m==2){
852  Nano->Height += HeiStep*(2.*Mat->Casuale()-1.);
853  if(Nano->Height < HeiMin || Nano->Height > HeiMax){
854  Nano->Height = OldNPos[4];
855  continue;
856  }
857  }
858  // recalc the energy
859  double Nrg = 0.;
860  for(int p=0;p<pNPart();p++){
861  double Dist2 = NanoDist2(Pm[p].Pos,0);
862  if(Dist2 > SQR(2.*Nano->Rad)) continue;
863  if(Pm[p].Typ == 0) Nrg += kElPhob*Dist2;
864  if(Pm[p].Typ == 1) Nrg += kElPhil*Dist2;
865  }
866  Nrg += (GainRad*SQR(Nano->Rad-MinRad) + GainHei*SQR(Nano->Height-MinHei));
867  // accept/remove
868  if(exp(OldNrg-Nrg) > Mat->Casuale()){
869  if(m==0){
870  printf("NewNrg \t%lf Rad ___ Hei ___ Pos %.2f %.2f %.2f\n",Nrg-OldNrg,Nano->Pos[0]-OldNPos[0],Nano->Pos[1]-OldNPos[1],Nano->Pos[2]-OldNPos[2]);
871  }
872  if(m==1){
873  printf("NewNrg \t%lf Rad %.3f Hei ___ Pos ___ ___ ___\n",Nrg-OldNrg,Nano->Rad);
874  }
875  if(m==2){
876  printf("NewNrg \t%lf Rad ___ Hei %.3f Pos ___ ___ ___\n",Nrg-OldNrg,Nano->Height);
877  }
878  OldNrg = Nrg;
879  }
880  else{
881  if(m==0){
882  for(int d=0;d<3;d++){
883  Nano->Pos[d] = OldNPos[d];
884  }
885  SetNanoBkf(0);
886  }
887  if(m==1){
888  Nano->Rad = OldNPos[3];
889  }
890  if(m==2){
891  Nano->Height = OldNPos[4];
892  }
893  }
894  //printf("%lf %lf %lf -> %lf %lf %lf\n",Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],OldNPos[0],OldNPos[1],OldNPos[2]);
895  //printf("%lf-> %lf %lf-> %lf\n",SRad,Nano->Rad,LRad,Nano->Height);
896  }
897  }
898  double Cm[3];
899  double CountCm = 0;
900  for(int p=0;p<pNPart();p++){
901  if(Pm[p].Typ == 1) continue;
902  double Pos[3];
903  for(int d=0;d<3;d++){
904  Pos[d] = Pm[p].Pos[d] - Nano->Pos[d];
905  Pos[d] -= floor(Pos[d]*pInvEdge(d))*pEdge(d);
906  }
907  double Rad = sqrt( SQR(Pos[CLat1]) + SQR(Pos[CLat2]) );
908  if(Rad > Nano->Height) continue;
909  for(int d=0;d<3;d++){
910  Cm[d] += Pm[p].Pos[d];
911  }
912  CountCm += 1.;
913  if(Pm[p].Pos[CNorm] < Nano->Pos[CNorm] - Nano->Rad || Pm[p].Pos[CNorm] > Nano->Pos[CNorm] + Nano->Rad) continue;
914  int vx = (int)(Pos[CLat1]/Nano->Height*NBin);
915  vx += NBin/2;
916  if(vx < 0 || vx >= NBin) continue;
917  int vy = (int)(Pos[CLat2]/Nano->Height*NBin);
918  vy += NBin/2;
919  if(vy < 0 || vy >= NBin) continue;
920  Count[vx*NBin+vy] += 1.;
921  }
922  double Area = 0.;
923  for(int vx=0;vx<NBin;vx++){
924  for(int vy=0;vy<NBin;vy++){
925  if(Count[vx*NBin+vy] < 1.) continue;
926  Area += 1.;
927  }
928  }
929  if(CountCm <= 0.){
930  printf("No particles in the torus\n");
931  return 1;
932  }
933  Nano->Area = SQR(Nano->Height)*Area/(double)(SQR(NBin));
934  for(int d=0;d<3;d++){
935  Cm[d] /= CountCm;
936  //Nano->Pos[d] = Cm[d];
937  CmStalk[d] = Nano->Pos[d];
938  }
939  SetNanoBkf(0);
940  printf("Pos %lf %lf %lf Area %lf Count %lf\n",Nano->Pos[0],Nano->Pos[1],Nano->Pos[2],Nano->Area,CountCm);
941  free(Count);
942  return 0;
943 }
944 #include <Cubo.h>
945 double VarData::NormalWeight(VAR_TRIANGLE *Triang,double *WeightL,int NGrid,int NTri){
946  double Edge[3] = {pEdge(0),pEdge(1),pEdge(2)};
947  NeiVertex VList(NTri,3,NGrid,Edge);
948  double CountStalk = 0.;
949  Vettore Ax(1.,0.,0.);
950  Vettore n1(3);
951  Vettore n2(3);
952  double Max = 0.;
953  for(int t=0;t<NTri;t++){
954  for(int v=0;v<3;v++){
955  int vCurr = Triang[t].v[v];
956  VList.Add(vCurr,t,Triang[t].p[v].x);
957  }
958  }
959  VList.Reorder();
960  VList.SetCounters();
961  for(int t=0;t<NTri;t++){
962  for(int v=0;v<3;v++){
963  int vCurr = Triang[t].v[v];
964  double Weight = 0.;
965  for(VList.SetCounters(vCurr);VList.IfItCell(vCurr);VList.IncrCurr(vCurr)){
966  int tt = VList.VertCurr(vCurr);
967  int p = tt*pNPCh();
968  for(int d=0;d<3;d++){
969  n1.Set(Pm[p].Vel[d],d);
970  }
971  if(n1.Norm() <= 0.) continue;
972  Weight += n2.Angle(&Ax,&n1);
973  }
974  WeightL[vCurr] = Weight;
975  if(Max < Weight) Max = Weight;
976  }
977  }
978  return Max;
979 }
981 void VarData::ConnectLineChain(VAR_LINE *Triang,int NGrid,int NTri){
982  SetNPart(2*NTri);
983  int *Exist = (int *)calloc(NGrid*NGrid,sizeof(int));
984  int NPart = 0;
985  DdLinkedList *Pc = new DdLinkedList(Gen->Edge,pNPart(),1.5);
986  double InvNGrid = 1./(double)NGrid;
987  for(int t=0;t<NTri;t++){
988  for(int v=0;v<2;v++){
989  int vx = (int)(Triang[t].p[v].x[0]*pEdge(0)*InvNGrid);
990  int vy = (int)(Triang[t].p[v].x[1]*pEdge(1)*InvNGrid);
991  int vv = vx*NGrid+vy;
992  Exist[vv] += 1;
993  if(Exist[vv] > 1) continue;
994  Pc->AddPart(NPart,Triang[t].p[v].x);
995  for(int d=0;d<3;d++) Pm[NPart].Pos[d] = Triang[t].p[v].x[d];
996  NPart++;
997  }
998  }
999  SetNPart(NPart);
1000  for(int p=0;p<pNPart();p++){
1001  Ln[p].Link[0] = 0;
1002  }
1003  double DistRel[4];
1004  char FName[60];
1005  int link = 0;
1006  double Pos[5] = {0.,.5*pCm(CLat2),pCm(CLat2),1.5*pCm(CLat2),pEdge(CLat2)};
1007  double pList[4];
1008  for(int p=0,c=0;p<pNPart();p++){
1009  if(Pm[p].Pos[CLat1] <= 0.9 && Pm[p].Pos[CLat2] > Pos[c] + 3.){
1010  Pos[c+1] = Pm[p].Pos[CLat2] + .1;
1011  pList[c] = p;
1012  c++;
1013  if(c==4) break;
1014  }
1015  }
1016  for(int p=pList[0],c=0;c<4;c++,p=pList[c]){
1017  sprintf(FName,"Chain%d.dat",c);
1018  FILE *FChain = fopen(FName,"w");
1019  fprintf(FChain,"%lf %lf\n",Pm[p].Pos[0],Pm[p].Pos[1]);
1020  link = p;
1021  for(int p1=0;p1<pNPart();p1++){
1022  for(Pc->SetCurr(p1);Pc->IfCurr();Pc->NextCurr()){
1023  if(p1 == Pc->p2Curr) continue;
1024  //printf("%d %d %d %d\n",p,p1,Pc->p2Curr,link);
1025  if(link == Pc->p2Curr){
1026  fprintf(FChain,"%lf %lf\n",Pm[p1].Pos[0],Pm[p1].Pos[1]);
1027  Ln[p1].Link[0] = link;
1028  link = p1;
1029  break;
1030  }
1031  }
1032  }
1033  }
1034 }
1036 void VarData::ConnectLineChain3(VAR_LINE *Triang,int NGrid,int NTri){
1037  SetNPart(2*NTri);
1038  int *Called = (int *)calloc(pNPart(),sizeof(int));
1039  int *Exist = (int *)calloc(NGrid*NGrid,sizeof(int));
1040  double *DirPrev = (double *)calloc(3*pNPart(),sizeof(int));
1041  int NChain = 0;
1042  int NPart = 0;
1043  double InvNGrid = 1./(double)NGrid;
1044  DdLinkedList *Pc = new DdLinkedList(Gen->Edge,pNPart(),2.5);
1045  for(int t=0;t<NTri;t++){
1046  for(int v=0;v<2;v++){
1047  int vx = (int)(Triang[t].p[v].x[0]*pEdge(0)*InvNGrid);
1048  int vy = (int)(Triang[t].p[v].x[1]*pEdge(1)*InvNGrid);
1049  int vv = vx*NGrid+vy;
1050  Exist[vv] += 1;
1051  if(Exist[vv] > 1) continue;
1052  Pc->AddPart(NPart,Triang[t].p[v].x);
1053  for(int d=0;d<3;d++) Pm[NPart].Pos[d] = Triang[t].p[v].x[d];
1054  NPart++;
1055  }
1056  }
1057  SetNPart(NPart);
1058  for(int p=0;p<pNPart();p++){
1059  Ln[p].NLink = 0;
1060  Ln[p].Link[0] = -1;
1061  Exist[p] = 0;
1062  Pm[p].CId = p;
1063  }
1064  // for(int p=0;p<pNPart();p++){
1065  // Ln[p].NLink = 1;
1066  // int link = Pc->FindClosest(p);
1067  // if(link == -1) Ln[p].NLink = 0;
1068  // Ln[p].Link[0] = link;
1069  // }
1070  double DistRel[4];
1071  for(int p=0;p<pNPart();p++){
1072  Ln[p].NLink = 1;
1073  double Closest = 1000000.;
1074  int pClosest = -1;
1075  for(Pc->SetCurr(p);Pc->IfCurr();Pc->NextCurr()){
1076  int link = Pc->p2Curr;
1077  if(link == p) continue;
1078  if(Ln[link].Link[0] == p) continue;
1079  Pc->Dist2Curr(DistRel);
1080  double Weight = DirPrev[link*3+0]*DistRel[0] + DirPrev[link*3+1]*DistRel[1];
1081  if(Weight >= 0.) Weight = 0.9;
1082  else if(Weight < 0.) Weight = 1.1;
1083  if(DistRel[3]*Weight < Closest){
1084  Closest = DistRel[3];
1085  pClosest = link;
1086  for(int d=0;d<3;d++){
1087  DirPrev[3*p+d] = DistRel[d];
1088  }
1089  }
1090  }
1091  if(pClosest == -1){
1092  Ln[p].NLink = 0;
1093  continue;
1094  }
1095  Ln[p].Link[0] = pClosest;
1096  Exist[pClosest] += 1;
1097  if(Exist[pClosest] > 1 && Ln[pClosest].NLink == 1) Ln[p].NLink = 0;
1098  }
1099  NChain = 0;
1100  for(int p=0;p<pNPart();p++){
1101  int link = Ln[p].Link[0];
1102  if(Pm[p].CId > NChain){
1103  NChain++;
1104  Pm[p].CId = NChain;
1105  }
1106  Pm[link].CId = Pm[p].CId;
1107  }
1108  NChain = 0;
1109  for(int p=pNPart();p>=0;p--){
1110  int link = Ln[p].Link[0];
1111  if(Pm[p].CId > NChain){
1112  NChain++;
1113  Pm[p].CId = NChain;
1114  }
1115  Pm[link].CId = Pm[p].CId;
1116  }
1117  NChain = 0;
1118  for(int p=0;p<pNPart();p++){
1119  int link = Ln[p].Link[0];
1120  if(Pm[p].CId > NChain){
1121  NChain++;
1122  Pm[p].CId = NChain;
1123  }
1124  Pm[link].CId = Pm[p].CId;
1125  }
1126  SetNChain(NChain);
1127  return;
1128  //isolate multiple connected points
1129  // for(int p=0;p<pNPart();p++){
1130  // int link = Ln[p].Link[0];
1131  // Called[link]++;
1132  // if(Called[link] > 1){
1133  // Ln[p].NLink = 0;
1134  // NPart--;
1135  // }
1136  // }
1137  for(int p=0;p<pNPart()-1;p++){
1138  if(Ln[p].NLink == 0) continue;
1139  int link = Ln[p].Link[0];
1140  int MemPos = -1;
1141  for(int pp=p+1;pp<pNPart();pp++){
1142  if(Pm[pp].Idx == link){
1143  MemPos = pp;
1144  break;
1145  }
1146  }
1147  if(MemPos == -1){
1148  Ln[p].NLink = 0;
1149  continue;
1150  }
1151  printf("%d %d %d %d\n",p,link,MemPos,Pm[p+1].Idx);
1152  //if(p < pNPart() - 2) SwapPart(p+1,p+2);
1153  SwapPart(p+1,MemPos);
1154  printf("%d %d\n",Pm[p+1].Idx,Pm[MemPos].Idx);
1155  }
1156  for(int p=0;p<pNPart()-1;p++){
1157  Ln[p].Link[0] = p+1;
1158  }
1159  // for(int p=0;p<pNPart();p++){
1160  // int link = Ln[p].Link[0];
1161  // for(int pp=0;pp<pNPart();pp++){
1162  // if(Pm[pp].Idx == link){
1163  // Ln[p].Link[0] = pp;
1164  // break;
1165  // }
1166  // }
1167  // }
1168  SetNPart(NPart);
1169  NChain = 0;
1170  for(int p=0;p<pNPart();p++){
1171  Pm[p].CId = NChain;
1172  if(Ln[p].NLink == 0) NChain++;
1173  }
1174  SetNChain(NChain);
1175  free(Exist);
1176  free(Called);
1177  free(DirPrev);
1178 }
1180 void VarData::ConnectLineChain2(VAR_LINE *Triang,int NGrid,int NTri){
1181  //can't follow degenerate triangles..
1182  double Edge[3] = {pEdge(0),pEdge(1),pEdge(2)};
1183  NeiVertex VList(NTri,2,NGrid,Edge);
1184  double CountStalk = 0.;
1185  Vettore Ax(1.,0.,0.);
1186  Vettore n1(3);
1187  Vettore n2(3);
1188  double Max = 0.;
1189  int NTria = 0;
1190  for(int t=0;t<NTri;t++){
1191  if(Triang[t].p[0].x[0] == Triang[t].p[1].x[0] && Triang[t].p[0].x[1] == Triang[t].p[1].x[1]) continue;
1192  for(int v=0;v<2;v++){
1193  int vCurr = Triang[t].v[v];
1194  VList.Add(vCurr,NTria,Triang[t].p[v].x);
1195  for(int d=0;d<3;d++) Pm[NTria*2+v].Pos[d] = Triang[t].p[v].x[d];
1196  NTria++;
1197  }
1198  }
1199  VList.Reorder();
1200  SetNPart(2*NTri);
1201  //VList.Print();
1202  VList.SetCounters();
1203  for(int p=0;p<pNPart();p++){
1204  Ln[p].NLink = 0;
1205  }
1206  for(int t=0;t<NTri;t++){
1207  for(int v=0;v<2;v++){
1208  int vCurr = Triang[t].v[v];
1209  for(VList.SetCounters(vCurr);VList.IfItCell(vCurr);VList.IncrCurr(vCurr)){
1210  int tt = VList.TriaCurr(vCurr);
1211  int p = tt*2+v;
1212  if(vCurr > tt){
1213  Ln[p].NLink = 1;
1214  Ln[p].Link[0] = vCurr;
1215  }
1216  }
1217  }
1218  }
1219  int nChain=0;
1220  for(int c=0;c<pNChain();c++){
1221  Ch[c].NPCh = 0;
1222  }
1223  for(int p=0;p<pNPart();p++){
1224  if(Pm[p].CId <= nChain) continue;
1225  Ch[nChain].InitBead = p;
1226  int pp1 = Ln[p].Link[0];
1227  int pp = 0;
1228  printf("------ %d %d\n",p,pp1);
1229  for(pp=pp1;Ln[pp].NLink > 0;pp=Ln[pp].Link[0]){
1230  printf(" %d %d\n",pp,Ch[nChain].NPCh);
1231  Pm[pp1].CId = nChain;
1232  Ch[nChain].NPCh++;
1233  pp = Ln[pp].Link[0];
1234  if(pp == pp1) break;
1235  }
1236  Ch[nChain].EndBead = pp;
1237  nChain++;
1238  }
1239  SetNChain(nChain);
1240  for(int c=0;c<pNChain();c++){
1241  printf("%d %d %d %d\n",c,Ch[c].InitBead,Ch[c].EndBead,Ch[c].NPCh);
1242  for(int p=Ch[c].InitBead,ppc=0;ppc<Ch[c].NPCh;p=Ln[p].Link[0]){
1243  printf("%d %d \n",c,p);
1244  }
1245  }
1246 }
1247 int VarData::StalkPos(double *OldPos){
1248  double CmStalk[3] = {OldPos[0],OldPos[1],OldPos[2]};
1249  //pPos(OldPos);
1250  //if(StalkPos4(OldPos,CmStalk)) return 1;
1251  int NTrials = 10;
1252  int n=0;
1253  while(StalkPos4(OldPos,CmStalk)){
1254  n++;
1255  if(n > NTrials) continue;
1256  }
1257  //pPos(CmStalk);
1258  for(int d=0;d<3;d++){
1259  //Nano->Pos[d] = CmStalk[d];// + OldPos[d];
1260  //Nano->PosBf[d] = Nano->Pos[d];// - floor(Nano->Pos[d]*pInvEdge(d))*pEdge(d);
1261  OldPos[d] = pNanoPos(0,d);;
1262  }
1263  OldPos[3] = Nano->Rad;
1264  OldPos[4] = Nano->Height;
1265  Nano->Axis[CLat1] = 0.;
1266  Nano->Axis[CLat2] = 0.;
1267  Nano->Axis[CNorm] = 1.;
1268  // Nano->Rad = .5;
1269  // Nano->Height = 4.;
1270  // Nano->Hamaker = 1.;
1271  // Nano->OffSet = 1.;
1272  // Nano->Shape = SHAPE_STALK;
1273  return 0;
1274 }
1276  const int NGrid = 26;
1277  double *Plot = (double *)calloc(NGrid*NGrid,sizeof(double));
1278  double *Plot2 = (double *)calloc(NGrid*NGrid,sizeof(double));
1279  double Cm[3] = {0.,0.,0.};
1280  double Cm1[3] = {0.,0.,0.};
1281  double Incr = 500.;
1282  double Threshold = SQR(SQR(1./Incr));
1283  int NReweight = 2;
1284  for(int gx=0;gx<NGrid;gx++){
1285  for(int gy=0;gy<NGrid;gy++){
1286  Plot[gx*NGrid+gy] = 1.;
1287  }
1288  }
1289  for(int p=0;p<pNPart();p++){
1290  if(Pm[p].Typ != 0 ) continue;
1291  int gx = (int)(Pm[p].Pos[CLat1]*pInvEdge(CLat1)*NGrid);
1292  int gy = (int)(Pm[p].Pos[CLat2]*pInvEdge(CLat2)*NGrid);
1293  if(gx < 0 || gx >= NGrid) continue;
1294  if(gy < 0 || gy >= NGrid) continue;
1295  Plot[gx*NGrid+gy] += Incr;
1296  }
1297  double Dx = .5*pEdge(CLat1)/(double)NGrid;
1298  double Dy = .5*pEdge(CLat2)/(double)NGrid;
1299  double Count = 0.;
1300  //FILE *Ciccia = fopen("PoreGrid.dat","w");
1301  double Area = 0.;
1302  //weigthing by neighbours
1303  for(int gx=0;gx<NGrid;gx++){
1304  for(int gy=0;gy<NGrid;gy++){
1305  Plot2[gx*NGrid+gy] = Plot[gx*NGrid+gy];
1306  }
1307  }
1308  for(int gx=0;gx<NGrid;gx++){
1309  for(int gy=0;gy<NGrid;gy++){
1310  if(gx == 0){ Plot[gx*NGrid+gy] = 1./SQR(SQR(Plot2[gx*NGrid+gy]));continue;}
1311  if(gy == 0){ Plot[gx*NGrid+gy] = 1./SQR(SQR(Plot2[gx*NGrid+gy]));continue;}
1312  if(gx == NGrid-1){ Plot[gx*NGrid+gy] = 1./SQR(SQR(Plot2[gx*NGrid+gy]));continue;}
1313  if(gy == NGrid-1){ Plot[gx*NGrid+gy] = 1./SQR(SQR(Plot2[gx*NGrid+gy]));continue;}
1314  Plot[gx*NGrid+gy] = 1./(Plot2[(gx-1)*NGrid+gy]*Plot2[(gx+1)*NGrid+gy]*Plot2[gx*NGrid+(gy-1)]*Plot2[gx*NGrid+(gy+1)]);
1315  }
1316  }
1317  //above the threshold
1318  for(int gx=0;gx<NGrid;gx++){
1319  double x = gx/(double)NGrid*pEdge(CLat1) + Dx;
1320  for(int gy=0;gy<NGrid;gy++){
1321  double y = gy/(double)NGrid*pEdge(CLat2) + Dy;
1322  //Plot[gx*NGrid+gy] /= 100.;
1323  //fprintf(Ciccia,"%lf %lf %lf 0\n",x,y,10.*Plot[gx*NGrid+gy]);
1324  Cm[CLat1] += x*Plot[gx*NGrid+gy];
1325  Cm[CLat2] += y*Plot[gx*NGrid+gy];
1326  Count += Plot[gx*NGrid+gy];
1327  if(Plot[gx*NGrid+gy] > Threshold ){
1328  Area += 1.;//Plot[gx*NGrid+gy];
1329  //fprintf(Ciccia,"%lf %lf %lf %d\n",x,y,1.,2);
1330  }
1331  //else fprintf(Ciccia,"%lf %lf %lf %d\n",x,y,.5,0);
1332  }
1333  }
1334  Area = Area/(double)SQR(NGrid)*pEdge(CLat1)*pEdge(CLat2);
1335  Area /= 1.;//Count;
1336  Nano->Rad = sqrt(Area/M_PI);
1337  Cm[0] /= Count;
1338  Cm[1] /= Count;
1339  //reweighting wrt to the former position
1340  for(int r=0;r<NReweight;r++){
1341  Count = 0.;
1342  Cm1[0] = Cm[0];Cm1[1] = Cm[1];
1343  Cm[0] = 0.;Cm[1] = 0.;
1344  Area = 0.;
1345  //printf("%lf %lf\n",Cm1[0],Cm1[1]);
1346  for(int gx=0;gx<NGrid;gx++){
1347  double x = gx/(double)NGrid*pEdge(CLat1) + Dx;
1348  for(int gy=0;gy<NGrid;gy++){
1349  double y = gy/(double)NGrid*pEdge(CLat2) + Dy;
1350  double Dist = (SQR(x-Cm1[CLat1]) + SQR(y-Cm1[CLat2]));
1351  if(Dist > SQR(1.1*Nano->Rad)) continue;
1352  double Weight = 1./pow(Dist,.2);
1353  //fprintf(Ciccia,"%lf %lf %lf 1\n",x,y,10.*Plot[gx*NGrid+gy]*Weight);
1354  Cm[CLat1] += x*Plot[gx*NGrid+gy]*Weight;
1355  Cm[CLat2] += y*Plot[gx*NGrid+gy]*Weight;
1356  Count += Plot[gx*NGrid+gy]*Weight;
1357  if(Plot[gx*NGrid+gy] > Threshold ) Area += 1.;
1358  }
1359  }
1360  Cm[0] /= Count;
1361  Cm[1] /= Count;
1362  Area = Area/(double)SQR(NGrid)*pEdge(CLat1)*pEdge(CLat2);
1363  //Nano->Rad = sqrt(Area/M_PI);
1364  }
1365  //assigning
1366  Cm[2] = pCm(2);
1367  for(int d=0;d<3;d++){
1368  Nano->Pos[d] = Cm[d];
1369  }
1370  SetNanoBkf(0);
1371  //asymmetry
1372  double AreaOut = 0.;
1373  double AreaIn = 0.;
1374  for(int gx=0;gx<NGrid;gx++){
1375  double x = gx/(double)NGrid*pEdge(CLat1) + Dx;
1376  for(int gy=0;gy<NGrid;gy++){
1377  double y = gy/(double)NGrid*pEdge(CLat2) + Dy;
1378  double Dist = SQR(x-Cm[0]) + SQR(y-Cm[1]);
1379  if(Dist <= SQR(Nano->Rad)){
1380  AreaIn += 1.;
1381  if(Plot[gx*NGrid+gy] < Threshold ) AreaOut += 1.;
1382  //fprintf(Ciccia,"%lf %lf %lf %d\n",x,y,0.,1);
1383  }
1384  else{
1385  if(Plot[gx*NGrid+gy] > Threshold ) AreaOut += 1.;
1386  }
1387  }
1388  }
1389  //fclose(Ciccia);
1390  free(Plot);
1391  return AreaOut/AreaIn;
1392 }
void AddPart(const int p, double *Pos)
Add a particle to the cell c.
Definition: Cubo.cpp:276
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
int CId
Chain Identifier.
Definition: VarData.h:224
int VertCurr(int v)
Current vertex iterator for the vertex v.
Definition: Cubo.cpp:1109
int InterBSpline1D(double *PlIn, double *PmOut, int NIn, int NOut)
1-d BSpline
double x[3]
Cartesian coordinates.
Definition: VarData.h:486
double NanoDist2(double *Pos, int n)
Pointer to a generic function.
Definition: VarData.h:764
int Type
Type of the chain (see list CHAIN_)
Definition: VarData.h:246
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
double Vel[4]
xyzr Velocity of the particle
Definition: VarData.h:220
LINKS * Ln
Array of linking between the particles.
Definition: VarData.h:1048
double Norm()
Return the norm of a Vettore.
double ShiftPos[3]
Reference position.
Definition: VarData.h:1058
NANO * Nano
Extra particle.
Definition: VarData.h:1044
Define a triangle.
Definition: VarData.h:509
double Bkf[3]
xyz Backfold distance
Definition: VarData.h:218
double Normalize()
Normlizes a Vettore.
int SetNChain(int NewNCh)
Set and reallocate the number of chains.
Definition: VarDataComm.cpp:78
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
double Edge[4]
xyzr edges of the simulation box
Definition: VarData.h:309
double Vel[3]
Chain velocity.
Definition: VarData.h:242
int IfItCell(int v)
End of the counters.
Definition: Cubo.cpp:1102
int BfDefChain()
Definition of the chain.
char Name[60]
BLock name.
Definition: VarData.h:257
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
double pCm(int d)
Center of mass of the system.
Definition: VarData.h:924
Connects the triangles by vertices.
Definition: Cubo.h:480
int CLat2
lateral coordinate
Definition: VarData.h:1078
Information for every block.
Definition: VarData.h:255
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.
XYZ p[3]
The three vertices.
Definition: VarData.h:500
double Val(int N)
Value of the N column.
double * x
Where the data are stored.
void BfPep()
Backfold the nano described as a cluster of monomers.
void Dist2Curr(double *DistRel)
Retrun the squared current interparticle distance.
Definition: Cubo.cpp:502
void ConvoluteMatrix(double *Plot, int NGrid, int NDim, int IfMinImConv)
Convolute with a matrix.
double m
y = m*x + q
int InitBead
Initial bead.
Definition: VarData.h:250
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
int v[3]
Reference for the vectors.
Definition: VarData.h:504
void SwapPart(int p1, int p2)
Swap two particle.
Definition: VarDataEl.cpp:378
int InitIdx
Initial particle position.
Definition: VarData.h:265
double Pos[3]
Position.
Definition: VarData.h:427
int Shape
0 none, 1 spherical, 2 cylindrical 3 wall
Definition: VarData.h:473
Linear interpolation.
double pNanoPos(int n, int d)
Return back folded nano position.
int SetNPart(int NewNPart)
Set and reallocate the number of particles.
Definition: VarDataComm.cpp:2
int StalkPos4(double *OldPos, double *CmStalk)
Find the position of the stalk forth method.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
int pNPCh()
Number of particle per chain.
void Add(int v, int t, double *Pos)
Add the triangle t at the vertex v.
Definition: Cubo.cpp:1036
double PorePos()
Find the position of the pore.
int pNBlock()
Number of blocks.
int StalkPos(double *OldPos)
Find the position of the stalk.
double Dir[3]
Directive of the chain.
Definition: VarData.h:240
void BackBone(double *Line, int NBin)
Describe the backbone of a filament.
int NChain
Number of chain.
Definition: VarData.h:347
void IncrCurr(int v)
Increment the counter.
Definition: Cubo.cpp:1106
double Casuale()
Random uniform number.
double Val(int row)
Returns a value in 1d.
void Reorder()
Reorder and fill the vertices.
Definition: Cubo.cpp:1058
void NextCurr()
Increase the iterator to the next couple.
Definition: Cubo.cpp:479
void Point2Shape(int iShape)
Point to the shape function.
Definition: VarDataEl.cpp:490
int NPCh
Number of particle per chain.
Definition: VarData.h:349
bool ShiftSys(int How)
Shift the system accordin to the SHIFT_ definitions.
void StalkPos2(double *OldPos, double *CmStalk)
Find the position of the stalk second method.
void ShiftRef(int BackFold)
Backfold the system wrt the reference position.
int SysType
Contains the definition of the system.
Definition: VarData.h:1086
void ConnectLineChain2(VAR_LINE *Triang, int NGrid, int NTri)
Connect the lines in a chain.
int pNNano()
Number of nanoparticles.
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
double Cm[3]
Center of mass of the system.
Definition: VarData.h:313
void SetCounters(int v)
Set counters to zero for the point v.
Definition: Cubo.cpp:1093
void WriteSurf(FILE *F2Write, double **Plot, int NSample, int OffSet)
Write the particle position as linked edges of squares.
double Pos[4]
xyzr Postion of the chain
Definition: VarData.h:238
int NPCh
particles per chain
Definition: VarData.h:259
void DefBlock(int *NChStep, int How)
Define four different blocks.
Matrice computes the algebric operations on matrices.
int TriaCurr(int v)
Current triangle for the vertex v.
Definition: Cubo.cpp:1112
void ConnectLineChain3(VAR_LINE *Triang, int NGrid, int NTri)
Connect the lines in a chain.
int EndBead
End bead.
Definition: VarData.h:252
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int NChain
chains
Definition: VarData.h:263
void StalkPos3(double *OldPos, double *CmStalk)
Find the position of the stalk third method.
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
int pNChain()
Number of chain.
int BfEdge()
Find the box size if missing.
double Angle(Vettore *u, Vettore *v)
Computes the angle between two Vetttore.
int CLat1
lateral coordinate
Definition: VarData.h:1076
int IfCurr()
Tell when the curr loop is over.
Definition: Cubo.cpp:498
bool BackFold(int How)
Backfold the particle position.
int p2Curr
Current particle.
Definition: Cubo.h:51
double Axis[3]
Rotation axis.
Definition: VarData.h:435
int NPart
Number of particle.
Definition: VarData.h:345
Domain decomposition as pointer to linked particles.
Definition: Cubo.h:162
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int NPCh
Number of particles per chain.
Definition: VarData.h:248
int NBlock
Number of blocks.
Definition: VarData.h:359
void VarMessage(const char *s,...)
If enabled call the function position.
Definition: VarData.cpp:29
void SetCurr(int p)
Gather information of the neighbouring cells.
Definition: Cubo.cpp:453
void Set(double Val, int Col)
Set the N column.
double q
y = m*x + q
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.
RETTA InterRett(double *Px, double *Py, int NMass)
Linear interpolation.