Skip to content

Commit

Permalink
CFbend: Formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
ax3l committed Jun 27, 2023
1 parent 4b6dfab commit 1956888
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 51 deletions.
99 changes: 49 additions & 50 deletions src/particles/elements/CFbend.H
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace impactx
using PType = ImpactXParticleContainer::ParticleType;

/** An combined-function bend, consisting of an ideal sector bend with
* an upright quadrupole focusing component.
* an upright quadrupole focusing component.
*
* @param ds Segment length in m.
* @param rc Radius of curvature in m.
Expand All @@ -42,8 +42,11 @@ namespace impactx
* k < 0 horizontal defocusing
* @param nslice number of slices used for the application of space charge
*/
CFbend( amrex::ParticleReal const ds, amrex::ParticleReal const rc, amrex::ParticleReal const k,
int const nslice)
CFbend(
amrex::ParticleReal const ds,
amrex::ParticleReal const rc,
amrex::ParticleReal const k,
int const nslice)
: Thick(ds, nslice), m_rc(rc), m_k(k)
{
}
Expand All @@ -65,8 +68,8 @@ namespace impactx
amrex::ParticleReal & AMREX_RESTRICT px,
amrex::ParticleReal & AMREX_RESTRICT py,
amrex::ParticleReal & AMREX_RESTRICT pt,
RefPart const & refpart) const {

RefPart const & refpart) const
{
using namespace amrex::literals; // for _rt and _prt

// access AoS data such as positions and cpu/id
Expand All @@ -91,75 +94,72 @@ namespace impactx
amrex::ParticleReal const gx = m_k + pow(m_rc,-2);
amrex::ParticleReal const omegax = sqrt(std::abs(gx));

if(gx > 0.0) {
// calculate expensive terms once
amrex::ParticleReal const sinx = sin(omegax * slice_ds);
amrex::ParticleReal const cosx = cos(omegax * slice_ds);
amrex::ParticleReal const r56 = slice_ds/betgam2
+ (sinx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2));

// advance position and momentum (focusing)
p.pos(RealAoS::x) = cosx*x + sinx/omegax*px - (1.0_prt - cosx)/(gx*bet*m_rc)*pt;
pxout = -omegax*sinx*x + cosx*px - sinx/(omegax*bet*m_rc)*pt;
if(gx > 0.0) {
// calculate expensive terms once
amrex::ParticleReal const sinx = sin(omegax * slice_ds);
amrex::ParticleReal const cosx = cos(omegax * slice_ds);
amrex::ParticleReal const r56 = slice_ds/betgam2
+ (sinx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2));

p.pos(RealAoS::t) = sinx/(omegax*bet*m_rc)*x + (1.0_prt - cosx)/(gx*bet*m_rc)*px
+ t + r56*pt;
ptout = pt;
// advance position and momentum (focusing)
p.pos(RealAoS::x) = cosx*x + sinx/omegax*px - (1.0_prt - cosx)/(gx*bet*m_rc)*pt;
pxout = -omegax*sinx*x + cosx*px - sinx/(omegax*bet*m_rc)*pt;

p.pos(RealAoS::t) = sinx/(omegax*bet*m_rc)*x + (1.0_prt - cosx)/(gx*bet*m_rc)*px
+ t + r56*pt;
ptout = pt;
} else {
// calculate expensive terms once
amrex::ParticleReal const sinhx = sinh(omegax * slice_ds);
amrex::ParticleReal const coshx = cosh(omegax * slice_ds);
amrex::ParticleReal const r56 = slice_ds/betgam2
+ (sinhx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2));

// advance position and momentum (defocusing)
p.pos(RealAoS::x) = coshx*x + sinhx/omegax*px - (1.0_prt - coshx)/(gx*bet*m_rc)*pt;
pxout = omegax*sinhx*x + coshx*px - sinhx/(omegax*bet*m_rc)*pt;

p.pos(RealAoS::t) = sinhx/(omegax*bet*m_rc)*x + (1.0_prt - coshx)/(gx*bet*m_rc)*px
+ t + r56*pt;
ptout = pt;

}
// calculate expensive terms once
amrex::ParticleReal const sinhx = sinh(omegax * slice_ds);
amrex::ParticleReal const coshx = cosh(omegax * slice_ds);
amrex::ParticleReal const r56 = slice_ds/betgam2
+ (sinhx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2));

// advance position and momentum (defocusing)
p.pos(RealAoS::x) = coshx*x + sinhx/omegax*px - (1.0_prt - coshx)/(gx*bet*m_rc)*pt;
pxout = omegax*sinhx*x + coshx*px - sinhx/(omegax*bet*m_rc)*pt;

p.pos(RealAoS::t) = sinhx/(omegax*bet*m_rc)*x + (1.0_prt - coshx)/(gx*bet*m_rc)*px
+ t + r56*pt;
ptout = pt;
}

// update vertical phase space variables
amrex::ParticleReal const gy = -m_k;
amrex::ParticleReal const omegay = sqrt(std::abs(gy));

if(gy > 0.0) {
// calculate expensive terms once
amrex::ParticleReal const siny = sin(omegay * slice_ds);
amrex::ParticleReal const cosy = cos(omegay * slice_ds);
if(gy > 0.0) {
// calculate expensive terms once
amrex::ParticleReal const siny = sin(omegay * slice_ds);
amrex::ParticleReal const cosy = cos(omegay * slice_ds);

// advance position and momentum (focusing)
p.pos(RealAoS::y) = cosy*y + siny/omegay*py;
pyout = -omegay*siny*y + cosy*py;
// advance position and momentum (focusing)
p.pos(RealAoS::y) = cosy*y + siny/omegay*py;
pyout = -omegay*siny*y + cosy*py;

} else {
// calculate expensive terms once
amrex::ParticleReal const sinhy = sinh(omegay * slice_ds);
amrex::ParticleReal const coshy = cosh(omegay * slice_ds);
// calculate expensive terms once
amrex::ParticleReal const sinhy = sinh(omegay * slice_ds);
amrex::ParticleReal const coshy = cosh(omegay * slice_ds);

// advance position and momentum (defocusing)
p.pos(RealAoS::y) = coshy*y + sinhy/omegay*py;
pyout = omegay*sinhy*y + coshy*py;
// advance position and momentum (defocusing)
p.pos(RealAoS::y) = coshy*y + sinhy/omegay*py;
pyout = omegay*sinhy*y + coshy*py;
}

// assign updated momenta
px = pxout;
py = pyout;
pt = ptout;

}

/** This pushes the reference particle.
*
* @param[in,out] refpart reference particle
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void operator() (RefPart & AMREX_RESTRICT refpart) const {

void operator() (RefPart & AMREX_RESTRICT refpart) const
{
using namespace amrex::literals; // for _rt and _prt

// assign input reference particle values
Expand Down Expand Up @@ -198,7 +198,6 @@ namespace impactx

// advance integrated path length
refpart.s = s + slice_ds;

}

private:
Expand Down
2 changes: 1 addition & 1 deletion src/python/elements.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -292,7 +292,7 @@ void init_elements(py::module& m)
.def(py::init<
amrex::ParticleReal const,
amrex::ParticleReal const,
amrex::ParticleReal const,
amrex::ParticleReal const,
int const>(),
py::arg("ds"), py::arg("rc"), py::arg("k"), py::arg("nslice") = 1,
"An ideal combined function bend (sector bend with quadrupole component)."
Expand Down

0 comments on commit 1956888

Please sign in to comment.