ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/SHAPES/visualizer.cpp
Revision: 1295
Committed: Thu Jun 24 15:31:52 2004 UTC (20 years ago) by gezelter
File size: 2136 byte(s)
Log Message:
renamed forcer, modified factorial, fixed makefile

File Contents

# Content
1 #include <iostream>
2 #include <fstream>
3 #include <string>
4 #include <vector>
5 #include <cmath>
6 #include "visualizerCmd.h"
7 #include "SHAPE.hpp"
8
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 double r, theta, phi, f, s, w, w0;
19 double fp, sp, wp, tot;
20 double rup, ru, rlow;
21 char* shapeFileName;
22 SHAPE* shape;
23 FILE* my_file;
24
25 //parse the command line options
26 if (cmdline_parser (argc, argv, &args_info) != 0)
27 exit(1) ;
28
29 if (args_info.shape_given){
30 shapeFileName = args_info.shape_arg;
31 }
32 else{
33 std::cerr << "Does not have shape file name" << endl;
34 exit(1);
35 }
36
37 shape = new SHAPE();
38 shape->readSHAPEfile(shapeFileName);
39
40
41 // grid has a default value (default=51), so it is always given
42 npts = args_info.grid_arg;
43
44
45
46
47 my_file = fopen("junk.dat", "w");
48 if (my_file == NULL) {
49 (void) fprintf(stderr, "No File\n");
50 exit(8);
51 }
52
53 npts = 51;
54
55 xmin = -5.0;
56 xmax = 5.0;
57
58 ymin = -5.0;
59 ymax = 5.0;
60
61 zmin = -5.0;
62 zmax = 5.0;
63
64 rup = 4.0;
65 ru = 3.35;
66 rlow = 2.75;
67 w0 = 0.07715;
68
69 for (i = 0; i < npts; i++) {
70 x = xmin + (xmax-xmin) * (double)i/(double)npts;
71
72 for (j = 0; j < npts; j++) {
73 y = ymin + (ymax-ymin) * (double)j/(double)npts;
74
75 for (k = 0; k < npts; k++) {
76 z = zmin + (zmax-zmin) * (double)k/(double)npts;
77
78 r = sqrt(x*x + y*y + z*z);
79 theta = acos(z/r);
80 phi = atan(y/x);
81
82 fp = pow((cos(theta) - 0.6),2) * pow((cos(theta) + 0.8),2);
83
84 sp = (rup + 2.0*r - 3.0*rlow) * pow((rup - r),2) / pow((rup -rlow),3);
85
86 if (r > rup) sp = 0.0;
87 if (r < rlow) sp = 1.0;
88
89 wp = sp*(fp - w0);
90
91 f = sin(theta)*sin(2.0*theta)*cos(2.0*phi);
92
93 s = (ru + 2.0*r - 3.0*rlow) * pow((ru - r),2) / pow((ru -rlow),3);
94
95 if (r > ru) s = 0.0;
96 if (r < rlow) s = 1.0;
97
98 w = f*s;
99
100 tot = w + wp;
101
102 write_size = fwrite(&tot,1,8,my_file);
103
104 }
105 }
106 }
107 fclose(my_file);
108 return(0);
109 }