# | Line 2 | Line 2 | |
---|---|---|
2 | #include <fstream> | |
3 | #include <string> | |
4 | #include <vector> | |
5 | + | #include <cmath> |
6 | #include "visualizerCmd.h" | |
7 | + | #include "SHAPE.hpp" |
8 | ||
7 | – | |
8 | – | int count_tokens(char *line, char *delimiters); |
9 | using namespace std; | |
10 | ||
11 | int main(int argc, char* argv[]){ | |
# | Line 15 | Line 15 | int main(int argc, char* argv[]){ | |
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 | < | FILE *my_file; |
18 | > | 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 | + | char* shapeFileName; |
26 | + | char* outputFileName; |
27 | + | SHAPE* shape; |
28 | + | FILE* outputFile; |
29 | + | |
30 | //parse the command line options | |
31 | if (cmdline_parser (argc, argv, &args_info) != 0) | |
32 | exit(1) ; | |
33 | ||
34 | < | if (args_info.input_given){ |
35 | < | shapeFileName = args_info.input_arg; |
34 | > | if (args_info.shape_given){ |
35 | > | shapeFileName = args_info.shape_arg; |
36 | } | |
37 | else{ | |
38 | std::cerr << "Does not have shape file name" << endl; | |
39 | exit(1); | |
40 | } | |
34 | – | |
35 | – | // grid has a default value (default=51), so it is always given |
36 | – | npts = args_info.grid_arg; |
41 | ||
42 | < | printf ("opening %s\n", shapeFileName); |
43 | < | shapeFile = fopen( shapeFileName, "r" ); |
42 | > | shape = new SHAPE(); |
43 | > | shape->readSHAPEfile(shapeFileName); |
44 | ||
45 | ||
46 | + | 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 | ||
54 | < | my_file = fopen("junk.dat", "w"); |
55 | < | if (my_file == NULL) { |
56 | < | (void) fprintf(stderr, "No File\n"); |
54 | > | outputFile = fopen(outputFileName, "w"); |
55 | > | |
56 | > | if (outputFile == NULL) { |
57 | > | (void) fprintf(stderr, "Unable to open outputFile %s!\n", outputFileName); |
58 | exit(8); | |
59 | } | |
60 | + | |
61 | + | // grid has a default value (default=51), so it is always given |
62 | + | npts = args_info.grid_arg; |
63 | ||
64 | < | npts = 51; |
64 | > | xmin = -10.0; |
65 | > | xmax = 10.0; |
66 | ||
67 | < | xmin = -5.0; |
68 | < | xmax = 5.0; |
67 | > | ymin = -10.0; |
68 | > | ymax = 10.0; |
69 | ||
70 | < | ymin = -5.0; |
71 | < | ymax = 5.0; |
70 | > | zmin = -10.0; |
71 | > | zmax = 10.0; |
72 | ||
73 | < | zmin = -5.0; |
74 | < | zmax = 5.0; |
73 | > | //sigmaProbe = 2.28; |
74 | > | //sProbe = 2.28; |
75 | > | //epsProbe = 0.020269601874; |
76 | > | sigmaProbe = 0.0; |
77 | > | sProbe = 0.0; |
78 | > | epsProbe = 1.0; |
79 | ||
60 | – | rup = 4.0; |
61 | – | ru = 3.35; |
62 | – | rlow = 2.75; |
63 | – | w0 = 0.07715; |
64 | – | |
80 | for (i = 0; i < npts; i++) { | |
81 | x = xmin + (xmax-xmin) * (double)i/(double)npts; | |
82 | ||
# | Line 72 | Line 87 | int main(int argc, char* argv[]){ | |
87 | z = zmin + (zmax-zmin) * (double)k/(double)npts; | |
88 | ||
89 | r = sqrt(x*x + y*y + z*z); | |
90 | < | theta = acos(z/r); |
90 | > | costheta = z/r; |
91 | phi = atan(y/x); | |
77 | – | |
78 | – | fp = pow((cos(theta) - 0.6),2) * pow((cos(theta) + 0.8),2); |
92 | ||
93 | < | sp = (rup + 2.0*r - 3.0*rlow) * pow((rup - r),2) / pow((rup -rlow),3); |
93 | > | sigmaShape = shape->getSigmaAt(costheta, phi); |
94 | > | sShape = shape->getSAt(costheta, phi); |
95 | > | epsShape = shape->getEpsAt(costheta, phi); |
96 | ||
97 | < | if (r > rup) sp = 0.0; |
83 | < | if (r < rlow) sp = 1.0; |
97 | > | // Lorentz-Berthellot combining rules: |
98 | ||
99 | < | wp = sp*(fp - w0); |
99 | > | sigma = (sigmaShape + sigmaProbe) / 2.0; |
100 | > | s = (sShape + sProbe) / 2.0; |
101 | > | eps = sqrt(epsShape * epsProbe); |
102 | ||
103 | < | f = sin(theta)*sin(2.0*theta)*cos(2.0*phi); |
88 | < | |
89 | < | s = (ru + 2.0*r - 3.0*rlow) * pow((ru - r),2) / pow((ru -rlow),3); |
103 | > | rTerm = s / (r - sigma + s); |
104 | ||
105 | < | if (r > ru) s = 0.0; |
106 | < | if (r < rlow) s = 1.0; |
107 | < | |
108 | < | w = f*s; |
109 | < | |
110 | < | tot = w + wp; |
111 | < | |
98 | < | write_size = fwrite(&tot,1,8,my_file); |
99 | < | |
105 | > | r6 = pow(rTerm, 6); |
106 | > | r12 = r6*r6; |
107 | > | |
108 | > | energy = 4.0 * eps * (r12 - r6); |
109 | > | |
110 | > | write_size = fwrite(&energy,1,8,outputFile); |
111 | > | |
112 | } | |
113 | } | |
114 | } | |
115 | < | fclose(my_file); |
115 | > | fclose(outputFile); |
116 | return(0); | |
117 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |