| 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 |
+ |
} |
