Allink  v0.1
VarDataEl.cpp
1 /***********************************************************************
2 VarDataEl: Elaboration functions for the VarData class. This functions
3 provides a simple manipulation of the data read by [Open]. The
4 options are provided to elaborate different system's shapes.
5 Copyright (C) 2008 by Giovanni Marelli <sabeiro@virgilio.it>
6 
7 
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or
11 (at your option) any later version.
12 
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 ***********************************************************************/
22 #include "../include/VarData.h"
23 
24 int VarData::RadDistr(int Values,double *Plot,double Border[2],int How){
25  return 1;
26 }
28  char *Info;
29  Info = (char *)malloc(256*sizeof(char));
30  Gen->Temp = 0.;
31  double VelSquare[3];
32  for(int d=0;d<3;d++){
33  Gen->Pre[d] = 0.;
34  VelSquare[d] = 0.;
35  }
36  for(int p=0;p<Gen->NPart;p++){
37  for(int d=0;d<3;d++){
38  Gen->Pre[d] += SQUARE(Pm[p].Vel[d]);
39  }
40  }
41  Gen->Temp = (Gen->Pre[0] + Gen->Pre[1] + Gen->Pre[2])/(3.*Gen->NPart);
42  for(int d=0;d<3;d++){
43  Gen->Pre[d] /= (Gen->Edge[0]*Gen->Edge[1]*Gen->Edge[2]);
44  }
45  Gen->SurfTens = Gen->Edge[2]*(Gen->Pre[2] - .5*(Gen->Pre[0] + Gen->Pre[1]) );
46  sprintf(Info,"Pre0\tPre1\tPre2\tSurfTens\tTemp\n%.4g\t%.4g\t%.4g\t%.4g\t%.4g\n",Gen->Pre[0],Gen->Pre[1],Gen->Pre[2],Gen->SurfTens,Gen->Temp);
47  return Info;
48 }
50  if(IfPlotMem){
51  MOMENTI m1;
52  m1.Num = 0;
53  }
54  PlotMem = (double *) calloc(SQR(NSample),sizeof(double));
55  double *PlotR = (double *) calloc(SQR(NEdge),sizeof(double));
56  //LoadDensFile(PlotR,NEdge);
57  MOMENTI m1 = SampleSurfacePart(PlotR,NEdge,0);
58  InterBSpline2D(PlotR,PlotMem,NEdge,NSample);
59  if(IfPlotMem){
60  free(PlotR);
61  }
62  IfPlotMem = 1;
63  return m1;
64 }
65 void VarData::SampleSurface(double *Plot,int NSample,int Type){
66  double Average=0.;
67  double NAverage=0;
68  double **Norma = (double **)calloc(NSample,sizeof(double));
69  for(int v=0;v<NSample;v++){
70  Norma[v] = (double *)calloc(NSample,sizeof(double));
71  for(int vv=0;vv<NSample;vv++){
72  Plot[v*NSample+vv] = 0.;
73  }
74  }
75  if(Type == 0){
76  for(int c=0;c<Gen->NChain;c++){
77  //if(!CHAIN_IF_TYPE(Ch[c].Type,NChType))continue;
78  Average += Ch[c].Pos[CNorm];NAverage += 1.;
79  int v = (int)(Ch[c].Pos[CLat1]*pInvEdge(CLat1)*NSample);
80  if( v < 0 || v >= NSample) continue;
81  int vv = (int)(Ch[c].Pos[CLat2]*pInvEdge(CLat2)*NSample);
82  if( vv < 0 || vv >= NSample) continue;
83  Plot[v*NSample+vv] += Ch[c].Pos[CNorm];
84  Norma[v][vv] += 1.;
85  //printf("%d %d %lf %lf\n",v,vv,Plot[v][vv],Norma[v][vv]);
86  }
87  }
88  else {
89  for(int p=0;p<Gen->NPart;p++){
90  Average += Pm[p].Pos[CNorm];NAverage += 1.;
91  int v = (int)(Pm[p].Pos[CLat1]*pInvEdge(CLat1)*NSample);
92  if( v < 0 || v >= NSample) continue;
93  int vv = (int)(Pm[p].Pos[CLat2]*pInvEdge(CLat2)*NSample);
94  if( vv < 0 || vv >= NSample) continue;
95  Plot[v*NSample+vv] += Pm[p].Pos[CNorm];
96  Norma[v][vv] += 1.;
97  //printf("%d %d %lf %lf\n",v,vv,Plot[v][vv],Norma[v][vv]);
98  }
99  }
100  Average /= (double)NAverage;
101  for(int v=0;v<NSample;v++){
102  for(int vv=0;vv<NSample;vv++){
103  if(Norma[v][vv] > 0.){
104  Plot[v*NSample+vv] /= Norma[v][vv];
105  }
106  else Plot[v*NSample+vv] = Average;
107  }
108  }
109  for(int v=0;v<NSample;v++)
110  free(Norma[v]);
111  free(Norma);
112 }
113 MOMENTI VarData::SampleSurface(Matrice *Plot,int NSample,int Type){
114  MOMENTI m1;
115  double Average=0.;
116  double NAverage=0;
117  double **Norma = (double **)calloc(NSample,sizeof(double));
118  for(int v=0;v<NSample;v++){
119  Norma[v] = (double *)calloc(NSample,sizeof(double));
120  for(int vv=0;vv<NSample;vv++){
121  Plot->Set(v,vv,0);
122  }
123  }
124  for(int c=0;c<Gen->NChain;c++){
125  if(!CHAIN_IF_TYPE(Ch[c].Type,NChType))continue;
126  Average += Ch[c].Pos[CNorm];NAverage += 1.;
127  int v = (int)(Ch[c].Pos[CLat1]*pInvEdge(CLat1)*NSample);
128  if( v < 0 || v >= NSample) continue;
129  int vv = (int)(Ch[c].Pos[CLat2]*pInvEdge(CLat2)*NSample);
130  if( vv < 0 || vv >= NSample) continue;
131  Plot->Add(v,vv,Ch[c].Pos[CNorm]);
132  Norma[v][vv] += 1.;
133  //printf("%d %d %'lf\n",v,vv,Plot[v][vv]);
134  }
135  Average /= (double)NAverage;
136  for(int v=0;v<NSample;v++){
137  for(int vv=0;vv<NSample;vv++){
138  if(Norma[v][vv] > 0.)
139  Plot->Set(v,vv,Plot->Val(v,vv)/Norma[v][vv]);
140  //pasticcio
141  else Plot->Set(v,vv,Average);
142  //printf("%d %d %'lf\n",v,vv,Plot[v][vv]);
143  }
144  }
145  free(Norma);
146  m1.Uno = Average;
147  m1.Num = SQR(NSample);
148  return m1;
149 }
150 void VarData::LoadDensFile(double **Plot,int NSample){
151  double Round = 0.001;
152  double *Count = (double *)calloc(3*NSample*NSample,sizeof(double));
153  for(int p=0;p<pNPart();p++){
154  int sr = (int)((pPos(p,0)+Round)*pInvEdge(0)*NSample);
155  if(sr < 0 || sr >= NSample) continue;
156  int sz = (int)((pPos(p,1)+Round)*pInvEdge(1)*NSample);
157  if(sz < 0 || sz >= NSample) continue;
158  int t = pType(p);
159  Plot[t][sr*NSample+sz] += pPos(p,2);
160  Count[(sr*NSample+sz)*3+t] += 1.;
161  }
162  Matrice Mask(5,5);
163  Mask.FillGaussian(.5,3.);
164  for(int t=0;t<3;t++){
165  for(int s=0;s<NSample*NSample;s++){
166  if(Count[s*3+t] > 0.)
167  Plot[t][s] /= Count[s*3+t];
168  }
169  }
170  int NDim = 2;
171  int IfMinImConv = 1;
172  for(int t=0;t<3;t++){
173  Mask.ConvoluteMatrix(Plot[t],NSample,NDim,IfMinImConv);
174  Mask.ConvoluteMatrix(Plot[t],NSample,NDim,IfMinImConv);
175  }
176  free(Count);
177 }
178 MOMENTI VarData::SampleSurfacePart(double *Plot,int NSample,int Type){
179  double Round = 0.001;
180  double Average=0.;
181  double NAverage=0;
182  double Min = 10000000.;
183  double Max =-10000000.;
184  double **Norma = (double **)calloc(NSample,sizeof(double));
185  MOMENTI m1;
186  for(int v=0;v<NSample;v++){
187  Norma[v] = (double *)calloc(NSample,sizeof(double));
188  for(int vv=0;vv<NSample;vv++){
189  Plot[v*NSample+vv] = 0;
190  }
191  }
192  for(int p=0;p<Gen->NPart;p++){
193  if(Pm[p].Typ != Type) continue;
194  int v = (int)((Pm[p].Pos[CLat1]+Round)*pInvEdge(CLat1)*NSample);
195  if( v < 0 || v >= NSample) continue;
196  int vv = (int)((Pm[p].Pos[CLat2]+Round)*pInvEdge(CLat2)*NSample);
197  if( vv < 0 || vv >= NSample) continue;
198  Plot[v*NSample+vv] += Pm[p].Pos[CNorm];
199  Norma[v][vv] += 1.;
200  Average += Pm[p].Pos[CNorm];NAverage += 1.;
201  }
202  Average /= (double)NAverage;
203  for(int v=0;v<NSample;v++){
204  for(int vv=0;vv<NSample;vv++){
205  if(Norma[v][vv] > 0.){
206  Plot[v*NSample+vv] /= Norma[v][vv];
207  }
208  else Plot[v*NSample+vv] = Average;
209  if(Plot[v*NSample+vv] < Min) Min = Plot[v*NSample+vv];
210  if(Plot[v*NSample+vv] > Max) Max = Plot[v*NSample+vv];
211  }
212  }
213  for(int v=0;v<NSample;v++)
214  free(Norma[v]);
215  free(Norma);
216  m1.Min = Min;
217  m1.Max = Max;
218  m1.Uno = Average;
219  m1.Num = SQR(NSample);
220  return m1;
221 }
222 int VarData::SpatialDerivative(Matrice *Surface,Matrice *Resp,SPLINE Weight,int NSample){
223  SampleSurface(Surface,NSample,NChType);
224  Matrice *Mask = new Matrice(Weight);
225  for(int h=0;h<Resp->Size();h++)
226  for(int w=0;w<Resp->Size();w++)
227  Resp->Set(h,w,0.);
228  Mat->ApplyFilter(Surface,Resp,Mask);
229  Mask->Transpose();
230  Mat->ApplyFilter(Surface,Resp,Mask);
231  delete Mask;
232  return 0;
233 }
235  Properties *Pr = new Properties();
236  double StructPhil = 0.;
237  double StructPhob = 0.;
238  for(int c=0;c<Gen->NChain;c++){
239  for(int cc=0;cc<Gen->NChain;cc++){
240  Pr->ChDiff += QUAD((Ch[c].Pos[CLat1]-Ch[cc].Pos[CLat1]));
241  Pr->ChDiff += QUAD((Ch[c].Pos[CLat2]-Ch[cc].Pos[CLat2]));
242  }
243  double CmPhil[3] = {0.,0.,0.};
244  double CmPhob[3] = {0.,0.,0.};
245  double RadPhil=0.;
246  double RadPhob=0.;
247  double NPhil=0.;
248  double NPhob=0.;
249  double q=1./Gen->Edge[3];
250  double Dist[3] = {0.,0.,0.};
251  for(int p=c*pNPCh();p<pNPCh()*c+Block[0].Asym;p++){//Phob
252  for(int d=0;d<3;d++){
253  CmPhob[d] += Pm[p].Pos[d];
254  for(int pp=c*Gen->NPCh;pp<Gen->NPCh*c+Block[0].Asym-1;pp++){//Phob
255  Dist[d] = Pm[p].Pos[d] - Pm[pp].Pos[d];
256  if(p == pp) StructPhob += 1.;
257  else
258  //StructPhob += Dist[d]*sin(q*Dist[d])/(q);
259  StructPhob += exp(-QUAD(( q*Dist[d] ))/6. );
260  }
261  if(p == Gen->NPCh*c+Block[0].Asym-1) continue;
262  Pr->RePhob += QUAD(( Pm[p].Pos[d] - Pm[p+1].Pos[d]));
263  }
264  }
265  for(int p=c*Gen->NPCh+Block[0].Asym;p<Gen->NPCh*(c+1)-1;p++){//Phil
266  for(int d=0;d<3;d++){
267  CmPhil[d] += Pm[p].Pos[d];
268  for(int pp=c*Gen->NPCh+Block[0].Asym;pp<Gen->NPCh*(c+1)-1;pp++){//Phil
269  Dist[d] = Pm[p].Pos[d] - Pm[pp].Pos[d];
270  if(p == pp) StructPhil += 1.;
271  else
272  // StructPhil += sin(q*Dist[d])/(q*Dist[d]);
273  StructPhil += exp(-QUAD(( q*Dist[d] ))/6. );
274  }
275  if( p == Gen->NPCh*(c+1)-1) continue;
276  Pr->RePhil += QUAD(( Pm[p].Pos[d] - Pm[p+1].Pos[d]));
277  }
278  }
279  for(int d=0;d<3;d++){
280  CmPhob[d] /= (double) Block[0].Asym;
281  CmPhil[d] /= (double) Gen->NPCh-Block[0].Asym;
282  }
283  for(int p=c*Gen->NPCh;p<Gen->NPCh*c+Block[0].Asym;p++){//Phob
284  for(int d=0;d<3;d++)
285  RadPhob += QUAD(( CmPhob[d] - Pm[p].Pos[d]));
286  }
287  for(int p=c*Gen->NPCh+Block[0].Asym;p<Gen->NPCh*(c+1);p++){//Phil
288  for(int d=0;d<3;d++)
289  RadPhil += QUAD(( CmPhil[d] - Pm[p].Pos[d]));
290  }
291  Pr->GyrPhob += RadPhob / (double)(Block[0].Asym*Gen->NChain);
292  Pr->GyrPhil += RadPhil / (double)((Gen->NPCh - Block[0].Asym)*Gen->NChain);
293  }
294  Pr->RePhob /= Gen->NChain;
295  Pr->RePhil /= Gen->NChain;
296  Pr->ChDiff /= (double)QUAD((Gen->NChain));
297  Pr->FactPhob += StructPhob/(double)(Gen->NChain);//*QUAD((Asym)));
298  Pr->FactPhil += StructPhil/(double)(Gen->NChain);//*QUAD((Gen->NPCh-Asym)));
299  Pr->Print();
300  return *Pr;
301 }
303  double *Interp1 = (double *)calloc(Gen->NPCh,sizeof(double));
304  double *Interp2 = (double *)calloc(Gen->NPCh,sizeof(double));
305  double *InterpN = (double *)calloc(Gen->NPCh,sizeof(double));
306  RETTA r1;
307  RETTA r2;
308  for(int c=0;c<Gen->NChain;c++){
309  double Comp[3]={0.,0.,0.};
310  for(int p= c*Gen->NPCh,pp=0;p<(c+1)*Gen->NPCh;p++,pp++){
311  Interp1[pp] = Pm[p].Pos[CLat1];
312  Interp2[pp] = Pm[p].Pos[CLat2];
313  InterpN[pp] = Pm[p].Pos[CNorm];
314  if(p == (c+1)*Gen->NPCh - 1) break;
315  double Segm[3]={0.,0.,0.};
316  double Rad=0.;
317  for(int d=0;d<3;d++){
318  Segm[d] = Pm[p].Pos[d] - Pm[p+1].Pos[d];
319  Rad += QUAD((Segm[d]));
320  }
321  Rad = sqrt(Rad);
322  for(int d=0;d<3;d++){
323  Comp[d] += Segm[d] / (Rad*Gen->NPCh);
324  }
325  }
326  if( ASS((Comp[CNorm])) < .7)
327  Ch[c].Type |= CHAIN_FLABBY;
328  else
329  Ch[c].Type |= CHAIN_STRETCH;
330  r1 = Mat->InterRett(Interp1,InterpN,Gen->NPCh);
331  r2 = Mat->InterRett(Interp2,InterpN,Gen->NPCh);
332  if( POS(r1.m) < .5 || POS(r2.m) < .5)
333  Ch[c].Type |= CHAIN_TILTED;
334  }
335  free(Interp1);
336  free(Interp2);
337  free(InterpN);
338  return 0;
339 }
340 void VarData::ChangeNChain(int NChain,int nBlock){
341  int NChainTemp = 0;
342  Block[nBlock].NChain = NChain;
343  Block[nBlock].EndIdx = Block[nBlock].InitIdx + Block[nBlock].NChain*Block[nBlock].NPCh;
344  NChainTemp = Block[0].NChain;
345  for(int b=1;b<Gen->NBlock;b++){
346  Block[b].InitIdx = Block[b-1].EndIdx;
347  Block[b].EndIdx = Block[b].InitIdx + Block[b].NChain*Block[b].NPCh;
348  NChainTemp += Block[b].NChain;
349  }
350  Gen->NChain = NChainTemp;
351 }
352 void VarData::SwapChain(int c1,int c2){
353  SwapChain(c1,c2,0);
354 }
355 void VarData::SwapChain(int c1,int c2,int b){
356  int p1 = Block[b].InitIdx+c1*Block[b].NPCh;
357  int p2 = Block[b].InitIdx+c2*Block[b].NPCh;
358  double Temp[3];
359  for(int p=0;p<Block[b].NPCh;p++){
360  for(int d=0;d<3;d++){
361  Temp[d] = Pm[p+p1].Pos[d];
362  Pm[p+p1].Pos[d] = Pm[p+p2].Pos[d];
363  Pm[p+p2].Pos[d] = Temp[d];
364  Temp[d] = Pm[p+p1].Vel[d];
365  Pm[p+p1].Vel[d] = Pm[p+p2].Vel[d];
366  Pm[p+p2].Vel[d] = Temp[d];
367  }
368  }
369  for(int d=0;d<3;d++){
370  Temp[d] = Ch[c1].Pos[d];
371  Ch[c1].Pos[d] = Ch[c2].Pos[d];
372  Ch[c2].Pos[d] = Temp[d];
373  Temp[d] = Ch[c1].Vel[d];
374  Ch[c1].Vel[d] = Ch[c2].Vel[d];
375  Ch[c2].Vel[d] = Temp[d];
376  }
377 }
378 void VarData::SwapPart(int p1,int p2){
379  double Temp[3];
380  for(int d=0;d<3;d++){
381  Temp[d] = Pm[p1].Pos[d];
382  Pm[p1].Pos[d] = Pm[p2].Pos[d];
383  Pm[p2].Pos[d] = Temp[d];
384  Temp[d] = Pm[p1].Vel[d];
385  Pm[p1].Vel[d] = Pm[p2].Vel[d];
386  Pm[p2].Vel[d] = Temp[d];
387  }
388  int tmp = Ln[p1].NLink;
389  Ln[p1].NLink = Ln[p2].NLink;
390  Ln[p2].NLink = tmp;
391  for(int l=0;l<MAX(Ln[p1].NLink,Ln[p2].NLink);l++){
392  tmp = Ln[p2].Link[l];
393  Ln[p2].Link[l] = Ln[p1].Link[l];
394  Ln[p1].Link[l] = tmp;
395  }
396  tmp = Pm[p1].CId;
397  Pm[p1].CId = Pm[p2].CId;
398  Pm[p2].CId = tmp;
399  tmp = Pm[p1].Typ;
400  Pm[p1].Typ = Pm[p2].Typ;
401  Pm[p2].Typ = tmp;
402  tmp = Pm[p1].Idx;
403  Pm[p1].Idx = Pm[p2].Idx;
404  Pm[p2].Idx = tmp;
405 }
406 double VarData::TwoPartDist2(int p1,int p2,double *DistRel){
407  return TwoPartDist2(Pm[p1].Pos,p2,DistRel);
408 }
409 double VarData::TwoPartDist2(double *Pos,int p2,double *DistRel){
410  for(int d=0;d<3;d++){
411  DistRel[d] = Pos[d] - Pm[p2].Pos[d];
412  DistRel[d] -= floor(DistRel[d]*pInvEdge(d) + .5)*pEdge(d);
413  }
414  return DistRel[3] = (SQR(DistRel[0])+SQR(DistRel[1])+SQR(DistRel[2]));
415 }
416 double VarData::TwoPartDist(double *Pos,int p2,double *DistRel){
417  return DistRel[3] = sqrt(TwoPartDist2(Pos,p2,DistRel));
418 }
419 double VarData::TwoPartDist(int p1,int p2,double *DistRel){
420  return TwoPartDist(Pm[p1].Pos,p2,DistRel);
421 }
422 int VarData::TwoPartDist(int p1,int p2,double *DistRel,double CutOff){
423  return TwoPartDist(Pm[p1].Pos,p2,DistRel,CutOff);
424 }
425 int VarData::TwoPartDist(double *Pos,int p2,double *DistRel,double CutOff){
426  for(int d=0;d<3;d++){
427  DistRel[d] = Pos[d] - Pm[p2].Pos[d];
428  DistRel[d] -= floor(DistRel[d]*pInvEdge(d) + .5)*pEdge(d);
429  }
430  DistRel[3] = (SQR(DistRel[0])+SQR(DistRel[1])+SQR(DistRel[2]));
431  if(DistRel[3] > SQR(CutOff)) return 0;
432  DistRel[3] = sqrt(DistRel[3]);
433  return 1;
434 }
435 void VarData::ShiftBlock(Vettore *Shift,int b){
436  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
437  if(!(Pm[p].CId%2)) continue;
438  for(int d=0;d<3;d++){
439  Pm[p].Pos[d] += Shift->Val(d);
440  Pm[p].Pos[d] -= floor(Pm[p].Pos[d]*pInvEdge(d))*pEdge(d);
441  }
442  }
443 }
444 void VarData::RotateBlock(Vettore *Axis,Vettore *Origin,int b){
445  double Norm = Axis->Normalize();
446  if(Norm <= 0.){printf("Cannot rotate block, the axis is null\n"); return;};
447  Vettore Proj(3);
448  Vettore Pos(3);
449  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
450  for(int d=0;d<3;d++){
451  Pos.Set(Pm[p].Pos[d]-Origin->Val(d),d);
452  }
453  Proj.Copy(&Pos);
454  //Proj.ProjOnAxis(Axis);
455  Proj.PerpTo(&Pos,Axis);
456  for(int d=0;d<3;d++){
457  Pm[p].Pos[d] = Origin->Val(d) + Pos[d] + 2.*Proj[d];
458  Pm[p].Pos[d] -= floor(Pm[p].Pos[d]*pInvEdge(d))*pEdge(d);
459  }
460  }
461 }
462 void VarData::MirrorBlock(Vettore *S1,Vettore *S2,Vettore *S3,int b){
463  Vettore PS(3);
464  Vettore Pos(3);
465  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
466  for(int d=0;d<3;d++){
467  Pos.Set(Pm[p].Pos[d],0);
468  }
469  PS.ProjOnSurf(S1,S2,S3,&Pos);
470  for(int d=0;d<3;d++){
471  Pm[p].Pos[d] = PS.Val(d);
472  }
473  }
474 }
475 void VarData::Transform(int b){
476  if(b>=pNBlock())return;
477  Vettore Axis(1.,0.,1.);
478  Axis.Normalize();
479  Vettore Shift(0.,0.,-.6*pEdge(2));
480  Vettore Px1(.5*pEdge(0),.5*pEdge(1),.1*pEdge(2));
481  Vettore Px2(.5*pEdge(0),.5*pEdge(1),.3*pEdge(2));
482  Vettore Px3(.5*pEdge(0),.5*pEdge(1),.7*pEdge(2));
483  int nNano = b-1;
484  Vettore Origin(pNanoPos(nNano,0),pNanoPos(nNano,1),pNanoPos(nNano,2));
485  Origin.Print();
486  //ShiftBlock(&Shift,b);
487  RotateBlock(&Axis,&Origin,b);
488  //MirrorBlock(&Px1,&Px2,&Px3,b);
489 }
490 void VarData::Point2Shape(int iShape){
491  if(VAR_IF_TYPE(iShape,SHAPE_NONE))
493  else if(VAR_IF_TYPE(iShape,SHAPE_SPH))
495  else if(VAR_IF_TYPE(iShape,SHAPE_CYL))
496  //Nano_Dist = &VarData::FieldCyl;
498  else if(VAR_IF_TYPE(iShape,SHAPE_CLUSTER))
500  else if(VAR_IF_TYPE(iShape,SHAPE_PILL))
502  else if(VAR_IF_TYPE(iShape,SHAPE_TILT))
503  //Nano_Dist = &VarData::FieldSphere;
505  else if(VAR_IF_TYPE(iShape,SHAPE_WALL))
507  else if(VAR_IF_TYPE(iShape,SHAPE_PORE))
509  else if(VAR_IF_TYPE(iShape,SHAPE_EXT))
511  else if(VAR_IF_TYPE(iShape,SHAPE_JANUS))
513  else if(VAR_IF_TYPE(iShape,SHAPE_STALK))
515  else if(VAR_IF_TYPE(iShape,SHAPE_TIP))
516  //Nano_Dist = &VarData::FieldParab;
518  else if(VAR_IF_TYPE(iShape,SHAPE_TORUS))
520  else if(VAR_IF_TYPE(iShape,SHAPE_HARM))
522  else if(VAR_IF_TYPE(iShape,SHAPE_UMBR))
524  else if(VAR_IF_TYPE(iShape,SHAPE_BOUND))
526  else{
528  //printf("No distance fuction assigned to the shape of the nano shape %d\n",iShape);
529  }
530 }
531 double VarData::NanoDist2(double x,double y,double z,int n){
532  double Pos[3] = {x,y,z};
533  return NanoDist2(Pos,n);
534 }
535 double VarData::FieldNo(double *Pos,int n){
536  return 100000.;
537 }
538 double VarData::FieldSphere(double *Pos,int n){
539  return (SQR(pNanoPos(n,0)-Pos[0])+SQR(pNanoPos(n,1)-Pos[1])+SQR(pNanoPos(n,2)-Pos[2]));
540 }
541 double VarData::FieldElips(double *Pos,int n){
542  double PosRel[3];
543  for(int d=0;d<3;d++){
544  PosRel[d] = pNanoPos(n,d) - Pos[d];
545  PosRel[d] -= floor(PosRel[d]*pInvEdge(d) + .5)*pEdge(d);
546  }
547  double Ratio = Nano[n].Height>0.?Nano[n].Rad/Nano[n].Height:.5;
548  double Lat = SQR(PosRel[0])+SQR(PosRel[1]);
549  double Norm = SQR(Ratio*(PosRel[2]));
550  return (Lat+Norm);
551 }
552 double VarData::FieldParab(double *Pos,int n){
553  double PosRel[3];
554  for(int d=0;d<3;d++){
555  PosRel[d] = pNanoPos(n,d) - Pos[d];
556  //PosRel[d] -= floor(PosRel[d]*pInvEdge(d) + .5)*pEdge(d);
557  }
558  double r2 = SQR(PosRel[0])+SQR(PosRel[1]);
559  return SQR(-4.*Nano[n].Height*Pos[CNorm]+r2+3.*Nano[n].Height+4.*Nano[n].Height*Nano[n].Pos[CNorm]);
560  double r = sqrt(r2);
561  double z = Pos[CNorm];//Nano[n].Pos[CNorm] - Pos[CNorm];
562  double a = -4.*Nano[n].Height;
563  double b = 1.;
564  double c = -2.*0.;
565  double d = 3.*Nano[n].Height+4.*Nano[n].Height*Nano[n].Pos[CNorm];
566  return SQR(a*z + b*r2 + c*r + d);
567 }
568 double VarData::FieldCyl(double *Pos,int n){
569  double Dist2 = SQR(pNanoPos(n,0)-Pos[0])+SQR(pNanoPos(n,1)-Pos[1]);
570  if(Pos[2] > Nano[n].Height*.5+pNanoPos(n,2))
571  Dist2 += SQR(Pos[2]-Nano[n].Height*.5-pNanoPos(n,2));
572  else if(Pos[2] < -Nano[n].Height*.5+Nano[n].Pos[2])
573  Dist2 += SQR(Pos[2]+Nano[n].Height*.5+pNanoPos(n,2));
574  return (Dist2);
575 }
576 double VarData::FieldTransMem(double *Pos,int n){
577  double PosRel[3];
578  for(int d=0;d<3;d++){
579  PosRel[d] = pNanoPos(n,d) - Pos[d];
580  PosRel[d] -= floor(PosRel[d]*pInvEdge(d) + .5)*pEdge(d);
581  }
582  // double RadLat = SQR(PosRel[CLat1]) + SQR(PosRel[CLat2]);
583  // if(PosRel[CNorm] > .5*Nano[n].Height + 1.) return 500.;
584  // if(PosRel[CNorm] < -.5*Nano[n].Height - 1.) return 500.;
585  // if(PosRel[CNorm] > .5*Nano[n].Height){
586  // PosRel[CNorm] -= .5*Nano[n].Height - .5*Nano[n].Rad;
587  // if(SQR(PosRel[CNorm]) > RadLat )
588  // return - RadLat + SQR(PosRel[CNorm]);
589  // else
590  // return RadLat - SQR(PosRel[CNorm]);
591  // }
592  // if(PosRel[CNorm] < -.5*Nano[n].Height){
593  // PosRel[CNorm] += .5*Nano[n].Height - .5*Nano[n].Rad;
594  // if(SQR(PosRel[CNorm]) > RadLat )
595  // return - RadLat + SQR(PosRel[CNorm]);
596  // else
597  // return RadLat - SQR(PosRel[CNorm]);
598  // }
599  // PosRel[CNorm] = 0.;
600  // return SQR(PosRel[CLat1]) + SQR(PosRel[CLat2]) + SQR(PosRel[CNorm]);
601  // if and how much above the cylinder (for the smoothing)
602  if(PosRel[CNorm] > Nano[n].Height*.5){
603  PosRel[CNorm] = (PosRel[CNorm] - (Nano[n].Height*.5+Nano[n].Rad));
604  for(int d=0;d<3;d++)
605  PosRel[d] *= .65;//*Nano[n].Rad;
606  if(PosRel[CNorm] > 0.) PosRel[CNorm] = 100.;
607  }
608  else if(PosRel[CNorm] < -Nano[n].Height*.5){
609  PosRel[CNorm] = (PosRel[CNorm] + (Nano[n].Height*.5+Nano[n].Rad));
610  for(int d=0;d<3;d++)
611  PosRel[d] *= .65;//*Nano[n].Rad;
612  if(PosRel[CNorm] < 0.) PosRel[CNorm] = 100.;
613  }
614  else{
615  PosRel[CNorm] = 0.;
616  }
617  double Dist2 = SQR(PosRel[0]) + SQR(PosRel[1]) + SQR(PosRel[2]);
618  return (Dist2);
619 }
620 double VarData::FieldTorus(double *Pos,int n){
621  double PosRel[3];
622  for(int d=0;d<3;d++){
623  PosRel[d] = pNanoPos(n,d) - Pos[d];
624  PosRel[d] -= floor(PosRel[d]*pInvEdge(d) + .5)*pEdge(d);
625  }
626  double Radxy = sqrt(SQR(PosRel[CLat1]) + SQR(PosRel[CLat2]));
627  double Temp = SQR(Nano[n].Height - Radxy) - SQR(Nano[n].Rad) + SQR(PosRel[CNorm]);
628  return SQR(Temp);
629 }
630 double VarData::FieldTilt(double *Pos,int n){
631  Vettore NanoAxis(Nano[n].Axis[0],Nano[n].Axis[1],Nano[n].Axis[2]);
632  Vettore PosRel(3);
633  Vettore Dist(3);
634  double dr[3];
635  for(int d=0;d<3;d++){
636  PosRel.Set(pNanoPos(n,d) - Pos[d],d);
637  }
638  double r2 = fabs(Dist.PerpTo3(&PosRel,&NanoAxis));
639  double HeiOnAxis = PosRel.ProjOnAxis(&NanoAxis);
640  if(fabs(HeiOnAxis) > .5*Nano[n].Height){
641  double Sign = HeiOnAxis > 0. ? 1. : -1.;
642  r2 = 0.;
643  for(int d=0;d<3;d++){
644  dr[d] = - (Sign*Nano[n].Height*.5*Nano[n].Axis[d]) + PosRel[d];
645  r2 += SQR(dr[d]);
646  }
647  }
648  return r2;
649 }
650 double VarData::FieldTiltWall(double *Pos,int n){
651  double a = Nano[n].Axis[0];
652  double b = Nano[n].Axis[1];
653  double c = Nano[n].Axis[2];
654  double d = -Nano[n].Axis[0]*Nano[n].Pos[0] - Nano[n].Axis[1]*Nano[n].Pos[1] - Nano[n].Axis[2]*Nano[n].Pos[2];
655  double Dist2 = a*Pos[0]+b*Pos[1]+c*Pos[2]+d;
656  Dist2 = SQR(Dist2)/(SQR(a)+SQR(b)+SQR(c));
657  return Dist2;
658 }
659 double VarData::FieldBound(double *Pos,int n){
660  double dr[3] = {Pos[0],Pos[1],Pos[2]};
661  dr[CNorm] = 0.;
662  if(Pos[CLat1] > .5*pEdge(CLat1)) dr[CLat1] = pEdge(CLat1) - Pos[CLat1];
663  if(Pos[CLat2] > .5*pEdge(CLat2)) dr[CLat2] = pEdge(CLat2) - Pos[CLat2];
664  if(dr[CLat1] < dr[CLat2] ) dr[CLat2] = 0.;
665  else dr[CLat1] = 0.;
666  return SQR(dr[0]) + SQR(dr[1]) + SQR(dr[2]);
667  // return MIN(SQR(dr[CLat1]),SQR(dr[CLat2]));
668 }
669 double VarData::FieldJanus(double *Pos,int n){
670  Vettore NanoAxis(Nano[n].Axis[0],Nano[n].Axis[1],Nano[n].Axis[2]);
671  Vettore PosRel(3);
672  Vettore Dist(3);
673  for(int d=0;d<3;d++){
674  PosRel.Set(pNanoPos(n,d) - Pos[d],d);
675  }
676  PosRel.Set(PosRel.Val(2)*2.,2);
677  double r2 = fabs(Dist.PerpTo3(&PosRel,&NanoAxis));
678  return SQR(r2);
679 }
680 
681 // # include "../Matematica/cvt.H"
682 
683 // int VarData::Voronoi(){
684 // int NDim = 2;
685 // int NCell = Gen->NChain;
686 // int Batch = 1000;
687 // int DontCreate = 4;
688 // int SampleType = 1;//Halton
689 // int NSample = 10000;
690 // int NItMax = 40;
691 // int NItFix = 1;
692 // int Seme = 45322643;
693 // int NIt = 0;
694 // double ItDiff = 0.;
695 // double Energy = 0.;
696 // double *Pos = (double *)calloc(NCell*NDim,sizeof(double));
697 // for(int c=0;c<NCell;c++){
698 // for(int d=0;d<NDim;d++){
699 // Pos[c*NDim+d] = Ch[c].Pos[d];
700 // }
701 // }
702 // cvt(NDim,NCell,Batch,DontCreate,SampleType,NSample,NItMax,NItFix,&Seme,Pos,&NIt,&ItDiff,&Energy);
703 // }
704 
705 // #ifndef VOROPP
706 // #include "../share/voro++/src/voro++.cc"
707 // #include "../../share/include/container.hh"
708 
709 // int VarData::Voronoi(){
710 // double xMin = 0.;
711 // double xMax = Gen->Edge[CLat1];
712 // double yMin = 0.;
713 // double yMax = Gen->Edge[CLat2];
714 // double zMin = .4*Gen->Edge[CNorm];
715 // double zMax = .6*Gen->Edge[CNorm];
716 // double LengthScale = 1./(xMax*2.6);
717 // int nxf = (int)( (xMax-xMin)*LengthScale ) + 1;
718 // int nyf = (int)( (yMax-yMin)*LengthScale ) + 1;
719 // int nzf = (int)( (zMax-zMin)*LengthScale ) + 1;
720 // bool xperiodic = true;
721 // bool yperiodic = true;
722 // bool zperiodic = false;
723 // const int memory=8;
724 // container con(xMin,xMax,yMin,yMax,zMin,zMax,nxf,nyf,nzf,xperiodic,yperiodic,zperiodic,memory);
725 // for(int c=0;c<Gen->NChain;c++){
726 // if(!CHAIN_IF_TYPE(Ch[c].Type,NChType)) continue;
727 // con.put(0,Ch[c].Pos[CLat1],Ch[c].Pos[CLat2],Ch[c].Pos[CNorm]);
728 // }
729 // //con.draw_cells_gnuplot("random_points_v.gnu");
730 // //con.draw_cells_pov("pack_six_cube_v.pov");
731 // //con.draw_particles_pov("pack_six_cube_p.pov");
732 // return 0;
733 // }
734 // #else
735 // int VarData::Voronoi(){
736 // printf("Voronoi library not supplied\n");
737 // return 1;
738 // }
739 // #endif
int RadDistr(int NSample, double *Plot, double Border[2], int How)
rzd representation of the system referring to
Definition: VarDataEl.cpp:24
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
int CId
Chain Identifier.
Definition: VarData.h:224
double FieldParab(double *Pos, int n)
Scalar field of a elipsoid.
Definition: VarDataEl.cpp:552
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
void Copy(Vettore *o)
Copy the vector.
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
NANO * Nano
Extra particle.
Definition: VarData.h:1044
double FieldTiltWall(double *Pos, int n)
Scalar field of a tilted cylinder.
Definition: VarDataEl.cpp:650
double Normalize()
Normlizes a Vettore.
double TwoPartDist2(int p1, int p2, double *RelDist)
Return the relative distance between two particles (wrapped)
Definition: VarDataEl.cpp:406
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
double Min
Minimum value.
double FieldJanus(double *Pos, int n)
Scalar field of a janus peptide.
Definition: VarDataEl.cpp:669
double FieldElips(double *Pos, int n)
Scalar field of a elipsoid.
Definition: VarDataEl.cpp:541
int Asym
Diblock limit of the chain.
Definition: VarData.h:273
double Pre[3]
Three components of the pressure.
Definition: VarData.h:317
double Temp
Temperature.
Definition: VarData.h:303
double Rad
Size.
Definition: VarData.h:445
double FieldNo(double *Pos, int n)
No field.
Definition: VarDataEl.cpp:535
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
Properties SysProperties()
Calculate some basis properties.
Definition: VarDataEl.cpp:234
double Height
Height of the cylinder.
Definition: VarData.h:449
int NChType
Type of chain selected.
Definition: VarData.h:1080
MOMENTI SampleSurfaceMem(int NSample)
Allocate and fill PlotMem with the particle average position.
Definition: VarDataEl.cpp:49
int pType(int p)
Return the type.
void ShiftBlock(Vettore *Shift, int b)
Shift a block wrt to Shift.
Definition: VarDataEl.cpp:435
int CLat2
lateral coordinate
Definition: VarData.h:1078
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
int IfPlotMem
If PlotMem is allocated and filled.
Definition: VarData.h:1094
double FieldTransMem(double *Pos, int n)
Scalar field of a transmembrane protein.
Definition: VarDataEl.cpp:576
void Print()
Prints the components.
NANO_DIST Nano_Dist
Pointer to a distance/field function.
Definition: VarData.h:762
double PerpTo3(Vettore *Pos, Vettore *Axis)
The vector perpendicolar in three dimension (faster)
double * PlotMem
Particle position/density on the square lattice.
Definition: VarData.h:1056
Some calculated properties of the system.
Definition: VarData.h:364
double Val(int N)
Value of the N column.
double m
y = m*x + q
double FieldCyl(double *Pos, int n)
Scalar field of a cylinder.
Definition: VarDataEl.cpp:568
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
Linear interpolation.
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 Idx
Particle identifier.
Definition: VarData.h:222
Moments of a distribution.
void Transpose()
Transpose the matrix.
double TwoPartDist(int p1, int p2, double *RelDist)
Return the relative distance between two particles (wrapped)
Definition: VarDataEl.cpp:419
int Folding()
Boh...
Definition: VarDataEl.cpp:302
int pNBlock()
Number of blocks.
void RotateBlock(Vettore *Axis, Vettore *Origin, int b)
Rotate a block wrt to the Axis from the Origin.
Definition: VarDataEl.cpp:444
void MirrorBlock(Vettore *Px1, Vettore *Px2, Vettore *Px3, int b)
Mirror the position wrt to a plane.
Definition: VarDataEl.cpp:462
int NChain
Number of chain.
Definition: VarData.h:347
int NEdge
Number of particles per edge.
Definition: VarData.h:1084
double ProjOnAxis(Vettore *a)
Projects along the axis.
double Val(int row)
Returns a value in 1d.
void Point2Shape(int iShape)
Point to the shape function.
Definition: VarDataEl.cpp:490
void ChangeNChain(int NChain, int b)
Update the new number of chains.
Definition: VarDataEl.cpp:340
int NPCh
Number of particle per chain.
Definition: VarData.h:349
bool Add(int row, int col, double Val)
Add the value of the coefficient to the previous one.
void SwapChain(int c1, int c2, int b)
Swap two chains.
Definition: VarDataEl.cpp:355
MOMENTI SampleSurfacePart(double *Plot, int NSample, int Type)
Define a normal coordinate for every patch.
Definition: VarDataEl.cpp:178
double Pos[4]
xyzr Postion of the chain
Definition: VarData.h:238
int NPCh
particles per chain
Definition: VarData.h:259
Matrice computes the algebric operations on matrices.
bool Set(int row, int column, double Val)
Set a coefficient.
double FieldSphere(double *Pos, int n)
Scalar field of a sphere.
Definition: VarDataEl.cpp:538
int Size()
Size of Matrice.
double Max
Maximum value.
double Uno
First moment.
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int NChain
chains
Definition: VarData.h:263
int EndIdx
End particle position.
Definition: VarData.h:267
int ApplyFilter(Matrice *Point, Matrice *Res, Matrice *Mask)
Applies the filter.
double pPos(int p, int d)
Return back folded position.
double RePhob
End2End distance of the.
Definition: VarData.h:373
void Transform(int block)
Transform a block.
Definition: VarDataEl.cpp:475
double FieldTilt(double *Pos, int n)
Scalar field of a tilted cylinder.
Definition: VarDataEl.cpp:630
int CLat1
lateral coordinate
Definition: VarData.h:1076
int SpatialDerivative(Matrice *Surface, Matrice *Resp, SPLINE Weight, int NSample)
Perform a spatial derivative on a surface.
Definition: VarDataEl.cpp:222
int Num
Number of points considered.
double ProjOnSurf(Vettore *S1, Vettore *S2, Vettore *S3, Vettore *P)
Project a point P on the point PS perpendicular to the surface described by the points S1...
Coefficient of a spline.
char * SysState()
Calculates fundamental quantities.
Definition: VarDataEl.cpp:27
double Axis[3]
Rotation axis.
Definition: VarData.h:435
void Print()
Print the current values.
Definition: VarData.h:369
double FieldTorus(double *Pos, int n)
Scalar field of a janus peptide.
Definition: VarDataEl.cpp:620
int NPart
Number of particle.
Definition: VarData.h:345
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int InterBSpline2D(double **PlIn, double **PmOut, int NIn, int NOut)
2-d BSpline
int Typ
Type.
Definition: VarData.h:226
int NBlock
Number of blocks.
Definition: VarData.h:359
double FieldBound(double *Pos, int n)
Scalar field of a hard wall at the box edges.
Definition: VarDataEl.cpp:659
void SampleSurface(double *Plot, int NSample, int Type)
Define a normal coordinate for every patch.
Definition: VarDataEl.cpp:65
double SurfTens
Surface tension.
Definition: VarData.h:321
void Set(double Val, int Col)
Set the N column.
void PerpTo(Vettore *o)
The vector perpendicolar.
int CNorm
Normal coordinate.
Definition: VarData.h:1074
int pNPart()
Number of particle.
RETTA InterRett(double *Px, double *Py, int NMass)
Linear interpolation.
void LoadDensFile(double **Plot, int NBin)
Load in the array Plot the density of the system.
Definition: VarDataEl.cpp:150