Allink  v0.1
VarDataContour.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::DensityProfile(int coord,int Values,int NType,double *DensityAv){
25  BfDefChain();
26  double VolumeInv = (double) (Values)/(Gen->Edge[0]*Gen->Edge[1]*Gen->Edge[2]);
27  int NPartition = 10;
28  double NPartInv = 1./(double)(NPartition*NPartition);
29  double *Density = (double *)calloc(NPartition*NPartition*NType*Values,sizeof(double));
30  double Pos[4];
31  char cLine[STRSIZE];
32  int t=0;
33  for(int b=0;b<Gen->NBlock;b++){
34  if(!strncmp(Block[b].Name,"PEP",3) ) continue;
35  if(!strcmp(Block[b].Name,"WATER") ) continue;
36  for(int p=Block[b].InitIdx,link=0;p<Block[b].EndIdx;p++){
37  t = Pm[p].Typ;
38  for(int d=0;d<3;d++)
39  Pos[d] = Pm[p].Pos[d];
40  int vx = (int)(NPartition*Pos[CLat1]/Gen->Edge[CLat2]);
41  if(vx < 0 || vx >= Values) continue;
42  int vy = (int)(NPartition*Pos[CLat2]/Gen->Edge[CLat2]);
43  if(vy < 0 || vy >= Values) continue;
44  int vz = (int)(Values*Pos[CNorm]/Gen->Edge[CNorm]);
45  if(vz < 0 || vz >= Values) continue;
46  if(CHAIN_IF_TYPE(Ch[Pm[p].CId].Type,CHAIN_ADDED))
47  t = 2;
48  Density[((vz*NType+t)*NPartition+vy)*NPartition+vx] += 1.;
49  }
50  }
51  for(int vx=0;vx<NPartition;vx++){
52  for(int vy=0;vy<NPartition;vy++){
53  double Media = 0.;
54  double Weight = 0.;
55  for(int v=0;v<Values;v++){
56  Weight += Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
57  Media += v*Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
58  }
59  int vMedia = Weight > 0 ? (int)(Media/Weight) : Values/2;
60  int vOffset = vMedia - Values/2;
61  for(int t=0;t<NType;t++){
62  for(int v=0;v<Values;v++){
63  int vOther = v + vOffset;
64  if(vOther < 0)
65  vOther = Values - vOther;
66  if(vOther >= Values)
67  vOther = vOther - Values;
68  DensityAv[v*NType+t] += Density[((vOther*NType+t)*NPartition+vy)*NPartition+vx]*VolumeInv;
69  }
70  }
71  }
72  }
73  free(Density);
74  return 0;
75 }
76 void VarData::Stalk(int Values,int NLevel,double **Plot,double Threshold){
77  double **Norma = (double **)calloc(NLevel,sizeof(double));
78  for(int l=0;l<NLevel;l++)
79  Norma[l] = (double *)calloc(Values*Values,sizeof(double));
80  for(int b=0;b<Gen->NBlock;b++){
81  for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
82  for(int v=0;v<Values;v++){
83  // if(Pm[p].Typ == 1 && Pm[p-1].Typ == 0)
84  {
85  int v = (int)(Pm[p].Pos[CLat1]/Gen->Edge[CLat1]*Values);
86  if( v < 0 || v >= Values) continue;
87  int vv = (int)(Pm[p].Pos[CLat2]/Gen->Edge[CLat2]*Values);
88  if( vv < 0 || vv >= Values) continue;
89  Plot[b][v*Values + vv] += Pm[p].Pos[CNorm];
90  Norma[b][v*Values + vv] += 1.;
91  }
92  }
93  }
94  }
95  for(int v=0;v<Values;v++)
96  for(int vv=0;vv<Values;vv++)
97  for(int l=0;l<NLevel;l++){
98  if(Norma[l][v*Values+vv] > Threshold)
99  Plot[l][v*Values+vv] /= Norma[l][v*Values+vv];
100  else
101  Plot[l][v*Values+vv] = 0.;
102  }
103  for(int l=0;l<NLevel;l++)
104  free(Norma[l]);
105  free(Norma);
106 }
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
int DensityProfile(int coord, int NSample, int NType, double *dDensity)
Calculate the density profile for the x, y, z, r coordinate.
void Stalk(int NSample, int NLevel, double **Plot, double Threshold)
Following the contour of a stalk.
int EndIdx
End particle position.
Definition: VarData.h:267
int CLat1
lateral coordinate
Definition: VarData.h:1076
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
int CNorm
Normal coordinate.
Definition: VarData.h:1074