Commit a3c52278 authored by Daniel Winklehner's avatar Daniel Winklehner
Browse files

Added overloaded space charge for cyclotron function -DW

parent c0afc394
......@@ -643,7 +643,7 @@ void PartBunch::calcGammas_cycl() {
bingamma_m[i] /= pbin_m->getTotalNumPerBin(i);
else
bingamma_m[i] = 0.0;
INFOMSG("Bin " << i << " : particle number=" << pbin_m->getTotalNumPerBin(i) << " gamma = " << bingamma_m[i] << endl);
INFOMSG("Bin " << i << " : particle number = " << pbin_m->getTotalNumPerBin(i) << " gamma = " << bingamma_m[i] << endl);
}
}
......@@ -1109,7 +1109,81 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
IpplTimings::stopTimer(selfFieldTimer_m);
}
void PartBunch::computeSelfFields_cycl(double gamma, Vector_t const meanR, Vektor<double, 4> const quaternion) {
IpplTimings::startTimer(selfFieldTimer_m);
/// set initial charge density to zero.
rho_m = 0.0;
/// set initial E field to zero
eg_m = Vector_t(0.0);
if(fs_m->hasValidSolver()) {
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
/// from charge to charge density.
double tmp2 = 1.0 / gamma / (hr_m[0] * hr_m[1] * hr_m[2]);
rho_m *= tmp2;
/// Lorentz transformation
/// In particle rest frame, the longitudinal length enlarged
Vector_t hr_scaled = hr_m ;
hr_scaled[1] *= gamma;
/// now charge density is in rho_m
/// calculate Possion equation (without coefficient: -1/(eps))
fs_m->solver_m->computePotential(rho_m, hr_scaled);
/// additional work of FFT solver
/// now the scalar potential is given back in rho_m
rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
/// retrive coefficient: -1/(eps)
rho_m *= getCouplingConstant();
/// calculate electric field vectors from field potential
eg_m = -Grad(rho_m, eg_m);
/// back Lorentz transformation
eg_m *= Vector_t(gamma, 1.0, gamma);
/*
//debug
// output field on the grid points
int m1 = (int)nr_m[0]-1;
int m2 = (int)nr_m[0]/2;
for (int i=0; i<m1; i++ )
*gmsg << "Field along x axis E = " << eg_m[i][m2][m2] << " Pot = " << rho_m[i][m2][m2] << endl;
for (int i=0; i<m1; i++ )
*gmsg << "Field along y axis E = " << eg_m[m2][i][m2] << " Pot = " << rho_m[m2][i][m2] << endl;
for (int i=0; i<m1; i++ )
*gmsg << "Field along z axis E = " << eg_m[m2][m2][i] << " Pot = " << rho_m[m2][m2][i] << endl;
// end debug
*/
/// interpolate electric field at particle positions.
Ef.gather(eg_m, this->R, IntrplCIC_t());
/// calculate coefficient
double betaC = sqrt(gamma * gamma - 1.0) / gamma / Physics::c;
/// calculate B field from E field
Bf(0) = betaC * Ef(2);
Bf(2) = -betaC * Ef(0);
}
// *gmsg<<"gamma ="<<gamma<<endl;
// *gmsg<<"dx,dy,dz =("<<hr_m[0]<<", "<<hr_m[1]<<", "<<hr_m[2]<<") [m] "<<endl;
// *gmsg<<"max of bunch is ("<<rmax_m(0)<<", "<<rmax_m(1)<<", "<<rmax_m(2)<<") [m] "<<endl;
// *gmsg<<"min of bunch is ("<<rmin_m(0)<<", "<<rmin_m(1)<<", "<<rmin_m(2)<<") [m] "<<endl;
IpplTimings::stopTimer(selfFieldTimer_m);
}
void PartBunch::computeSelfFields_cycl(int bin) {
IpplTimings::startTimer(selfFieldTimer_m);
......@@ -1191,6 +1265,19 @@ void PartBunch::computeSelfFields_cycl(int bin) {
IpplTimings::stopTimer(selfFieldTimer_m);
}
inline void PartBunch::rotateWithQuaternion(Vector_t v, Vektor<double, 4> const quaternion) {
// rotates a Vector_t (3 elements) counter-clockwise using the quaternion 'quaternion'.
// Flip direction of rotation by quaternionVectorcomponent *= -1
Vector_t const quaternionVectorComponent = Vector_t(quaternion(1), quaternion(2), quaternion(3));
double const quaternionScalarComponent = quaternion(0);
v = 2.0f * dot(quaternionVectorComponent, v) * quaternionVectorComponent
+ (quaternionScalarComponent * quaternionScalarComponent -
dot(quaternionVectorComponent, quaternionVectorComponent)) * v
+ 2.0f * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}
void PartBunch::setBCAllOpen() {
for(int i = 0; i < 2 * 3; ++i) {
bc_m[i] = new ZeroFace<double, 3, Mesh_t, Center_t>(i);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment