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

  shape = new SHAPE();
  shape->readSHAPEfile(shapeFileName);

  if (args_info.inputs_num > 1) {
      outputFileName = args_info.inputs[1];
  } else {
    std::cerr << "No output file name specified." << endl;
    outputFileName = strdup(shapeFileName);
    strcat(outputFileName, ".viz");
    std::cerr << "Using: " << outputFileName << endl;
  }

  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;

  xmin = -range;
  xmax = range;

  ymin = -range;
  ymax = range;

  zmin = -range;
  zmax = range;

  //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:	1361
Committed:	Tue Jul 20 20:22:43 2004 UTC (19 years, 11 months ago) by gezelter
File size:	2728 byte(s)
Log Message:	using gengetopts --unamed-opts="SHAPEFILE [OUTPUTFILE]"
#	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	chrisfen	1315	double range;
25	gezelter	1299
26	gezelter	1295	char* shapeFileName;
27	gezelter	1299	char* outputFileName;
28	gezelter	1295	SHAPE* shape;
29	gezelter	1299	FILE* outputFile;
30	gezelter	1290
31			//parse the command line options
32			if (cmdline_parser (argc, argv, &args_info) != 0)
33			exit(1) ;
34
35	gezelter	1361
36			if (args_info.inputs_num > 0) {
37			shapeFileName = args_info.inputs[0];
38			} else {
39	gezelter	1290	std::cerr << "Does not have shape file name" << endl;
40			exit(1);
41			}
42	gezelter	1295
43			shape = new SHAPE();
44			shape->readSHAPEfile(shapeFileName);
45
46	gezelter	1361	if (args_info.inputs_num > 1) {
47			outputFileName = args_info.inputs[1];
48			} else {
49			std::cerr << "No output file name specified." << endl;
50			outputFileName = strdup(shapeFileName);
51			strcat(outputFileName, ".viz");
52			std::cerr << "Using: " << outputFileName << endl;
53	gezelter	1299	}
54	gezelter	1290
55	gezelter	1299	outputFile = fopen(outputFileName, "w");
56	gezelter	1290
57	gezelter	1299	if (outputFile == NULL) {
58			(void) fprintf(stderr, "Unable to open outputFile %s!\n", outputFileName);
59	gezelter	1290	exit(8);
60			}
61	gezelter	1299
62			// grid has a default value (default=51), so it is always given
63			npts = args_info.grid_arg;
64	chrisfen	1315	range = args_info.range_arg;
65	gezelter	1290
66	chrisfen	1315	xmin = -range;
67			xmax = range;
68	gezelter	1290
69	chrisfen	1315	ymin = -range;
70			ymax = range;
71	gezelter	1290
72	chrisfen	1315	zmin = -range;
73			zmax = range;
74	gezelter	1290
75	gezelter	1300	//sigmaProbe = 2.28;
76			//sProbe = 2.28;
77			//epsProbe = 0.020269601874;
78			sigmaProbe = 0.0;
79			sProbe = 0.0;
80			epsProbe = 1.0;
81
82	gezelter	1290	for (i = 0; i < npts; i++) {
83	chrisfen	1311	x = xmin + (xmax-xmin) * (double)i/(double)(npts-1);
84	gezelter	1290
85			for (j = 0; j < npts; j++) {
86	chrisfen	1311	y = ymin + (ymax-ymin) * (double)j/(double)(npts-1);
87	gezelter	1290
88			for (k = 0; k < npts; k++) {
89	chrisfen	1311	z = zmin + (zmax-zmin) * (double)k/(double)(npts-1);
90	gezelter	1290
91			r = sqrt(xx + yy + z*z);
92	gezelter	1300	costheta = z/r;
93	gezelter	1307	phi = atan2(y,x);
94	gezelter	1290
95	gezelter	1299	sigmaShape = shape->getSigmaAt(costheta, phi);
96			sShape = shape->getSAt(costheta, phi);
97			epsShape = shape->getEpsAt(costheta, phi);
98	gezelter	1290
99	gezelter	1299	// Lorentz-Berthellot combining rules:
100	gezelter	1290
101	gezelter	1299	sigma = (sigmaShape + sigmaProbe) / 2.0;
102			s = (sShape + sProbe) / 2.0;
103			eps = sqrt(epsShape * epsProbe);
104	gezelter	1290
105	gezelter	1299	rTerm = s / (r - sigma + s);
106
107			r6 = pow(rTerm, 6);
108			r12 = r6*r6;
109	gezelter	1290
110	gezelter	1299	energy = 4.0 * eps * (r12 - r6);
111
112			write_size = fwrite(&energy,1,8,outputFile);
113
114	gezelter	1290	}
115			}
116			}
117	gezelter	1299	fclose(outputFile);
118	gezelter	1290	return(0);
119			}