Allink  v0.1
VarDataProfile.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 ;
25 int VarData::DensityProfile(int coord,int NBin,int NType,double *DensityAv){
26  BfDefChain();
27  int NPartition = 5;
28  double *Density = (double *)calloc(NPartition*NPartition*NType*NBin,sizeof(double));
29  double Pos[4];
30  char cLine[STRSIZE];
31  for(int b=0;b<Gen->NBlock;b++){
32  if(!strncmp(Block[b].Name,"PEP",3) ) continue;
33  if(!strcmp(Block[b].Name,"WATER") ) continue;
34  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
35  int t = Pm[p].Typ;
36  for(int d=0;d<3;d++)
37  Pos[d] = Pm[p].Pos[d];
38  int vx = (int)(NPartition*Pos[CLat1]/Gen->Edge[CLat2]);
39  if(vx < 0 || vx >= NBin) continue;
40  int vy = (int)(NPartition*Pos[CLat2]/Gen->Edge[CLat2]);
41  if(vy < 0 || vy >= NBin) continue;
42  int vz = (int)(NBin*Pos[CNorm]/Gen->Edge[CNorm]);
43  if(vz < 0 || vz >= NBin) continue;
44  if(CHAIN_IF_TYPE(Ch[Pm[p].CId].Type,CHAIN_ADDED))
45  t = 2;
46  Density[((vz*NType+t)*NPartition+vy)*NPartition+vx] += 1.;
47  }
48  }
49  for(int vx=0;vx<NPartition;vx++){
50  for(int vy=0;vy<NPartition;vy++){
51  double Media = 0.;
52  double Weight = 0.;
53  for(int v=0;v<NBin;v++){
54  Weight += Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
55  Media += v*Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
56  }
57  int vMedia = Weight > 0 ? (int)(Media/Weight) : NBin/2;
58  int vOffset = vMedia - NBin/2;
59  for(int t=0;t<NType;t++){
60  for(int v=0;v<NBin;v++){
61  int vOther = v + vOffset;
62  if(vOther < 0)
63  vOther = NBin - vOther;
64  if(vOther >= NBin)
65  vOther = vOther - NBin;
66  DensityAv[v*NType+t] += Density[((vOther*NType+t)*NPartition+vy)*NPartition+vx];
67  }
68  }
69  }
70  }
71  free(Density);
72  return 0;
73 }
74 // Volume element in cylindrical symmetry
75 void VarData::VolumeCircSlab(double *VolContr,int NBin){
76  if(pEdge(3) <= 0.){printf("Lateral radius not defined\n");};
77  double LenMin = .5*MIN(Gen->Edge[CLat1],Gen->Edge[CLat2]);
78  double LenMax = .5*MAX(Gen->Edge[CLat1],Gen->Edge[CLat2]);
79  double RapCorner = LenMin/LenMax*NBin;
80  double RapMin = LenMin/Gen->Edge[3]*NBin;
81  double Volume2 = 0.;
82  for(int v=0;v<NBin;v++){
83  double Rad2 = (v+1);
84  double Volume1 = PI*QUAD(Rad2);
85  if(Rad2 > RapMin){
86  double Pro2 = sqrt(QUAD(v+1)-QUAD(RapMin));
87  double Angle2 = atan(RapMin/Pro2);
88  Volume1 = 2.*QUAD(Rad2)*Angle2 + 2.*Pro2*RapMin;
89  if(Rad2 > RapCorner){
90  double theta = acos(RapCorner/Rad2);
91  double b = Rad2*sin(theta);
92  //double b = sqrt( SQR(Rad2) - SQR(RapCorner));
93  Volume1 -= 4.*(Rad2*Rad2*theta - b*RapCorner);
94  //printf("%lf %lf %lf %lf %lf %lf %lf\n",Rad2,RapCorner,theta,b,Volume1-Volume2,Rad2*Rad2*theta ,b*RapCorner);
95  }
96  }
97  double Volume = (Volume1 - Volume2)*Gen->Edge[CNorm]*SQR(Gen->Edge[3])/(double)CUB(NBin);
98  VolContr[v] = Volume;
99  Volume2 = Volume1;
100  }
101 }
102 int VarData::Worm(int Partition,int NBin,double *Border,double *dPoint){
103  int cLong = CLat2;
104  int cLat = CLat1;
105  double *Cm;
106  RETTA *r1;
107  r1 = (RETTA *)calloc(Partition,sizeof(RETTA));
108  Cm = (double *)calloc(3*Partition,sizeof(double));
109  double* Dis0 = (double *)calloc(Gen->NPart,sizeof(double));
110  double* Dis1 = (double *)calloc(Gen->NPart,sizeof(double));
111  for(int pa = 0,cc=0;pa<Partition;pa++){
112  for(int c=0;c<Gen->NChain;c++){
113  int v = (int)(Ch[c].Pos[cLong]*Partition/Gen->Edge[cLong]);
114  if(v != pa) continue;
115  *(Dis0 + cc) = Ch[c].Pos[cLat];
116  *(Dis1 + cc) = Ch[c].Pos[cLong];
117  //printf("%lf %lf %lf\n",Ch[c].Pos[0],Ch[c].Pos[1],Ch[c].Pos[2]);
118  for(int d=0;d<3;d++)
119  *(Cm + 3*v + d) += Ch[c].Pos[d];
120  cc++;
121  }
122  r1[pa] = Mat->InterRett(Dis0,Dis1,cc);
123  for(int d=0;d<3;d++)
124  *(Cm + 3*pa + d) /= cc;
125  cc = 0;
126  }
127  double pos0;
128  double pos1;
129  double qcm;
130  double qpoint = 0.;
131  double dist;
132  int pp0 = 0;
133  int pp1 = 0;
134  //FILE *Ciccia = fopen("Ciccia.dat","w");
135  for(int pa=0;pa<Partition;pa++){
136  //printf("m %lf Cm %lf %lf %lf \n",r1[pa].m,*(Cm+3*pa+0),*(Cm+3*pa+1),*(Cm+3*pa+2));
137  for(int p=0;p<Gen->NPart;p++){
138  int v = (int) (Pm[p].Pos[cLong]*Partition/Gen->Edge[cLong]);
139  if(v != pa) continue;
140  qcm = *(Cm + 3*pa + cLat) + *(Cm + 3*pa + cLong)/r1[pa].m;
141  qpoint = Pm[p].Pos[cLat] - r1[pa].m*Pm[p].Pos[cLong];
142  pos0 = (qcm-qpoint)*r1[pa].m/(QUAD(r1[pa].m) + 1.);
143  pos1 = pos0*r1[pa].m + qpoint;
144  dist = QUAD((pos0-Cm[3*pa+cLong])) + QUAD((pos1-Cm[3*pa+cLat]));
145  if(Pm[p].Typ == 0){
146  Dis0[pp0] = pow(dist,.5);
147  if(Pm[p].Pos[cLat] < *(Cm+3*pa+cLat))
148  Dis0[pp0] *= -1.;
149  //printf("Dis0 %lf\n",Dis0[pp0]);
150  pp0++;
151  }
152  else if(Pm[p].Typ == 1){
153  Dis1[pp1] = pow(dist,.5);
154  if(Pm[p].Pos[cLat] < *(Cm+3*pa+cLat))
155  Dis1[pp1] *= -1.;
156  pp1++;
157  }
158  if(pa == 2){
159  //fprintf(Ciccia,"Ch %d %lf %lf Pos %lf %lf pos %lf %lf Dis1 %lf Dis2 %lf dist %lf\n",Pm[p].CId,Ch[Pm[p].CId].Pos[0],Ch[Pm[p].CId].Pos[2],Pm[p].Pos[0],Pm[p].Pos[2],pos0,pos1,Dis0[pp0],Dis1[pp1],pow(dist,.5));
160  }
161  }
162  // fclose(Ciccia);
163  for(int v=0;v<NBin;v++)
164  dPoint[v*2] = 0.;
165  for(int i=0;i<pp0;i++){
166  int v = (int) ((Dis0[i]-Border[0])/(Border[1] - Border[0])*NBin);
167  if(v>= NBin) continue;
168  dPoint[v*2] += 1.;
169  }
170  for(int v=0;v<NBin;v++){
171  dPoint[v*2] *= (NBin*Partition)/(Gen->Edge[0]*Gen->Edge[1]*Gen->Edge[2]);
172  }
173  for(int v=0;v<NBin;v++)
174  dPoint[v*2] = 0.;
175  for(int i=0;i<pp1;i++){
176  int v = (int) ((Dis1[i]-Border[0])/(Border[1] - Border[0])*NBin);
177  if(v>= NBin) continue;
178  dPoint[v*2] += 1.;
179  }
180  for(int v=0;v<NBin;v++){
181  dPoint[v*2+1] /= (double) (NBin*Partition)/(Gen->Edge[0]*Gen->Edge[1]*Gen->Edge[2]);
182  //if(pa==0)printf("%d %lf %lf\n",v,dPoint[v],dPoint1[v]);
183  }
184  pp0=0;
185  pp1=0;
186  }
187  for(int v=0;v<NBin;v++){
188  dPoint[v*2] /= (double) NBin;
189  dPoint[v*2+1] /= (double)NBin;
190  }
191  return 1;
192 }
193 void VarData::Stalk(int NBin,int NLevel,double **Plot,double Threshold){
194  double **Norma = (double **)calloc(NLevel,sizeof(double));
195  for(int l=0;l<NLevel;l++)
196  Norma[l] = (double *)calloc(NBin*NBin,sizeof(double));
197  for(int b=0;b<Gen->NBlock;b++){
198  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
199  for(int v=0;v<NBin;v++){
200  // if(Pm[p].Typ == 1 && Pm[p-1].Typ == 0)
201  {
202  int v = (int)(Pm[p].Pos[CLat1]/Gen->Edge[CLat1]*NBin);
203  if( v < 0 || v >= NBin) continue;
204  int vv = (int)(Pm[p].Pos[CLat2]/Gen->Edge[CLat2]*NBin);
205  if( vv < 0 || vv >= NBin) continue;
206  Plot[b][v*NBin + vv] += Pm[p].Pos[CNorm];
207  Norma[b][v*NBin + vv] += 1.;
208  }
209  }
210  }
211  }
212  for(int v=0;v<NBin;v++)
213  for(int vv=0;vv<NBin;vv++)
214  for(int l=0;l<NLevel;l++){
215  if(Norma[l][v*NBin+vv] > Threshold)
216  Plot[l][v*NBin+vv] /= Norma[l][v*NBin+vv];
217  else
218  Plot[l][v*NBin+vv] = 0.;
219  }
220  for(int l=0;l<NLevel;l++)
221  free(Norma[l]);
222  free(Norma);
223 }
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
BLOCK * Block
Information for every block.
Definition: VarData.h:1054
double Edge[4]
xyzr edges of the simulation box
Definition: VarData.h:309
int BfDefChain()
Definition of the chain.
int CLat2
lateral coordinate
Definition: VarData.h:1078
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
double m
y = m*x + q
Linear interpolation.
int DensityProfile(int coord, int NSample, int NType, double *dDensity)
Calculate the density profile for the x, y, z, r coordinate.
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
int NChain
Number of chain.
Definition: VarData.h:347
int Worm(int Partition, int NSample, double *Border, double *dPoint)
Density profile along a worm like micelle.
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
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int EndIdx
End particle position.
Definition: VarData.h:267
int CLat1
lateral coordinate
Definition: VarData.h:1076
int NPart
Number of particle.
Definition: VarData.h:345
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int Typ
Type.
Definition: VarData.h:226
int NBlock
Number of blocks.
Definition: VarData.h:359
void VolumeCircSlab(double *VolContr, int NSample)
Fill an array of.
int CNorm
Normal coordinate.
Definition: VarData.h:1074
RETTA InterRett(double *Px, double *Py, int NMass)
Linear interpolation.