--- trunk/SHAPES/visualizer.cpp	2004/06/23 20:53:17	1290
+++ trunk/SHAPES/visualizer.cpp	2004/06/28 23:05:25	1315
@@ -2,10 +2,10 @@
 #include <fstream>
 #include <string>
 #include <vector>
+#include <cmath>
 #include "visualizerCmd.h"
+#include "SHAPE.hpp"
 
-
-int count_tokens(char *line, char *delimiters);
 using namespace std;
 
 int main(int argc, char* argv[]){
@@ -15,91 +15,105 @@ int main(int argc, char* argv[]){
   double xmin, xmax, x;
   double ymin, ymax, y;
   double zmin, zmax, z;
-  double r, theta, phi, f, s, w, w0;
-  double fp, sp, wp, tot;
-  double rup, ru, rlow;
-  FILE *my_file;
+  double r, theta, phi, costheta;
+  double sigmaShape, sShape, epsShape;
+  double sigmaProbe, sProbe, epsProbe;
+  double sigma, s, eps;
+  double rTerm, r6, r12;
+  double energy;
+  double range;
 
+  char* shapeFileName;
+  char* outputFileName;
+  SHAPE* shape;
+  FILE* outputFile;
+
   //parse the command line options
   if (cmdline_parser (argc, argv, &args_info) != 0)
     exit(1) ;
 
-  if (args_info.input_given){
-    shapeFileName = args_info.input_arg;
+  if (args_info.shape_given){
+    shapeFileName = args_info.shape_arg;
   }
   else{
     std::cerr << "Does not have shape file name" << endl;
     exit(1);
   }
-  
-  // grid has a default value (default=51), so it is always given
-  npts = args_info.grid_arg;
 
-  printf ("opening %s\n", shapeFileName);
-  shapeFile = fopen( shapeFileName, "r" );
+  shape = new SHAPE();
+  shape->readSHAPEfile(shapeFileName);
 
 
+  if (args_info.output_given){
+    outputFileName = args_info.output_arg;
+  }
+  else{
+    std::cerr << "Does not have output file name" << endl;
+    exit(1);
+  }
 
-  my_file = fopen("junk.dat", "w");
-  if (my_file == NULL) {
-    (void) fprintf(stderr, "No File\n");
+  outputFile = fopen(outputFileName, "w");
+
+  if (outputFile == NULL) {
+    (void) fprintf(stderr, "Unable to open outputFile %s!\n", outputFileName);
     exit(8);
   }
+    
+  // grid has a default value (default=51), so it is always given
+  npts = args_info.grid_arg;
+  range = args_info.range_arg;
 
-  npts = 51;
+  xmin = -range;
+  xmax = range;
 
-  xmin = -5.0;
-  xmax = 5.0;
+  ymin = -range;
+  ymax = range;
 
-  ymin = -5.0;
-  ymax = 5.0;
+  zmin = -range;
+  zmax = range;
 
-  zmin = -5.0;
-  zmax = 5.0;
+  //sigmaProbe = 2.28;
+  //sProbe = 2.28;
+  //epsProbe = 0.020269601874;    
+  sigmaProbe = 0.0;
+  sProbe = 0.0;
+  epsProbe = 1.0;    
 
-  rup = 4.0;
-  ru = 3.35;
-  rlow = 2.75;
-  w0 = 0.07715;
-
   for (i = 0; i < npts; i++) {
-    x = xmin + (xmax-xmin) * (double)i/(double)npts;
+    x = xmin + (xmax-xmin) * (double)i/(double)(npts-1);
 
     for (j = 0; j < npts; j++) {
-      y = ymin + (ymax-ymin) * (double)j/(double)npts;
+      y = ymin + (ymax-ymin) * (double)j/(double)(npts-1);
 
       for (k = 0; k < npts; k++) {
-        z = zmin + (zmax-zmin) * (double)k/(double)npts;
+        z = zmin + (zmax-zmin) * (double)k/(double)(npts-1);
 
         r = sqrt(x*x + y*y + z*z);
-        theta = acos(z/r);
-        phi = atan(y/x);
-          
-        fp = pow((cos(theta) - 0.6),2)  * pow((cos(theta) + 0.8),2);
+        costheta = z/r;
+        phi = atan2(y,x);
 
-        sp = (rup + 2.0*r - 3.0*rlow) * pow((rup - r),2) / pow((rup -rlow),3);
+        sigmaShape = shape->getSigmaAt(costheta, phi);
+        sShape = shape->getSAt(costheta, phi);
+        epsShape = shape->getEpsAt(costheta, phi);
 
-        if (r > rup) sp = 0.0;
-        if (r < rlow) sp = 1.0;
+        // Lorentz-Berthellot combining rules:
 
-        wp = sp*(fp - w0);
+        sigma = (sigmaShape + sigmaProbe) / 2.0;
+        s = (sShape + sProbe) / 2.0;
+        eps = sqrt(epsShape * epsProbe);
 
-        f = sin(theta)*sin(2.0*theta)*cos(2.0*phi);
-        
-        s = (ru + 2.0*r - 3.0*rlow) * pow((ru - r),2) / pow((ru -rlow),3);
+        rTerm = s / (r - sigma + s);
 
-        if (r > ru) s = 0.0;
-        if (r < rlow) s = 1.0;
-
-        w = f*s;
-
-        tot = w + wp;
-
-        write_size = fwrite(&tot,1,8,my_file);
-
+        r6 = pow(rTerm, 6);
+        r12 = r6*r6;
+        
+        energy = 4.0 * eps * (r12 - r6);
+        
+        write_size = fwrite(&energy,1,8,outputFile);
+        
       }
     }
   }
-  fclose(my_file);
+  fclose(outputFile);
   return(0);
 }