Allink  v0.1
VarDataExp.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::PairCorrelation(double *dPoint,int NSample,int How,int Type){
25  BfDefChain();
26  //printf("%d %d\n",NSample,How);
27  double dNSample = 1./(double)NSample;
28  for(int v=0;v<NSample;v++)
29  dPoint[v] = 0.;
30  if(How == 0){//Monomer
31  for(int p=0;p<Gen->NPart;p++){
32  if(Pm[p].Typ != Type)continue;
33  Pm[p].Pos[3] = 0.;
34  for(int pp=0;pp<Gen->NPart;pp++){
35  if(Pm[pp].Typ != Type)continue;
36  if( p == pp)continue;
37  double Dist2 = QUAD(Pm[p].Pos[0] - Pm[pp].Pos[0]);
38  Dist2 += QUAD(Pm[p].Pos[1] - Pm[pp].Pos[1]);
39  Dist2 += QUAD(Pm[p].Pos[2] - Pm[pp].Pos[2]);
40  Dist2 = pow(Dist2,.5);
41  int v = (int)(Dist2*pInvEdge(3)*NSample);
42  if( v < 0 || v >= NSample) continue;
43  dPoint[v] += 1.;
44  }
45  }
46  // printf("Processing: %d/%d %.1f %% \r",p,HowMany,p/(double)HowMany/100.);
47  for(int v=0;v<NSample;v++){
48  dPoint[v] /= Gen->NPart*(Gen->NPart-1);
49  }
50  }
51  //printf("\n");
52  if(How == 1){//Chain
53  for(int c=0;c<Gen->NChain;c++){
54  Ch[c].Pos[3] = 0.;
55  if(Ch[c].Type != Type)continue;
56  for(int cc=0;cc<Gen->NChain;cc++){
57  if(Ch[cc].Type != Type)continue;
58  if( c == cc)continue;
59  Ch[c].Pos[3] += QUAD((Ch[c].Pos[CLat1] - Ch[cc].Pos[CLat1]));
60  Ch[c].Pos[3] += QUAD((Ch[c].Pos[CLat2] - Ch[cc].Pos[CLat2]));
61  }
62  Ch[c].Pos[3] = pow(Ch[c].Pos[3]/(double)Gen->NChain,.5);
63  int v = (int)(Ch[c].Pos[3] / Gen->Edge[3]*NSample);
64  if( v < 0 || v >= NSample) continue;
65  dPoint[v] += 1.;
66  }
67  }
68  for(int v=0;v<NSample;v++){
69  dPoint[v] /= DUE_PI*( QUAD((Gen->Edge[3]*(v+1)*dNSample)) - QUAD((Gen->Edge[3]*v*dNSample)) );
70  }
71  return 0;
72 }
73 int VarData::PairCorrelationSquare(double **dPoint,int NSample,int Type){
74  BfDefChain();
75  double dNSample = 1./(double)NSample;
76  double InvNChain=1./(double)Gen->NChain;
77  int vRef = (int)(NSample/2.);
78  for(int c=0;c<Gen->NChain;c++){
79  if(!CHAIN_IF_TYPE(Ch[c].Type,NChType)) continue;
80  double ChX1 = (Ch[c].Pos[CLat1] - Gen->Edge[CLat1]*.5);
81  double ChY1 = (Ch[c].Pos[CLat2] - Gen->Edge[CLat2]*.5);
82  int vx1 = (int)(ChX1 / (Gen->Edge[CLat1])*NSample);
83  int vy1 = (int)( (ChY1) / (Gen->Edge[CLat2])*NSample);
84  for(int cc=0;cc<Gen->NChain;cc++){
85  if(!CHAIN_IF_TYPE(Ch[cc].Type,NChType)) continue;
86  if( c == cc) continue;
87  double ChX2 = (Ch[cc].Pos[CLat1] - Gen->Edge[CLat1]*.5);
88  double ChY2 = (Ch[cc].Pos[CLat2] - Gen->Edge[CLat2]*.5);
89  int vx2 = (int)(ChX2 / (Gen->Edge[CLat1])*NSample);
90  int vy2 = (int)( (ChY2) / (Gen->Edge[CLat2])*NSample);
91  int vx = vx1 - vx2;
92  // if(vx < -vRef) vx += vRef;
93  // if(vx > vRef ) vx -= vRef;
94  if( vx + vRef < 0 || vx+vRef >= NSample) continue;
95  int vy = vy1 - vy2;
96  // if(vy < -vRef) vy += vRef;
97  // if(vy > vRef) vy -= vRef;
98  if( vy + vRef < 0 || vy+vRef >= NSample) continue;
99  // printf("%d %d %d %lf %lf\n",c,vx,vy,ChX1,ChX2);
100  dPoint[vx+vRef][vy+vRef] += InvNChain;
101  }
102  }
103  return 0;
104 }
105 int VarData::PairCorrelationPep(double **dPoint,int NSample,int Type){
106  BfDefChain();
107  double dNSample = 1./(double)NSample;
108  double InvNChain=1./(double)Gen->NChain;
109  int vRef = (int)(NSample/2.);
110  for(int c=0;c<Gen->NChain;c++){
111  if(!CHAIN_IF_TYPE(Ch[c].Type,NChType))continue;
112  double ChX = remainder(Ch[c].Pos[CLat1] - pNanoPos(0,CLat1),Gen->Edge[CLat1]);
113  double ChY = remainder(Ch[c].Pos[CLat2] - pNanoPos(0,CLat2),Gen->Edge[CLat2]);
114  int vx = (int)(ChX/Gen->Edge[CLat1]*NSample);
115  int vy = (int)(ChY/Gen->Edge[CLat2]*NSample);
116  if( vx + vRef < 0 || vx+vRef >= NSample) continue;
117  if( vy + vRef < 0 || vy+vRef >= NSample) continue;
118  dPoint[vx+vRef][vy+vRef] += InvNChain;
119  }
120  return 0;
121 }
122 int VarData::PairCorrelationRound(double **dPoint,int NSample,int Type){
123  BfDefChain();
124  double dNSample = 1./(double)NSample;
125  double InvNChain=1./(double)Gen->NChain;
126  for(int c=0;c<Gen->NChain;c++){
127  double ChRad=0.;
128  double ChAngle=0.;
129  if(!CHAIN_IF_TYPE(Ch[c].Type,NChType))continue;
130  for(int cc=0;cc<Gen->NChain;cc++){
131  if(!CHAIN_IF_TYPE(Ch[cc].Type,NChType))continue;
132  if( c == cc)continue;
133  double ChX = remainder(Ch[c].Pos[CLat1] - Ch[cc].Pos[CLat1],Gen->Edge[CLat1]);
134  // if(ChX > .5*Gen->Edge[CLat1])
135  // ChX = Ch[c].Pos[CLat1] + Ch[cc].Pos[CLat1] - Gen->Edge[CLat1];
136  // else if (ChX < -.5*Gen->Edge[CLat1])
137  // ChX = -Ch[c].Pos[CLat1] - Ch[cc].Pos[CLat1] + Gen->Edge[CLat1];
138  double ChY = remainder(Ch[c].Pos[CLat2] - Ch[cc].Pos[CLat2],Gen->Edge[CLat2]);
139  // if(ChY > .5*Gen->Edge[CLat2])
140  // ChY = Ch[c].Pos[CLat2] + Ch[cc].Pos[CLat2] - Gen->Edge[CLat2];
141  // else if (ChY < -.5*Gen->Edge[CLat2])
142  // ChY = -Ch[c].Pos[CLat2] - Ch[cc].Pos[CLat2] + Gen->Edge[CLat2];
143  ChRad = sqrt( QUAD((ChX)) + QUAD((ChY)) );
144  ChAngle = acos(ChX / ChRad);
145  if(ChY < 0)
146  ChAngle = DUE_PI - ChAngle;
147  int v = (int)(ChRad / (Gen->Edge[3])*NSample);
148  if( v < 0 || v >= NSample) continue;
149  int vv = (int)( (ChAngle) / (DUE_PI)*NSample);
150  if( vv < 0 || vv >= NSample) continue;
151  //printf("%d %d %lf %lf %lf %lf\n",v,vv,ChRad,ChX,ChY,ChAngle);
152  dPoint[v][vv] += InvNChain;
153  }
154  }
155  return 0;
156 }
157 int VarData::Scattering2d(double **Plot,int NSample,int Type){
158  BfDefChain();
159  double dNSample = 1./(double)NSample;
160  double InvNChain=1./(double)Gen->NChain;
161  int vRef = (int)(NSample/2.);
162  double *CosSin = (double *)calloc(2*SQR(NSample),sizeof(double));
163  for(int p=0;p<pNPart();p++){
164  //if(!CHAIN_IF_TYPE(Ch[Pm[p].CId].Type,NChType))continue;
165  // double ChX = (Ch[c].Pos[CLat1] - Nano->PosBf[CLat1]);
166  // double ChY = (Ch[c].Pos[CLat2] - Nano->PosBf[CLat2]);
167  double ChX = (Pm[p].Pos[CLat1] - pNanoPos(0,CLat1));
168  double ChY = (Pm[p].Pos[CLat2] - pNanoPos(0,CLat2));
169  ChX -= floor(ChX*pInvEdge(CLat1))*pEdge(CLat1);
170  ChY -= floor(ChY*pInvEdge(CLat2))*pEdge(CLat2);
171  ChX *= pInvEdge(CLat1);
172  ChY *= pInvEdge(CLat2);
173  // printf("%lf %lf %lf %lf\n",vvx,ChX,vvx*ChX*DUE_PI,cos(vvx*ChX*DUE_PI));
174  for(int vx=0;vx<NSample;vx++){
175  for(int vy=0;vy<NSample;vy++){
176  double qx = vx*dNSample;
177  double qy = vy*dNSample;
178  CosSin[(vx*NSample+vy)*2 ] += cos( (qx*ChX+qy*ChY)*DUE_PI );
179  CosSin[(vx*NSample+vy)*2+1] += sin( (qx*ChX+qy*ChY)*DUE_PI );
180  }
181  }
182  }
183  for(int vx=0;vx<NSample;vx++){
184  for(int vy=0;vy<NSample;vy++){
185  Plot[vx][vy] = (SQR(CosSin[(vx*NSample+vy)*2 ]) + SQR(CosSin[(vx*NSample+vy)*2+1]))*SQR(InvNChain);
186  }
187  }
188  return 0;
189 }
190 int VarData::Scattering2D(double **dPoint,int NSample,int Type){
191  BfDefChain();
192  double dNSample = 1./(double)NSample;
193  double InvNChain=1./(double)Gen->NChain;
194  for(int c=0;c<Gen->NChain;c++){
195  if(!CHAIN_IF_TYPE(Ch[c].Type,NChType))continue;
196  for(int cc=0;cc<Gen->NChain;cc++){
197  if(!CHAIN_IF_TYPE(Ch[cc].Type,NChType))continue;
198  if( c == cc)continue;
199  double ChX = (Ch[c].Pos[CLat1] - Ch[cc].Pos[CLat1]);
200  if(ChX > .5*Gen->Edge[CLat1]) ChX = Ch[c].Pos[CLat1] + Ch[cc].Pos[CLat1] - Gen->Edge[CLat1];
201  else if (ChX < -.5*Gen->Edge[CLat1]) ChX = -Ch[c].Pos[CLat1] - Ch[cc].Pos[CLat1] + Gen->Edge[CLat1];
202  double ChY = (Ch[c].Pos[CLat2] - Ch[cc].Pos[CLat2]);
203  if(ChY > .5*Gen->Edge[CLat2]) ChY = Ch[c].Pos[CLat2] + Ch[cc].Pos[CLat2] - Gen->Edge[CLat2];
204  else if (ChY < -.5*Gen->Edge[CLat2]) ChY = -Ch[c].Pos[CLat2] - Ch[cc].Pos[CLat2] + Gen->Edge[CLat2];
205  int v = NSample/2 + (int)( 2.*(ChX / Gen->Edge[CLat1])*(NSample/2) );
206  int vv = NSample/2 + (int)(2.*ChY / (Gen->Edge[CLat2])*(NSample/2));
207  //if(ChX > Gen->Edge[CLat1]*.5 || ChY > Gen->Edge[CLat2]*.5 ) printf("%d %d %lf %lf\n",v,vv,ChX,ChY);
208  // printf("%d %d %lf %lf\n",v,vv,ChX,ChY);
209  if( vv < 0 || vv >= NSample) continue;
210  if( v < 0 || v >= NSample) continue;
211  dPoint[vv][v] += InvNChain;//*sin(ChRec*ChRad)/ChRec;
212  }
213  }
214  return 0;
215 }
216 void VarData::Spettro2d(double *Points,int NSample,int Type){
217  double *Plot = (double *)calloc(SQR(NSample),sizeof(double));
218  double *InPoints = (double *)calloc(NSample*NSample,sizeof(double));
219  SampleSurface(Plot,NSample,Type);
220  for(int v=0;v<SQR(NSample);v++){
221  InPoints[v] = Plot[v];
222  }
223  Mat->Spettro2d(InPoints,Points,NSample);
224  free(Plot);
225  free(InPoints);
226 }
227 void VarData::Spettro2d(double *Plot,int NSample){
228  double *Points = (double *)calloc(NSample*NSample,sizeof(double));
229  for(int v=0;v<SQR(NSample);v++){
230  Points[v] = Plot[v];
231  }
232  Mat->Spettro2d(Points,Plot,NSample);
233  free(Points);
234 }
CHAIN * Ch
Information on all chains.
Definition: VarData.h:1050
int PairCorrelationPep(double **Point, int NSample, int Type)
2-d pair correlation on a square fererring to the pep position
Definition: VarDataExp.cpp:105
double Edge[4]
xyzr edges of the simulation box
Definition: VarData.h:309
int BfDefChain()
Definition of the chain.
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
int Scattering2d(double **Point, int NSample, int Type)
2-d Scattering
Definition: VarDataExp.cpp:157
int NChType
Type of chain selected.
Definition: VarData.h:1080
int CLat2
lateral coordinate
Definition: VarData.h:1078
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
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 NChain
Number of chain.
Definition: VarData.h:347
void Spettro2d(double *st, double *sw, int NMass)
Compute the 2d spectrum of.
double Pos[4]
xyzr Postion of the chain
Definition: VarData.h:238
int PairCorrelation(double *Point, int NSample, int How, int Type)
1-d pair correlation
Definition: VarDataExp.cpp:24
void Spettro2d(double *Points, int NSample, int Type)
1-d spectrum of a surface
Definition: VarDataExp.cpp:216
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
int Scattering2D(double **Point, int NSample, int Type)
2-d scattering
Definition: VarDataExp.cpp:190
int CLat1
lateral coordinate
Definition: VarData.h:1076
int PairCorrelationSquare(double **Point, int NSample, int Type)
2-d pair correlation on a square
Definition: VarDataExp.cpp:73
int NPart
Number of particle.
Definition: VarData.h:345
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
void SampleSurface(double *Plot, int NSample, int Type)
Define a normal coordinate for every patch.
Definition: VarDataEl.cpp:65
int PairCorrelationRound(double **Point, int NSample, int Type)
Circular 2-d pair correlation.
Definition: VarDataExp.cpp:122
int pNPart()
Number of particle.