--- trunk/OOPSE/staticProps/GofRtheta.cpp	2003/10/21 19:33:19	811
+++ trunk/OOPSE/staticProps/GofRtheta.cpp	2003/12/16 20:49:11	879
@@ -12,18 +12,24 @@ GofRtheta::GofRtheta( char* key1, char* key2, int theN
 GofRtheta::GofRtheta( char* key1, char* key2, int theNatoms, int theNbins ):
   PairCorrType(  key1, key2, theNatoms, theNbins )
 {
-  int i;
+  int i, j;
   
   strcpy( corrType, "GofRtheta" );
   
-  currHist = new int[nBins];
-  currGofRtheta = new double[nBins];
-  avgGofRtheta  = new double[nBins];
-
+  currCosHist   = new int*[nBins];
+  currGofRtheta = new double*[nBins];
+  avgGofRtheta  = new double*[nBins];
+  
   for(i=0;i<nBins;i++){
-    currHist[i] = 0;
-    currGofRtheta[i] = 0.0;
-    avgGofRtheta[i]  = 0.0;
+    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;
@@ -31,25 +37,59 @@ GofRtheta::~GofRtheta(){
 
 GofRtheta::~GofRtheta(){
 
-  delete[] currHist;
+  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 bin;
-
+			   double uHatI[3], double uHatJ[3] ){
+  int binR, binTheta, i;
+  double dot, cosTheta;
+  double halfNbins;
+  
   if( correlateMe ){
-    
-    bin = (int)( dist / delR );
-    if( bin < nBins )currHist[bin] += 2;
-    
+
+    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;
+  int i,j;
   double rLower, rUpper, volSlice;
   double nIdeal;
 
@@ -63,17 +103,20 @@ void GofRtheta::accumulateFrame( void ){
     volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
     nIdeal = volSlice * pairConstant;
 
-    currGofRtheta[i] = currHist[i] / nIdeal;
-    currHist[i] = 0;
+    for(j=0;j<nBins;j++){
+      
+      currGofRtheta[i][j] = (currCosHist[i][j] / nIdeal);
+      currCosHist[i][j]   = 0;
 
-    avgGofRtheta[i] += currGofRtheta[i];    
+      avgGofRtheta[i][j] += currGofRtheta[i][j];    
+    }
   }
 }
 
 
 void GofRtheta::writeCorr( char* outPrefix ){
 
-  double rValue, corrValue;
+  double rValue, cosValue, deltaCos, corrValue;
   double rLower, rUpper, volSlice;
   double nIdeal;
   int i;
@@ -95,14 +138,22 @@ void GofRtheta::writeCorr( char* outPrefix ){
     simError();
   }
 
-  outStream << "#rValue\tcorrValue\n";
+  outStream << "#rValue; cosValue; corrValue\n";
 
-  for(i=0;i<nBins;i++){
 
+  deltaCos = 2.0 / nBins;
+  for(i=0;i<nBins;i++){
+    
     rValue = ( i + 0.5 ) * delR;
-    corrValue = avgGofRtheta[i] / nFrames;
-
-    outStream << rValue << "\t" << corrValue << "\n";
+    
+    for(j=0;j<nBins;j++){
+     
+      cosValue = -1 + (j * deltaCos);
+	
+      corrValue = avgGofRtheta[i][j] / nFrames;
+      
+      outStream << rValue << "; " << cosValue << "; " << corrValue << "\n";
+    }
   }
 
   outStream.close();