Skip to content

Commit

Permalink
Solenoid: apply radiation after multipolar kick
Browse files Browse the repository at this point in the history
  • Loading branch information
@szymonlopaciuk
szymonlopaciuk committed Jun 13, 2024
1 parent 513d6ae commit ae07076
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 49 deletions.
2 changes: 1 addition & 1 deletion tests/test_slicing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -280,7 +280,7 @@ def test_slicing_strategy_matching():
if name == 'keep_this':
assert isinstance(element, xt.Quadrupole)
elif name == 'keep_drifts':
assert (element, xt.Drift)
assert isinstance(element, xt.Drift)
elif name.startswith('drift_'):
assert 'DriftSlice' in type(element).__name__
elif name == 'keep_thin':
Expand Down
2 changes: 1 addition & 1 deletion tests/test_spacecharge_in_ring.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ def test_ring_with_spacecharge(test_context, mode):
else:
assert isinstance(line[0], xf.SpaceChargeBiGaussian)

if mode is not 'frozen':
if mode != 'frozen':
assert line.iscollective
else:
assert not line.iscollective
Expand Down
60 changes: 53 additions & 7 deletions xtrack/beam_elements/elements_src/solenoid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,12 @@ void Solenoid_track_local_particle(SolenoidData el, LocalParticle* part0) {
#endif

#ifndef XTRACK_SOLENOID_NO_SYNRAD // does this do the right thing?
double dp_record_entry = 0.;
double dpx_record_entry = 0.;
double dpy_record_entry = 0.;
double dp_record_exit = 0.;
double dpx_record_exit = 0.;
double dpy_record_exit = 0.;
double dp_record_entry = 0.;
double dpx_record_entry = 0.;
double dpy_record_entry = 0.;
double dp_record_exit = 0.;
double dpx_record_exit = 0.;
double dpy_record_exit = 0.;
#endif

int64_t num_multipole_kicks = SolenoidData_get_num_multipole_kicks(el);
Expand All @@ -39,13 +39,59 @@ void Solenoid_track_local_particle(SolenoidData el, LocalParticle* part0) {

for (int ii = 0; ii < num_multipole_kicks; ii++) {
//start_per_particle_block (part0->part)
Solenoid_thick_track_single_particle(part, slice_length, ks, radiation_flag,
#ifndef XTRACK_SOLENOID_NO_SYNRAD
double const old_px = LocalParticle_get_px(part);
double const old_py = LocalParticle_get_py(part);
double const old_ax = LocalParticle_get_ax(part);
double const old_ay = LocalParticle_get_ay(part);
double const old_zeta = LocalParticle_get_zeta(part);
double const rvv = LocalParticle_get_rvv(part);
#endif

track_solenoid_thick_single_particle(part, slice_length, ks, radiation_flag,
&dp_record_entry, &dpx_record_entry, &dpy_record_entry,
&dp_record_exit, &dpx_record_exit, &dpy_record_exit);

track_multipolar_kick_bend(
part, order, inv_factorial_order, knl, ksl, factor_knl_ksl,
kick_weight, 0, 0, 0, 0);

#ifndef XTRACK_SOLENOID_NO_SYNRAD
double l_path, curv;
if (radiation_flag > 0 && length > 0){
double const new_ax = LocalParticle_get_ax(part);
double const new_ay = LocalParticle_get_ay(part);

double const old_px_mech = old_px - old_ax;
double const old_py_mech = old_py - old_ay;

double const new_px_mech = LocalParticle_get_px(part) - new_ax;
double const new_py_mech = LocalParticle_get_py(part) - new_ay;

double const dpx = new_px_mech - old_px_mech;
double const dpy = new_py_mech - old_py_mech;

curv = sqrt(dpx * dpx + dpy * dpy) / length;

// Path length for radiation
double const dzeta = LocalParticle_get_zeta(part) - old_zeta;
l_path = rvv * (length - dzeta);

LocalParticle_add_to_px(part, -new_ax);
LocalParticle_add_to_py(part, -new_ay);

if (radiation_flag == 1){
synrad_average_kick(part, curv, l_path,
dp_record_exit, dpx_record_exit, dpy_record_exit);
}
else if (radiation_flag == 2){
synrad_emit_photons(part, curv, l_path, NULL, NULL);
}

LocalParticle_add_to_px(part, new_ax);
LocalParticle_add_to_py(part, new_ay);
}
#endif
//end_per_particle_block
}

Expand Down
82 changes: 42 additions & 40 deletions xtrack/beam_elements/elements_src/track_solenoid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,11 @@
#define IS_ZERO(X) (fabs(X) < 1e-9)

/*gpufun*/
void Solenoid_thick_track_single_particle(
void track_solenoid_thick_single_particle(
LocalParticle* part,
double length,
double ks,
int64_t radiation_flag,
double* dp_record_entry, double* dpx_record_entry, double* dpy_record_entry,
double* dp_record_exit, double* dpx_record_exit, double* dpy_record_exit
int64_t radiation_flag
) {
const double sk = ks / 2;

Expand Down Expand Up @@ -74,64 +72,68 @@ void Solenoid_thick_track_single_particle(
double const new_ax = -0.5 * Bz * new_y * q0 * QELEM / P0_J;
double const new_ay = 0.5 * Bz * new_x * q0 * QELEM / P0_J;

LocalParticle_set_x(part, new_x);
LocalParticle_set_px(part, new_px);
LocalParticle_set_y(part, new_y);
LocalParticle_set_py(part, new_py);
LocalParticle_add_to_zeta(part, add_to_zeta);
LocalParticle_add_to_s(part, length);
LocalParticle_set_ax(part, new_ax);
LocalParticle_set_ay(part, new_ay);
}

/*gpufun*/
void Solenoid_thick_track_single_particle(
LocalParticle* part,
double length,
double ks,
int64_t radiation_flag,
double* dp_record_entry, double* dpx_record_entry, double* dpy_record_entry,
double* dp_record_exit, double* dpx_record_exit, double* dpy_record_exit
) {
#ifndef XTRACK_SOLENOID_NO_SYNRAD
double const old_px = LocalParticle_get_px(part);
double const old_py = LocalParticle_get_py(part);
double const old_ax = LocalParticle_get_ax(part);
double const old_ay = LocalParticle_get_ay(part);
double const old_zeta = LocalParticle_get_zeta(part);
double const rvv = LocalParticle_get_rvv(part);
#endif

track_solenoid_thick_single_particle(part, length, ks, radiation_flag);

#ifndef XTRACK_SOLENOID_NO_SYNRAD
double l_path, curv;
if (radiation_flag > 0 && length > 0){
double const new_ax = LocalParticle_get_ax(part);
double const new_ay = LocalParticle_get_ay(part);

double const old_ax = LocalParticle_get_ax(part);
double const old_ay = LocalParticle_get_ay(part);

double const old_px_mech = px - old_ax;
double const old_py_mech = py - old_ay;
double const old_px_mech = old_px - old_ax;
double const old_py_mech = old_py - old_ay;

double const new_px_mech = new_px - new_ax;
double const new_py_mech = new_py - new_ay;
double const new_px_mech = LocalParticle_get_px(part) - new_ax;
double const new_py_mech = LocalParticle_get_py(part) - new_ay;

double const dpx = new_px_mech - old_px_mech;
double const dpy = new_py_mech - old_py_mech;

curv = sqrt(dpx*dpx + dpy*dpy) / length;
curv = sqrt(dpx * dpx + dpy * dpy) / length;

// Path length for radiation
double const dzeta = add_to_zeta;
double const rvv = LocalParticle_get_rvv(part);
double const dzeta = LocalParticle_get_zeta(part) - old_zeta;
l_path = rvv * (length - dzeta);

// LocalParticle_add_to_px(part, -old_ax);
// LocalParticle_add_to_py(part, -old_ay);
// if (radiation_flag == 1){
// synrad_average_kick(part, curv, l_path / 2,
// dp_record_entry, dpx_record_entry, dpy_record_entry);
// }
// else if (radiation_flag == 2){
// synrad_emit_photons(part, curv, l_path / 2, NULL, NULL);
// }
// LocalParticle_add_to_px(part, old_ax);
// LocalParticle_add_to_py(part, old_ay);
}
#endif

LocalParticle_set_x(part, new_x);
LocalParticle_add_to_px(part, new_px - px);
LocalParticle_set_y(part, new_y);
LocalParticle_add_to_py(part, new_py - py);
LocalParticle_add_to_zeta(part, add_to_zeta);
LocalParticle_add_to_s(part, length);
LocalParticle_set_ax(part, new_ax);
LocalParticle_set_ay(part, new_ay);

#ifndef XTRACK_SOLENOID_NO_SYNRAD

if (radiation_flag > 0 && length > 0){
LocalParticle_add_to_px(part, -new_ax);
LocalParticle_add_to_py(part, -new_ay);

if (radiation_flag == 1){
synrad_average_kick(part, curv, l_path,
dp_record_exit, dpx_record_exit, dpy_record_exit);
}
else if (radiation_flag == 2){
synrad_emit_photons(part, curv, l_path, NULL, NULL);
}

LocalParticle_add_to_px(part, new_ax);
LocalParticle_add_to_py(part, new_ay);
}
Expand Down

0 comments on commit ae07076

Please sign in to comment.