# | 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; |
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* my_file; |
28 | > | FILE* outputFile; |
29 | ||
30 | //parse the command line options | |
31 | if (cmdline_parser (argc, argv, &args_info) != 0) | |
# | Line 37 | Line 42 | int main(int argc, char* argv[]){ | |
42 | shape = new SHAPE(); | |
43 | shape->readSHAPEfile(shapeFileName); | |
44 | ||
40 | – | |
41 | – | // grid has a default value (default=51), so it is always given |
42 | – | npts = args_info.grid_arg; |
45 | ||
46 | < | |
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 | + | outputFile = fopen(outputFileName, "w"); |
55 | ||
56 | < | my_file = fopen("junk.dat", "w"); |
57 | < | if (my_file == NULL) { |
49 | < | (void) fprintf(stderr, "No File\n"); |
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 | ||
64 | – | rup = 4.0; |
65 | – | ru = 3.35; |
66 | – | rlow = 2.75; |
67 | – | w0 = 0.07715; |
68 | – | |
80 | for (i = 0; i < npts; i++) { | |
81 | < | x = xmin + (xmax-xmin) * (double)i/(double)npts; |
81 | > | x = xmin + (xmax-xmin) * (double)i/(double)(npts-1); |
82 | ||
83 | for (j = 0; j < npts; j++) { | |
84 | < | y = ymin + (ymax-ymin) * (double)j/(double)npts; |
84 | > | y = ymin + (ymax-ymin) * (double)j/(double)(npts-1); |
85 | ||
86 | for (k = 0; k < npts; k++) { | |
87 | < | z = zmin + (zmax-zmin) * (double)k/(double)npts; |
87 | > | z = zmin + (zmax-zmin) * (double)k/(double)(npts-1); |
88 | ||
89 | r = sqrt(x*x + y*y + z*z); | |
90 | < | theta = acos(z/r); |
91 | < | phi = atan(y/x); |
81 | < | |
82 | < | fp = pow((cos(theta) - 0.6),2) * pow((cos(theta) + 0.8),2); |
90 | > | costheta = z/r; |
91 | > | phi = atan2(y,x); |
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; |
87 | < | 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); |
92 | < | |
93 | < | 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 | < | |
102 | < | write_size = fwrite(&tot,1,8,my_file); |
103 | < | |
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 |