OOPSE/staticProps/GofRtheta.cpp

#include <iostream>
#include <fstream>

#include <cstring>
#include <cmath>

#include "simError.h"
#include "PairCorrType.hpp"

using namespace std;

GofRtheta::GofRtheta( char* key1, char* key2, int theNatoms, int theNbins ):
  PairCorrType(  key1, key2, theNatoms, theNbins )
{
  int i, j;
  
  strcpy( corrType, "GofRtheta" );
  
  currCosHist   = new int*[nBins];
  currGofRtheta = new double*[nBins];
  avgGofRtheta  = new double*[nBins];
  
  for(i=0;i<nBins;i++){
    currCosHist[i]   = new int[nBins];
    currGofRtheta[i] = new double[nBins];
    avgGofRtheta[i]  = new double[nBins];
    
    for(j=0;j<nBins;j++){
      currCosHist[i][j]   = 0;
      currGofRtheta[i][j] = 0.0;
      avgGofRtheta[i][j]  = 0.0;
    }
  }
    
  nFrames = 0;
}

GofRtheta::~GofRtheta(){

  int i;

  for(i=0;i<nBins;i++){
    delete[] currCosHist[i];
    delete[] currGofRtheta[i];
    delete[] avgGofRtheta[i];
  }

  delete[] currCosHist;
  delete[] currGofRtheta;
  delete[] avgGofRtheta;
}

void GofRtheta::correlate( double Rij[3], double dist, 
                           double uHatI[3], double uHatJ[3] ){
  int binR, binTheta, i;
  double dot, cosTheta;
  double halfNbins;
  
  if( correlateMe ){

    if( uHatI != NULL ){
      
      dot= 0.0;
      for(i=0;i<3;i++)
        dot += Rij[i] * uHatI[i];
      
      cosTheta = dot / dist;
      
      binR = (int)( dist / delR );
      if( binR < nBins ){
        
        halfNbins = (nBins-1) * 0.5;
        binTheta = (int)( (cosTheta * halfNbins) + halfNbins );
        
        currCosHist[binR][binTheta] += 1;
      }
    } 
    else{
      
      sprintf( painCave.errMsg,
               "GofRtheta error: atom \"%s\" does not have a "
               "directional component\n",
               atom1
               );
      painCave.isFatal = 1;
      simError();
    }    
  }
}

void GofRtheta::accumulateFrame( void ){
  int i,j;
  double rLower, rUpper, volSlice;
  double nIdeal;

  nFrames++;
  
  for(i=0;i<nBins;i++){

    rLower = i * delR;
    rUpper = rLower + delR;

    volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
    nIdeal = volSlice * pairConstant;

    for(j=0;j<nBins;j++){
      
      currGofRtheta[i][j] = (currCosHist[i][j] / nIdeal);
      currCosHist[i][j]   = 0;

      avgGofRtheta[i][j] += currGofRtheta[i][j];    
    }
  }
}


