trunk/SHAPES/visualizer.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cmath>
#include "visualizerCmd.h"
#include "SHAPE.hpp"

using namespace std;

int main(int argc, char* argv[]){

  gengetopt_args_info args_info;
  int npts, i, j, k, write_size;
  double xmin, xmax, x;
  double ymin, ymax, y;
  double zmin, zmax, z;
  double r, theta, phi, costheta;
  double sigmaShape, sShape, epsShape;
  double sigmaProbe, sProbe, epsProbe;
  double sigma, s, eps;
  double rTerm, r6, r12;
  double energy;

  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.shape_given){
    shapeFileName = args_info.shape_arg;
  }
  else{
    std::cerr << "Does not have shape file name" << endl;
    exit(1);
  }

  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);
  }

  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;

  xmin = -10.0;
  xmax = 10.0;

  ymin = -10.0;
  ymax = 10.0;

  zmin = -10.0;
  zmax = 10.0;

  //sigmaProbe = 2.28;
  //sProbe = 2.28;
  //epsProbe = 0.020269601874;    
  sigmaProbe = 0.0;
  sProbe = 0.0;
  epsProbe = 1.0;    

  for (i = 0; i < npts; i++) {
    x = xmin + (xmax-xmin) * (double)i/(double)(npts-1);

    for (j = 0; j < npts; j++) {
      y = ymin + (ymax-ymin) * (double)j/(double)(npts-1);

      for (k = 0; k < npts; k++) {
        z = zmin + (zmax-zmin) * (double)k/(double)(npts-1);

        r = sqrt(x*x + y*y + z*z);
        costheta = z/r;
        phi = atan2(y,x);

        sigmaShape = shape->getSigmaAt(costheta, phi);
        sShape = shape->getSAt(costheta, phi);
        epsShape = shape->getEpsAt(costheta, phi);

        // Lorentz-Berthellot combining rules:

        sigma = (sigmaShape + sigmaProbe) / 2.0;
        s = (sShape + sProbe) / 2.0;
        eps = sqrt(epsShape * epsProbe);

        rTerm = s / (r - sigma + s);

        r6 = pow(rTerm, 6);
        r12 = r6*r6;
        
        energy = 4.0 * eps * (r12 - r6);
        
        write_size = fwrite(&energy,1,8,outputFile);
        
      }
    }
  }
  fclose(outputFile);
  return(0);
}
Revision:	1311
Committed:	Fri Jun 25 22:16:39 2004 UTC (20 years ago) by chrisfen
File size:	2546 byte(s)
Log Message:	Recentering of the visualized potentials, adding tolerance option
#	User	Rev	Content
1	gezelter	1290	#include <iostream>
2			#include <fstream>
3			#include <string>
4			#include <vector>
5	gezelter	1295	#include <cmath>
6	gezelter	1290	#include "visualizerCmd.h"
7	gezelter	1295	#include "SHAPE.hpp"
8	gezelter	1290
9			using namespace std;
10
11			int main(int argc, char* argv[]){
12
13			gengetopt_args_info args_info;
14			int npts, i, j, k, write_size;
15			double xmin, xmax, x;
16			double ymin, ymax, y;
17			double zmin, zmax, z;
18	gezelter	1299	double r, theta, phi, costheta;
19			double sigmaShape, sShape, epsShape;
20			double sigmaProbe, sProbe, epsProbe;
21			double sigma, s, eps;
22			double rTerm, r6, r12;
23			double energy;
24
25	gezelter	1295	char* shapeFileName;
26	gezelter	1299	char* outputFileName;
27	gezelter	1295	SHAPE* shape;
28	gezelter	1299	FILE* outputFile;
29	gezelter	1290
30			//parse the command line options
31			if (cmdline_parser (argc, argv, &args_info) != 0)
32			exit(1) ;
33
34	gezelter	1295	if (args_info.shape_given){
35			shapeFileName = args_info.shape_arg;
36	gezelter	1290	}
37			else{
38			std::cerr << "Does not have shape file name" << endl;
39			exit(1);
40			}
41	gezelter	1295
42			shape = new SHAPE();
43			shape->readSHAPEfile(shapeFileName);
44
45	gezelter	1290
46	gezelter	1299	if (args_info.output_given){
47			outputFileName = args_info.output_arg;
48			}
49			else{
50			std::cerr << "Does not have output file name" << endl;
51			exit(1);
52			}
53	gezelter	1290
54	gezelter	1299	outputFile = fopen(outputFileName, "w");
55	gezelter	1290
56	gezelter	1299	if (outputFile == NULL) {
57			(void) fprintf(stderr, "Unable to open outputFile %s!\n", outputFileName);
58	gezelter	1290	exit(8);
59			}
60	gezelter	1299
61			// grid has a default value (default=51), so it is always given
62			npts = args_info.grid_arg;
63	gezelter	1290
64	gezelter	1300	xmin = -10.0;
65			xmax = 10.0;
66	gezelter	1290
67	gezelter	1300	ymin = -10.0;
68			ymax = 10.0;
69	gezelter	1290
70	gezelter	1300	zmin = -10.0;
71			zmax = 10.0;
72	gezelter	1290
73	gezelter	1300	//sigmaProbe = 2.28;
74			//sProbe = 2.28;
75			//epsProbe = 0.020269601874;
76			sigmaProbe = 0.0;
77			sProbe = 0.0;
78			epsProbe = 1.0;
79
80	gezelter	1290	for (i = 0; i < npts; i++) {
81	chrisfen	1311	x = xmin + (xmax-xmin) * (double)i/(double)(npts-1);
82	gezelter	1290
83			for (j = 0; j < npts; j++) {
84	chrisfen	1311	y = ymin + (ymax-ymin) * (double)j/(double)(npts-1);
85	gezelter	1290
86			for (k = 0; k < npts; k++) {
87	chrisfen	1311	z = zmin + (zmax-zmin) * (double)k/(double)(npts-1);
88	gezelter	1290
89			r = sqrt(xx + yy + z*z);
90	gezelter	1300	costheta = z/r;
91	gezelter	1307	phi = atan2(y,x);
92	gezelter	1290
93	gezelter	1299	sigmaShape = shape->getSigmaAt(costheta, phi);
94			sShape = shape->getSAt(costheta, phi);
95			epsShape = shape->getEpsAt(costheta, phi);
96	gezelter	1290
97	gezelter	1299	// Lorentz-Berthellot combining rules:
98	gezelter	1290
99	gezelter	1299	sigma = (sigmaShape + sigmaProbe) / 2.0;
100			s = (sShape + sProbe) / 2.0;
101			eps = sqrt(epsShape * epsProbe);
102	gezelter	1290
103	gezelter	1299	rTerm = s / (r - sigma + s);
104
105			r6 = pow(rTerm, 6);
106			r12 = r6*r6;
107	gezelter	1290
108	gezelter	1299	energy = 4.0 * eps * (r12 - r6);
109
110			write_size = fwrite(&energy,1,8,outputFile);
111
112	gezelter	1290	}
113			}
114			}
115	gezelter	1299	fclose(outputFile);
116	gezelter	1290	return(0);
117			}