Allink: src/VarData/VarDataProfile.cpp Source File

00001 /***********************************************************************
00002 VarDataEl: Elaboration functions for the VarData class. This functions
00003 provides a simple manipulation of the data read by [Open]. The
00004 options are provided to elaborate different system's shapes.
00005 Copyright (C) 2008 by Giovanni Marelli <sabeiro@virgilio.it>
00006 
00007 
00008 This program is free software; you can redistribute it and/or modify
00009 it under the terms of the GNU General Public License as published by
00010 the Free Software Foundation; either version 2 of the License, or
00011 (at your option) any later version.
00012 
00013 This program is distributed in the hope that it will be useful,
00014 but WITHOUT ANY WARRANTY; without even the implied warranty of
00015 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00016 GNU General Public License for more details.
00017 
00018 You should have received a copy of the GNU General Public License
00019 along with this program; if not, write to the Free Software
00020 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00021 ***********************************************************************/
00022 #include "../include/VarData.h"
00023 ;            
00025 int VarData::DensityProfile(int coord,int NBin,int NType,double *DensityAv){
00026   BfDefChain();
00027   int NPartition = 5;
00028   double *Density = (double *)calloc(NPartition*NPartition*NType*NBin,sizeof(double));
00029   double Pos[4];
00030   char cLine[STRSIZE];
00031   for(int b=0;b<Gen->NBlock;b++){
00032     if(!strncmp(Block[b].Name,"PEP",3) ) continue;
00033     if(!strcmp(Block[b].Name,"WATER") ) continue;
00034     for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
00035       int t = Pm[p].Typ;
00036       for(int d=0;d<3;d++)
00037    Pos[d] = Pm[p].Pos[d];
00038       int vx = (int)(NPartition*Pos[CLat1]/Gen->Edge[CLat2]);
00039       if(vx < 0 || vx >= NBin) continue;
00040       int vy = (int)(NPartition*Pos[CLat2]/Gen->Edge[CLat2]);
00041       if(vy < 0 || vy >= NBin) continue;
00042       int vz = (int)(NBin*Pos[CNorm]/Gen->Edge[CNorm]);
00043       if(vz < 0 || vz >= NBin) continue;
00044       if(CHAIN_IF_TYPE(Ch[Pm[p].CId].Type,CHAIN_ADDED))
00045    t = 2;
00046       Density[((vz*NType+t)*NPartition+vy)*NPartition+vx] += 1.;
00047     }
00048   }
00049   for(int vx=0;vx<NPartition;vx++){
00050     for(int vy=0;vy<NPartition;vy++){
00051       double Media = 0.;
00052       double Weight = 0.;
00053       for(int v=0;v<NBin;v++){
00054    Weight += Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
00055    Media += v*Density[((v*NType+0)*NPartition+vy)*NPartition+vx];
00056       }
00057       int vMedia = Weight > 0 ? (int)(Media/Weight) : NBin/2;
00058       int vOffset = vMedia - NBin/2;
00059       for(int t=0;t<NType;t++){
00060    for(int v=0;v<NBin;v++){
00061      int vOther = v + vOffset;
00062      if(vOther < 0)
00063        vOther = NBin - vOther;
00064      if(vOther >= NBin)
00065        vOther = vOther - NBin;
00066      DensityAv[v*NType+t] += Density[((vOther*NType+t)*NPartition+vy)*NPartition+vx];
00067    }
00068       }
00069     }
00070   }
00071   free(Density);
00072   return 0;
00073 }
00074 // Volume element in cylindrical symmetry
00075 void VarData::VolumeCircSlab(double *VolContr,int NBin){
00076   if(pEdge(3) <= 0.){printf("Lateral radius not defined\n");};
00077   double LenMin = .5*MIN(Gen->Edge[CLat1],Gen->Edge[CLat2]);
00078   double LenMax = .5*MAX(Gen->Edge[CLat1],Gen->Edge[CLat2]);
00079   double RapCorner = LenMin/LenMax*NBin;
00080   double RapMin = LenMin/Gen->Edge[3]*NBin;
00081   double Volume2 = 0.;
00082   for(int v=0;v<NBin;v++){
00083     double Rad2 = (v+1);
00084     double Volume1 = PI*QUAD(Rad2);
00085     if(Rad2 > RapMin){
00086       double Pro2 = sqrt(QUAD(v+1)-QUAD(RapMin));
00087       double Angle2 = atan(RapMin/Pro2);
00088       Volume1 = 2.*QUAD(Rad2)*Angle2 + 2.*Pro2*RapMin;
00089       if(Rad2 > RapCorner){
00090    double theta = acos(RapCorner/Rad2);
00091    double b = Rad2*sin(theta);
00092    //double b = sqrt( SQR(Rad2) - SQR(RapCorner));
00093    Volume1 -= 4.*(Rad2*Rad2*theta - b*RapCorner);
00094    //printf("%lf %lf %lf %lf %lf %lf %lf\n",Rad2,RapCorner,theta,b,Volume1-Volume2,Rad2*Rad2*theta ,b*RapCorner);
00095       }
00096     }
00097     double Volume = (Volume1 - Volume2)*Gen->Edge[CNorm]*SQR(Gen->Edge[3])/(double)CUB(NBin);
00098     VolContr[v] = Volume;
00099     Volume2 = Volume1;
00100   }
00101 }
00102 int VarData::Worm(int Partition,int NBin,double *Border,double *dPoint){
00103   int cLong = CLat2;
00104   int cLat = CLat1;
00105   double *Cm;
00106   RETTA *r1;
00107   r1 = (RETTA *)calloc(Partition,sizeof(RETTA));
00108   Cm = (double *)calloc(3*Partition,sizeof(double));
00109   double* Dis0 = (double *)calloc(Gen->NPart,sizeof(double));
00110   double* Dis1  = (double *)calloc(Gen->NPart,sizeof(double));
00111   for(int pa = 0,cc=0;pa<Partition;pa++){
00112     for(int c=0;c<Gen->NChain;c++){
00113       int v = (int)(Ch[c].Pos[cLong]*Partition/Gen->Edge[cLong]);
00114       if(v != pa) continue;
00115       *(Dis0 + cc) = Ch[c].Pos[cLat];
00116       *(Dis1 + cc) = Ch[c].Pos[cLong];
00117       //printf("%lf %lf %lf\n",Ch[c].Pos[0],Ch[c].Pos[1],Ch[c].Pos[2]);
00118       for(int d=0;d<3;d++)
00119    *(Cm + 3*v + d) += Ch[c].Pos[d];
00120       cc++;
00121     }
00122     r1[pa] = Mat->InterRett(Dis0,Dis1,cc);
00123     for(int d=0;d<3;d++)
00124       *(Cm + 3*pa + d) /= cc;
00125     cc = 0;
00126   }
00127   double pos0;
00128   double pos1;
00129   double qcm;
00130   double qpoint = 0.;
00131   double dist;
00132   int pp0 = 0;
00133   int pp1 = 0;
00134   //FILE *Ciccia = fopen("Ciccia.dat","w");
00135   for(int pa=0;pa<Partition;pa++){
00136     //printf("m %lf Cm %lf %lf %lf \n",r1[pa].m,*(Cm+3*pa+0),*(Cm+3*pa+1),*(Cm+3*pa+2));
00137     for(int p=0;p<Gen->NPart;p++){
00138       int v = (int) (Pm[p].Pos[cLong]*Partition/Gen->Edge[cLong]);
00139       if(v != pa) continue;
00140       qcm = *(Cm + 3*pa + cLat) + *(Cm + 3*pa + cLong)/r1[pa].m;
00141       qpoint = Pm[p].Pos[cLat] - r1[pa].m*Pm[p].Pos[cLong];
00142       pos0 = (qcm-qpoint)*r1[pa].m/(QUAD(r1[pa].m) + 1.);
00143       pos1 = pos0*r1[pa].m + qpoint;
00144       dist = QUAD((pos0-Cm[3*pa+cLong])) + QUAD((pos1-Cm[3*pa+cLat]));
00145       if(Pm[p].Typ == 0){
00146    Dis0[pp0] = pow(dist,.5);
00147    if(Pm[p].Pos[cLat] < *(Cm+3*pa+cLat))
00148      Dis0[pp0] *= -1.;
00149    //printf("Dis0 %lf\n",Dis0[pp0]);
00150    pp0++;
00151       }
00152       else if(Pm[p].Typ == 1){
00153    Dis1[pp1] = pow(dist,.5);
00154    if(Pm[p].Pos[cLat] < *(Cm+3*pa+cLat))
00155      Dis1[pp1] *= -1.;
00156    pp1++;
00157       }
00158       if(pa == 2){
00159    //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));
00160       }
00161     }
00162     //      fclose(Ciccia);
00163     for(int v=0;v<NBin;v++)
00164       dPoint[v*2] = 0.;
00165     for(int i=0;i<pp0;i++){
00166       int v = (int) ((Dis0[i]-Border[0])/(Border[1] - Border[0])*NBin);
00167       if(v>= NBin) continue;
00168       dPoint[v*2] += 1.;
00169     }
00170     for(int v=0;v<NBin;v++){
00171       dPoint[v*2] *= (NBin*Partition)/(Gen->Edge[0]*Gen->Edge[1]*Gen->Edge[2]);
00172     }
00173     for(int v=0;v<NBin;v++)
00174       dPoint[v*2] = 0.;
00175     for(int i=0;i<pp1;i++){
00176       int v = (int) ((Dis1[i]-Border[0])/(Border[1] - Border[0])*NBin);
00177       if(v>= NBin) continue;
00178       dPoint[v*2] += 1.;
00179     }
00180     for(int v=0;v<NBin;v++){
00181       dPoint[v*2+1] /= (double) (NBin*Partition)/(Gen->Edge[0]*Gen->Edge[1]*Gen->Edge[2]);
00182       //if(pa==0)printf("%d %lf %lf\n",v,dPoint[v],dPoint1[v]);
00183     }
00184     pp0=0;
00185     pp1=0;
00186   }
00187   for(int v=0;v<NBin;v++){
00188     dPoint[v*2] /= (double) NBin;
00189     dPoint[v*2+1] /= (double)NBin;
00190   }
00191   return 1;
00192 }
00193 void VarData::Stalk(int NBin,int NLevel,double **Plot,double Threshold){
00194   double **Norma = (double **)calloc(NLevel,sizeof(double));
00195   for(int l=0;l<NLevel;l++)
00196     Norma[l] = (double *)calloc(NBin*NBin,sizeof(double));
00197   for(int b=0;b<Gen->NBlock;b++){
00198     for(int p=Block[b].InitIdx;p<Block[b].EndIdx;p++){
00199       for(int v=0;v<NBin;v++){
00200    // if(Pm[p].Typ == 1 && Pm[p-1].Typ == 0)
00201      {
00202      int v = (int)(Pm[p].Pos[CLat1]/Gen->Edge[CLat1]*NBin);
00203      if( v < 0 || v >= NBin) continue;
00204      int vv = (int)(Pm[p].Pos[CLat2]/Gen->Edge[CLat2]*NBin);
00205      if( vv < 0 || vv >= NBin) continue;
00206      Plot[b][v*NBin + vv] += Pm[p].Pos[CNorm];
00207      Norma[b][v*NBin + vv] += 1.;
00208    }
00209       }
00210     }
00211   }
00212   for(int v=0;v<NBin;v++)
00213     for(int vv=0;vv<NBin;vv++)
00214       for(int l=0;l<NLevel;l++){
00215    if(Norma[l][v*NBin+vv] > Threshold)
00216      Plot[l][v*NBin+vv] /= Norma[l][v*NBin+vv];
00217    else
00218      Plot[l][v*NBin+vv] = 0.;
00219       }
00220   for(int l=0;l<NLevel;l++)
00221     free(Norma[l]);
00222   free(Norma);
00223 }