void GofRtheta::writeCorr( char* outPrefix ){

  double rValue, cosValue, deltaCos, corrValue;
  double rLower, rUpper, volSlice;
  double nIdeal;
  int i;
  char outName[200];
  
  sprintf( outName,
           "%s-%s-%s.GofRtheta",
           outPrefix,
           atom1,
           atom2 );
  ofstream outStream( outName );

  if( !outStream ){

    sprintf( painCave.errMsg,
             "Error opening \"%s\" for output.\n",
             outName );
    painCave.isFatal = 1;
    simError();
  }

  outStream << "#rValue; cosValue; corrValue\n";


  deltaCos = 2.0 / nBins;
  for(i=0;i<nBins;i++){
    
    rValue = ( i + 0.5 ) * delR;
    
    for(j=0;j<nBins;j++){
     
      cosValue = -1 + (j * deltaCos);
        
      corrValue = avgGofRtheta[i][j] / nFrames;
      
      outStream << rValue << "; " << cosValue << "; " << corrValue << "\n";
    }
  }

  outStream.close();
}
Revision:	879
Committed:	Tue Dec 16 20:49:11 2003 UTC (21 years, 10 months ago) by mmeineke
File size:	3093 byte(s)
Log Message:	finished gofRtheta and added gofRomega. both need to be debugged and tested.
#	User	Rev	Content
1	mmeineke	811	#include <iostream>
2			#include <fstream>
3
4			#include <cstring>
5			#include <cmath>
6
7			#include "simError.h"
8			#include "PairCorrType.hpp"
9
10			using namespace std;
11
12			GofRtheta::GofRtheta( char* key1, char* key2, int theNatoms, int theNbins ):
13			PairCorrType( key1, key2, theNatoms, theNbins )
14			{
15	mmeineke	873	int i, j;
16	mmeineke	811
17			strcpy( corrType, "GofRtheta" );
18
19	mmeineke	873	currCosHist = new int*[nBins];
20			currGofRtheta = new double*[nBins];
21			avgGofRtheta = new double*[nBins];
22
23	mmeineke	811	for(i=0;i<nBins;i++){
24	mmeineke	873	currCosHist[i] = new int[nBins];
25			currGofRtheta[i] = new double[nBins];
26			avgGofRtheta[i] = new double[nBins];
27
28			for(j=0;j<nBins;j++){
29			currCosHist[i][j] = 0;
30			currGofRtheta[i][j] = 0.0;
31			avgGofRtheta[i][j] = 0.0;
32			}
33	mmeineke	811	}
34
35			nFrames = 0;
36			}
37
38			GofRtheta::~GofRtheta(){
39
40	mmeineke	873	int i;
41
42			for(i=0;i<nBins;i++){
43			delete[] currCosHist[i];
44			delete[] currGofRtheta[i];
45			delete[] avgGofRtheta[i];
46			}
47
48	mmeineke	879	delete[] currCosHist;
49	mmeineke	811	delete[] currGofRtheta;
50			delete[] avgGofRtheta;
51			}
52
53			void GofRtheta::correlate( double Rij[3], double dist,
54	mmeineke	873	double uHatI[3], double uHatJ[3] ){
55	mmeineke	879	int binR, binTheta, i;
56	mmeineke	873	double dot, cosTheta;
57			double halfNbins;
58
59			if( correlateMe ){
60	mmeineke	811
61	mmeineke	873	if( uHatI != NULL ){
62
63			dot= 0.0;
64			for(i=0;i<3;i++)
65	mmeineke	879	dot += Rij[i] * uHatI[i];
66	mmeineke	873
67			cosTheta = dot / dist;
68
69	mmeineke	879	binR = (int)( dist / delR );
70			if( binR < nBins ){
71	mmeineke	873
72			halfNbins = (nBins-1) * 0.5;
73	mmeineke	879	binTheta = (int)( (cosTheta * halfNbins) + halfNbins );
74	mmeineke	873
75	mmeineke	879	currCosHist[binR][binTheta] += 1;
76	mmeineke	873	}
77	mmeineke	879	}
78			else{
79
80			sprintf( painCave.errMsg,
81			"GofRtheta error: atom \"%s\" does not have a "
82			"directional component\n",
83			atom1
84			);
85			painCave.isFatal = 1;
86			simError();
87	mmeineke	873	}
88	mmeineke	811	}
89			}
90
91			void GofRtheta::accumulateFrame( void ){
92	mmeineke	879	int i,j;
93	mmeineke	811	double rLower, rUpper, volSlice;
94			double nIdeal;
95
96			nFrames++;
97
98			for(i=0;i<nBins;i++){
99
100			rLower = i * delR;
101			rUpper = rLower + delR;
102
103			volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
104			nIdeal = volSlice * pairConstant;
105
106	mmeineke	873	for(j=0;j<nBins;j++){
107
108	mmeineke	879	currGofRtheta[i][j] = (currCosHist[i][j] / nIdeal);
109	mmeineke	873	currCosHist[i][j] = 0;
110	mmeineke	811
111	mmeineke	873	avgGofRtheta[i][j] += currGofRtheta[i][j];
112			}
113	mmeineke	811	}
114			}
115
116
117			void GofRtheta::writeCorr( char* outPrefix ){
118
119	mmeineke	873	double rValue, cosValue, deltaCos, corrValue;
120	mmeineke	811	double rLower, rUpper, volSlice;
121			double nIdeal;
122			int i;
123			char outName[200];
124
125			sprintf( outName,
126			"%s-%s-%s.GofRtheta",
127			outPrefix,
128			atom1,
129			atom2 );
130			ofstream outStream( outName );
131
132			if( !outStream ){
133
134			sprintf( painCave.errMsg,
135			"Error opening \"%s\" for output.\n",
136			outName );
137			painCave.isFatal = 1;
138			simError();
139			}
140
141	mmeineke	873	outStream << "#rValue; cosValue; corrValue\n";
142	mmeineke	811
143	mmeineke	873
144			deltaCos = 2.0 / nBins;
145	mmeineke	811	for(i=0;i<nBins;i++){
146	mmeineke	873
147	mmeineke	811	rValue = ( i + 0.5 ) * delR;
148	mmeineke	873
149			for(j=0;j<nBins;j++){
150
151			cosValue = -1 + (j * deltaCos);
152
153	mmeineke	879	corrValue = avgGofRtheta[i][j] / nFrames;
154	mmeineke	873
155			outStream << rValue << "; " << cosValue << "; " << corrValue << "\n";
156			}
157	mmeineke	811	}
158
159			outStream.close();
160			}