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

# 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     double r, theta, phi, f, s, w, w0;
19     double fp, sp, wp, tot;
20     double rup, ru, rlow;
21 gezelter 1295 char* shapeFileName;
22     SHAPE* shape;
23     FILE* my_file;
24 gezelter 1290
25     //parse the command line options
26     if (cmdline_parser (argc, argv, &args_info) != 0)
27     exit(1) ;
28    
29 gezelter 1295 if (args_info.shape_given){
30     shapeFileName = args_info.shape_arg;
31 gezelter 1290 }
32     else{
33     std::cerr << "Does not have shape file name" << endl;
34     exit(1);
35     }
36 gezelter 1295
37     shape = new SHAPE();
38     shape->readSHAPEfile(shapeFileName);
39    
40 gezelter 1290
41     // grid has a default value (default=51), so it is always given
42     npts = args_info.grid_arg;
43    
44 gezelter 1295
45 gezelter 1290
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     }