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[] currHist;
  delete[] currGofRtheta;
  delete[] avgGofRtheta;
}

void GofRtheta::correlate( double Rij[3], double dist, 
                           double uHatI[3], double uHatJ[3] ){
  int bin, 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;
      
      bin = (int)( dist / delR );
      if( bin < nBins ){
        
        halfNbins = (nBins-1) * 0.5;
        binTheta = (int)( (dot * halfNbins) + halfNbins );
        
        currCosHist[bin][binTheta] += 1;
      }
    }    
  }
}

void GofRtheta::accumulateFrame( void ){
  int i;
  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) * nBins;
      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] / nFrames;
      
      outStream << rValue << "; " << cosValue << "; " << corrValue << "\n";
    }
  }

  outStream.close();
}
Revision:	873
Committed:	Fri Nov 21 20:07:14 2003 UTC (20 years, 7 months ago) by mmeineke
File size:	2866 byte(s)
Log Message:	begun work on add ing in the GofR,CosTheta
#	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	811	delete[] currHist;
49			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			int bin, binTheta, i;
56			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			dot += rij[i] * uHatI[i];
66
67			cosTheta = dot / dist;
68
69			bin = (int)( dist / delR );
70			if( bin < nBins ){
71
72			halfNbins = (nBins-1) * 0.5;
73			binTheta = (int)( (dot * halfNbins) + halfNbins );
74
75			currCosHist[bin][binTheta] += 1;
76			}
77			}
78	mmeineke	811	}
79			}
80
81			void GofRtheta::accumulateFrame( void ){
82			int i;
83			double rLower, rUpper, volSlice;
84			double nIdeal;
85
86			nFrames++;
87
88			for(i=0;i<nBins;i++){
89
90			rLower = i * delR;
91			rUpper = rLower + delR;
92
93			volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
94			nIdeal = volSlice * pairConstant;
95
96	mmeineke	873	for(j=0;j<nBins;j++){
97
98			currGofRtheta[i][j] = (currCosHist[i][j] / nIdeal) * nBins;
99			currCosHist[i][j] = 0;
100	mmeineke	811
101	mmeineke	873	avgGofRtheta[i][j] += currGofRtheta[i][j];
102			}
103	mmeineke	811	}
104			}
105
106
107			void GofRtheta::writeCorr( char* outPrefix ){
108
109	mmeineke	873	double rValue, cosValue, deltaCos, corrValue;
110	mmeineke	811	double rLower, rUpper, volSlice;
111			double nIdeal;
112			int i;
113			char outName[200];
114
115			sprintf( outName,
116			"%s-%s-%s.GofRtheta",
117			outPrefix,
118			atom1,
119			atom2 );
120			ofstream outStream( outName );
121
122			if( !outStream ){
123
124			sprintf( painCave.errMsg,
125			"Error opening \"%s\" for output.\n",
126			outName );
127			painCave.isFatal = 1;
128			simError();
129			}
130
131	mmeineke	873	outStream << "#rValue; cosValue; corrValue\n";
132	mmeineke	811
133	mmeineke	873
134			deltaCos = 2.0 / nBins;
135	mmeineke	811	for(i=0;i<nBins;i++){
136	mmeineke	873
137	mmeineke	811	rValue = ( i + 0.5 ) * delR;
138	mmeineke	873
139			for(j=0;j<nBins;j++){
140
141			cosValue = -1 + (j * deltaCos);
142
143			corrValue = avgGofRtheta[i] / nFrames;
144
145			outStream << rValue << "; " << cosValue << "; " << corrValue << "\n";
146			}
147	mmeineke	811	}
148
149			outStream.close();
150			}