1 |
– |
#include "SRI.hpp" |
2 |
– |
#include "Atom.hpp" |
1 |
|
|
2 |
|
#include <math.h> |
3 |
|
#include <iostream> |
4 |
|
#include <stdlib.h> |
5 |
|
|
6 |
< |
void GhostBend::set_atoms( Atom &a, Atom &b ){ |
6 |
> |
#include <simError.h> |
7 |
> |
#include "SRI.hpp" |
8 |
> |
#include "Atom.hpp" |
9 |
|
|
10 |
+ |
|
11 |
+ |
|
12 |
+ |
GhostBend::GhostBend( Atom &a, Atom &b ){ |
13 |
+ |
|
14 |
|
c_p_a = &a; |
15 |
< |
c_p_b = &b; |
15 |
> |
|
16 |
> |
if( !b.isDirectional() ){ |
17 |
> |
|
18 |
> |
// if atom b is not directional, then bad things will happen |
19 |
> |
|
20 |
> |
sprintf( painCave.errMsg, |
21 |
> |
" Ghost Bend error: Atom # %d of type \"%s\" is not " |
22 |
> |
"directional.\n", |
23 |
> |
c_p_b->getIndex, |
24 |
> |
c_p_b->getType() ); |
25 |
> |
painCave.isFatal = 1; |
26 |
> |
simError(); |
27 |
> |
} |
28 |
|
|
29 |
+ |
atomB = ( DirectionalAtom* ) &b; |
30 |
+ |
|
31 |
|
c_potential_E = 0.0; |
32 |
+ |
|
33 |
+ |
#ifdef IS_MPI |
34 |
+ |
sprintf( checkPointMsg, |
35 |
+ |
"Succesfull creation of Ghost bend\n" ); |
36 |
+ |
MPIcheckPoint(); |
37 |
+ |
#endif // is_mpi |
38 |
+ |
|
39 |
|
} |
40 |
|
|
41 |
|
|
54 |
|
double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj; |
55 |
|
double u[3]; |
56 |
|
|
57 |
< |
dx = c_p_a->getX() - c_p_b->getX(); |
58 |
< |
dy = c_p_a->getY() - c_p_b->getY(); |
59 |
< |
dz = c_p_a->getZ() - c_p_b->getZ(); |
57 |
> |
dx = c_p_a->getX() - atomB->getX(); |
58 |
> |
dy = c_p_a->getY() - atomB->getY(); |
59 |
> |
dz = c_p_a->getZ() - atomB->getZ(); |
60 |
|
|
61 |
< |
c_p_b->getU(u); |
61 |
> |
atomB->getU(u); |
62 |
|
|
63 |
|
gx = u[0]; |
64 |
|
gy = u[1]; |
136 |
|
c_p_a->addFy(dayi); |
137 |
|
c_p_a->addFz(dazi); |
138 |
|
|
139 |
< |
c_p_b->addFx(daxj + daxk); |
140 |
< |
c_p_b->addFy(dayj + dayk); |
141 |
< |
c_p_b->addFz(dazj + dazk); |
139 |
> |
atomB->addFx(daxj + daxk); |
140 |
> |
atomB->addFy(dayj + dayk); |
141 |
> |
atomB->addFz(dazj + dazk); |
142 |
|
|
143 |
< |
c_p_b->addTx(gy*dazk - gz*dayk); |
144 |
< |
c_p_b->addTy(gz*daxk - gx*dazk); |
145 |
< |
c_p_b->addTz(gx*dayk - gy*daxk); |
143 |
> |
atomB->addTx(gy*dazk - gz*dayk); |
144 |
> |
atomB->addTy(gz*daxk - gx*dazk); |
145 |
> |
atomB->addTz(gx*dayk - gy*daxk); |
146 |
|
|
147 |
|
return; |
148 |
|
} |
149 |
+ |
|
150 |
+ |
void GhostBend::setConstants( double the_c1, double the_c2, double the_c3, |
151 |
+ |
double the_Th0 ){ |
152 |
+ |
c1 = the_c1; |
153 |
+ |
c2 = the_c2; |
154 |
+ |
c3 = the_c3; |
155 |
+ |
theta0 = the_Th0; |
156 |
+ |
} |
157 |
+ |
|
158 |
+ |
|
159 |
+ |
double GhostBend::bend_force( double theta ){ |
160 |
+ |
|
161 |
+ |
double dt, dt2; |
162 |
+ |
double force; |
163 |
+ |
|
164 |
+ |
dt = ( theta - theta0 ) * M_PI / 180.0; |
165 |
+ |
dt2 = dt * dt; |
166 |
+ |
|
167 |
+ |
c_potential_E = ( c1 * dt2 ) + ( c2 * dt ) + c3; |
168 |
+ |
force = -( ( 2.0 * c1 * dt ) + c2 ); |
169 |
+ |
return force; |
170 |
+ |
} |