From 2346e962041cbfd4b92dff851587fc5ba1d482e0 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Wed, 29 Oct 2025 15:50:33 -0300 Subject: [PATCH 01/24] add field3d element --- include/trackcpp/auxiliary.h | 8 ++++++-- include/trackcpp/elements.h | 7 ++++++- src/elements.cpp | 10 ++++++++++ 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 3cdfdc5..43edb5d 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -37,7 +37,8 @@ class PassMethodsClass { static const int pm_kickmap_pass = 8; static const int pm_matrix_pass = 9; static const int pm_drift_g2l_pass = 10; - static const int pm_nr_pms = 11; // counter for number of passmethods + static const int pm_field3d_pass = 11; + static const int pm_nr_pms = 12; // counter for number of passmethods PassMethodsClass() { passmethods.push_back("identity_pass"); passmethods.push_back("drift_pass"); @@ -50,6 +51,7 @@ class PassMethodsClass { passmethods.push_back("kicktable_pass"); passmethods.push_back("matrix_pass"); passmethods.push_back("drift_g2l_pass"); + passmethods.push_back("field3d_pass"); } int size() const { return passmethods.size(); } std::string operator[](const int i) const { return passmethods[i]; } @@ -71,7 +73,8 @@ struct PassMethod { pm_kickmap_pass = 8, pm_matrix_pass = 9, pm_drift_g2l_pass = 10, - pm_nr_pms = 11, + pm_field3d_pass = 11, + pm_nr_pms = 12, }; }; @@ -89,6 +92,7 @@ const std::vector pm_dict = { "kicktable_pass", "matrix_pass", "drift_g2l_pass", + "pm_field3d_pass", }; struct RadiationState { diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 94942d9..3639fec 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -63,7 +63,10 @@ class Element { double phase_lag = 0; // [rad] int kicktable_idx = -1; // index of kickmap object in kicktable_list double rescale_kicks = 1.0; // for kickmaps - + double ks = 0; + double kx = 0; + std::vector> coefs = std::vector>(5, std::vector(5, 0.0));; + std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; Matrix matrix66 = Matrix(6); @@ -111,6 +114,7 @@ class Element { static Element sextupole (const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); + static Element field3d (const std::string& fam_name_, const double& length_, const std::vector>& coefs_, const double& step_size_ = 0.2); bool operator==(const Element& o) const; bool operator!=(const Element& o) const { return !(*this == o); }; @@ -132,5 +136,6 @@ void initialize_quadrupole(Element& element, const double& K, const int& nr_step void initialize_sextupole(Element& element, const double& S, const int& nr_steps); void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); +void initialize_field3d(Element& element, const std::vector>& coefs_, const double& step_size_ = 0.2); #endif diff --git a/src/elements.cpp b/src/elements.cpp index 5d96fb1..9ab8ad7 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -156,6 +156,12 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt return e; } +Element Element::field3d (const std::string& fam_name_, const double& length_, const std::vector>& coefs_, const double& step_size_ = 0.2) { + Element e = Element(fam_name_, length_); + initialize_field3d(e, coefs_, step_size_); + return e; +} + void print_polynom(std::ostream& out, const std::string& label, const std::vector& polynom) { int order = 0; @@ -308,3 +314,7 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n element.kicktable_idx = kicktable_idx; element.rescale_kicks = rescale_kicks; } + +void initialize_field3d(Element& element, const std::vector>& coefs_, const double& step_size_) { + element.pass_method = PassMethod::pm_kickmap_pass; +} From ed79562821c97d145719846fbefc4004eff994cb Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 30 Oct 2025 07:31:34 -0300 Subject: [PATCH 02/24] update elements with kx and ks --- include/trackcpp/elements.h | 7 +++++-- src/elements.cpp | 12 ++++++++---- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 3639fec..f9c7988 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -65,6 +65,9 @@ class Element { double rescale_kicks = 1.0; // for kickmaps double ks = 0; double kx = 0; + // double E = 0; + // double P0 = 0; only for tests + // double beta0 = 0; std::vector> coefs = std::vector>(5, std::vector(5, 0.0));; std::vector polynom_a = default_polynom; @@ -114,7 +117,7 @@ class Element { static Element sextupole (const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); - static Element field3d (const std::string& fam_name_, const double& length_, const std::vector>& coefs_, const double& step_size_ = 0.2); + static Element field3d (const std::string& fam_name_, const double& length_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_ = 40); bool operator==(const Element& o) const; bool operator!=(const Element& o) const { return !(*this == o); }; @@ -136,6 +139,6 @@ void initialize_quadrupole(Element& element, const double& K, const int& nr_step void initialize_sextupole(Element& element, const double& S, const int& nr_steps); void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); -void initialize_field3d(Element& element, const std::vector>& coefs_, const double& step_size_ = 0.2); +void initialize_field3d(Element& element, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_); #endif diff --git a/src/elements.cpp b/src/elements.cpp index 9ab8ad7..03c9cba 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -156,9 +156,9 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt return e; } -Element Element::field3d (const std::string& fam_name_, const double& length_, const std::vector>& coefs_, const double& step_size_ = 0.2) { +Element Element::field3d (const std::string& fam_name_, const double& length_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_) { Element e = Element(fam_name_, length_); - initialize_field3d(e, coefs_, step_size_); + initialize_field3d(e, kx_, ks_, coefs_, nr_steps_); return e; } @@ -315,6 +315,10 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n element.rescale_kicks = rescale_kicks; } -void initialize_field3d(Element& element, const std::vector>& coefs_, const double& step_size_) { - element.pass_method = PassMethod::pm_kickmap_pass; +void initialize_field3d(Element& element, const double& kx, const double& ks, const std::vector>& coefs, const int nr_steps) { + element.pass_method = PassMethod::pm_field3d_pass; + element.nr_steps = nr_steps; + element.coefs = coefs; + element.kx = kx; + element.ks = ks; } From 224344c9ace85563fcf1dc708681d2a2b31f9ca2 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 30 Oct 2025 11:42:36 -0300 Subject: [PATCH 03/24] add s0 property to element --- include/trackcpp/elements.h | 12 +++++------- src/elements.cpp | 7 ++++--- src/field3d.cpp | 29 +++++++++++++++++++++++++++++ 3 files changed, 38 insertions(+), 10 deletions(-) create mode 100644 src/field3d.cpp diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index f9c7988..4735deb 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -63,11 +63,9 @@ class Element { double phase_lag = 0; // [rad] int kicktable_idx = -1; // index of kickmap object in kicktable_list double rescale_kicks = 1.0; // for kickmaps - double ks = 0; - double kx = 0; - // double E = 0; - // double P0 = 0; only for tests - // double beta0 = 0; + double ks = 0; // [1/m] + double kx = 0; // [1/m] + double s_init = 0; // [m] std::vector> coefs = std::vector>(5, std::vector(5, 0.0));; std::vector polynom_a = default_polynom; @@ -117,7 +115,7 @@ class Element { static Element sextupole (const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); - static Element field3d (const std::string& fam_name_, const double& length_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_ = 40); + static Element field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_ = 40); bool operator==(const Element& o) const; bool operator!=(const Element& o) const { return !(*this == o); }; @@ -139,6 +137,6 @@ void initialize_quadrupole(Element& element, const double& K, const int& nr_step void initialize_sextupole(Element& element, const double& S, const int& nr_steps); void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); -void initialize_field3d(Element& element, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_); +void initialize_field3d(Element& element, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_); #endif diff --git a/src/elements.cpp b/src/elements.cpp index 03c9cba..b5c7dc7 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -156,9 +156,9 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt return e; } -Element Element::field3d (const std::string& fam_name_, const double& length_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_) { +Element Element::field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_) { Element e = Element(fam_name_, length_); - initialize_field3d(e, kx_, ks_, coefs_, nr_steps_); + initialize_field3d(e, s0_, kx_, ks_, coefs_, nr_steps_); return e; } @@ -315,10 +315,11 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n element.rescale_kicks = rescale_kicks; } -void initialize_field3d(Element& element, const double& kx, const double& ks, const std::vector>& coefs, const int nr_steps) { +void initialize_field3d(Element& element, const double& s0, const double& kx, const double& ks, const std::vector>& coefs, const int nr_steps) { element.pass_method = PassMethod::pm_field3d_pass; element.nr_steps = nr_steps; element.coefs = coefs; element.kx = kx; element.ks = ks; + element.s_init = s0; } diff --git a/src/field3d.cpp b/src/field3d.cpp new file mode 100644 index 0000000..05c5653 --- /dev/null +++ b/src/field3d.cpp @@ -0,0 +1,29 @@ +// TRACKCPP - Particle tracking code +// Copyright (C) 2015 LNLS Accelerator Physics Group +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef _FIELD3D_H +#define _FIELD3D_H + +#include +#include +#include +#include +#include + + + + +#endif From 836294d30035ee3f3430cb88ecb5c2c04da0c8b9 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 30 Oct 2025 11:53:53 -0300 Subject: [PATCH 04/24] add field3d.h and field3d.cpp --- include/trackcpp/field3d.h | 98 +++++++++++++++++++ src/field3d.cpp | 191 ++++++++++++++++++++++++++++++++++++- 2 files changed, 286 insertions(+), 3 deletions(-) create mode 100644 include/trackcpp/field3d.h diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h new file mode 100644 index 0000000..ce69081 --- /dev/null +++ b/include/trackcpp/field3d.h @@ -0,0 +1,98 @@ +// TRACKCPP - Particle tracking code +// Copyright (C) 2015 LNLS Accelerator Physics Group +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef _FIELD3D_H +#define _FIELD3D_H + +#include "auxiliary.h" +#include +#include + + +template +T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) + + +template +T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) + + +template +T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) + + +template +T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) + + +template +T calc_D(const double& beta0, const T& delta) + +template +void exp_h1_z(const double& beta0, Pos& map, double step) + + +void exp_h1_s(double& s, double step) + +template +void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) + + +template +void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) + +template +void exp_h2_y(const double& beta0, Pos& map, double step) + + +template +void exp_h2_z(const double& beta0, Pos& map, double step) + + +template +void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) + + +template +void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) + + +template +void exp_h3_x(const double& beta0, Pos& map, double step) + +template +void exp_h3_z(const double& beta0, Pos& map, double step) + +template +void prop_h1(const double& beta0, Pos& map, double& s, double step) + +template +void prop_h2(const double& beta0, Pos& map, double step) + +template +void prop_h3(const double& beta0, Pos& map, double step) + +template +void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) + +template +void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) + +template +void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double& s, double step) + + +#endif diff --git a/src/field3d.cpp b/src/field3d.cpp index 05c5653..f6bcce8 100644 --- a/src/field3d.cpp +++ b/src/field3d.cpp @@ -14,9 +14,6 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -#ifndef _FIELD3D_H -#define _FIELD3D_H - #include #include #include @@ -24,6 +21,194 @@ #include +template +T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +{ + T ay_ = 0.0; + int M = coefs.size(); + int N = coefs[0].size(); + + for (int m = 1; m <= M; ++m) { + for (int n = 1; n <= N; ++n) { + double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); + double fac = coefs[m - 1][n - 1] * (m * kx) / (n * ks * ky); + ay_ += fac * sin(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); + } + } + + return -1 * ay_/brho; +} + +template +T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +{ + T ax_ = 0.0; + const int M = coefs.size(); + const int N = coefs[0].size(); + + for (int m = 1; m <= M; ++m) { + for (int n = 1; n <= N; ++n) { + double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); + double fac = coefs[m - 1][n - 1] / (n * ks); + ax_ += fac * cos(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); + } + } + + return -1 * ax_/brho; +} + +template +T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +{ + T day_dx = 0.0; + int M = coefs.size(); + int N = coefs[0].size(); + + for (int m = 1; m <= M; ++m) { + for (int n = 1; n <= N; ++n) { + double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); + double fac = coefs[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); + day_dx += fac * cos(m * kx * x)* std::cos(n * ks * s)* (cosh(ky * y) - 1.0); + } + } + + return -1 * day_dx/brho; +} + +template +T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +{ + T dax_dy = 0.0; + int M = coefs.size(); + int N = coefs[0].size(); + + for (int m = 1; m <= M; ++m) { + for (int n = 1; n <= N; ++n) { + double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); + double fac = coefs[m - 1][n - 1] * ky / (n * ks * m * kx); + dax_dy += fac * sin(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); + } + } + + return -1* dax_dy/brho; +} + +template +T calc_D(const double& beta0, const T& delta) { + return sqrt(1.0 + 2.0 * delta / beta0 + delta * delta); +} + +template +void exp_h1_z(const double& beta0, Pos& map, double step) { + T d = calc_D(beta0, map.de); + T factor = (1.0 / beta0 - (1.0 / beta0 + map.de) / d); + map.dl += factor * step / 2.0; +} + +void exp_h1_s(double& s, double step) { + s += step / 2.0; +} + +template +void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { + T factor = inty_day_dx(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + map.px += factor; +} + + +template +void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { + T factor = ay(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + map.py += factor; +} + +template +void exp_h2_y(const double& beta0, Pos& map, double step) { + T d = calc_D(beta0, map.de); + map.ry += map.py/d*step/2.0; +} + + +template +void exp_h2_z(const double& beta0, Pos& map, double step) { + T d = calc_D(beta0, map.de); + T factor = (pow(map.py, 2)*(1/beta0+map.de)/(2*pow(d, 3))); + map.dl -= factor*step/2; +} + + +template +void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { + T factor = ax(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + map.px += factor; +} + + +template +void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { + T factor = intx_dax_dy(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + map.py += factor; +} + + +template +void exp_h3_x(const double& beta0, Pos& map, double step) { + T d = calc_D(beta0, map.de); + map.rx += map.px/d*step; +} + +template +void exp_h3_z(const double& beta0, Pos& map, double step) { + T d = calc_D(beta0, map.de); + T factor = (pow(map.px, 2)*(1/beta0+map.de)/(2*pow(d, 3))); + map.dl -= factor*step; +} + +template +void prop_h1(const double& beta0, Pos& map, double& s, double step) { + exp_h1_z(beta0, map, step); + exp_h1_s(s, step); +} + +template +void prop_h2(const double& beta0, Pos& map, double step) { + exp_h2_y(beta0, map, step); + exp_h2_z(beta0, map, step); +} + +template +void prop_h3(const double& beta0, Pos& map, double step) { + exp_h3_x(beta0, map, step); + exp_h3_z(beta0, map, step); +} + +template +void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { + exp_ix_px(brho, kx, ks, coefs, map, s, sign, step); + exp_ix_py(brho, kx, ks, coefs, map, s, sign, step); + +} + +template +void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { + exp_iy_px(brho, kx, ks, coefs, map, s, sign, step); + exp_iy_py(brho, kx, ks, coefs, map, s, sign, step); + +} +template +void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double& s, double step) { + prop_h1(beta0, map, s, step); + prop_iy(brho, kx, ks, coefs, map, s, +1, step); + prop_h2(beta0, map, step); + prop_iy(brho, kx, ks, coefs, map, s, -1, step); + prop_ix(brho, kx, ks, coefs, map, s, +1, step); + prop_h3(beta0, map, step); + prop_ix(brho, kx, ks, coefs, map, s, -1, step); + prop_iy(brho, kx, ks, coefs, map, s, +1, step); + prop_h2(beta0, map, step); + prop_iy(brho, kx, ks, coefs, map, s, -1, step); + prop_h1(beta0, map, s, step); +} #endif From 4d6fc60ce2e072bfbeb66066b7a3271a8cd254cc Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 30 Oct 2025 13:57:34 -0300 Subject: [PATCH 05/24] add field3d passmethod --- include/trackcpp/accelerator.h | 1 + include/trackcpp/elements.h | 2 +- include/trackcpp/passmethods.h | 1 + include/trackcpp/passmethods.hpp | 17 +++++++++++++++++ src/elements.cpp | 2 +- 5 files changed, 21 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 0653507..044d425 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -18,6 +18,7 @@ #define _ACCELERATOR_H #include "kicktable.h" +#include "field3d.h" #include "elements.h" #include #include diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 4735deb..d75ace8 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -65,7 +65,7 @@ class Element { double rescale_kicks = 1.0; // for kickmaps double ks = 0; // [1/m] double kx = 0; // [1/m] - double s_init = 0; // [m] + double s0 = 0; // [m] std::vector> coefs = std::vector>(5, std::vector(5, 0.0));; std::vector polynom_a = default_polynom; diff --git a/include/trackcpp/passmethods.h b/include/trackcpp/passmethods.h index 51a38aa..35f9e67 100644 --- a/include/trackcpp/passmethods.h +++ b/include/trackcpp/passmethods.h @@ -66,6 +66,7 @@ template Status::type pm_thinsext_pass (Pos &pos, c template Status::type pm_kickmap_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_matrix_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_drift_g2l_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_field3d_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); #include "passmethods.hpp" diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 9b315ea..33d984c 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -113,6 +113,7 @@ Status::type kicktablethinkick(Pos& pos, const int& kicktable_idx, return status; } + template void matthinkick(Pos &pos, const Matrix &m) { @@ -468,6 +469,22 @@ Status::type pm_kickmap_pass(Pos &pos, const Element &elem, return status; } +template +Status::type pm_field3d_pass(Pos &pos, const Element &elem, + const Accelerator& accelerator) { + + global_2_local(pos, elem); + const double brho = get_magnetic_rigidity(accelerator.energy); + const double gamma = energy / electron_rest_energy_eV; + const double beta0 = sqrt(1 - 1/(gamma*gamma)); + double step = elem.length / float(elem.nr_steps); + for (int i=0; i Status::type pm_matrix_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { diff --git a/src/elements.cpp b/src/elements.cpp index b5c7dc7..1e3d748 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -321,5 +321,5 @@ void initialize_field3d(Element& element, const double& s0, const double& kx, co element.coefs = coefs; element.kx = kx; element.ks = ks; - element.s_init = s0; + element.s0 = s0; } From 6902e9ea38eafe0df769a2561655b6290c113fe3 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 30 Oct 2025 14:09:13 -0300 Subject: [PATCH 06/24] fix semicolons errors --- include/trackcpp/elements.h | 2 +- include/trackcpp/field3d.h | 44 ++++++++++++++++---------------- include/trackcpp/passmethods.hpp | 4 +-- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index d75ace8..5a00a6f 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -66,7 +66,7 @@ class Element { double ks = 0; // [1/m] double kx = 0; // [1/m] double s0 = 0; // [m] - std::vector> coefs = std::vector>(5, std::vector(5, 0.0));; + std::vector> coefs = std::vector>(5, std::vector(5, 0.0)); std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h index ce69081..73f8904 100644 --- a/include/trackcpp/field3d.h +++ b/include/trackcpp/field3d.h @@ -20,79 +20,79 @@ #include "auxiliary.h" #include #include - +#include "pos.h" template -T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); template -T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); template -T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); template -T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); template -T calc_D(const double& beta0, const T& delta) +T calc_D(const double& beta0, const T& delta); template -void exp_h1_z(const double& beta0, Pos& map, double step) +void exp_h1_z(const double& beta0, Pos& map, double step); -void exp_h1_s(double& s, double step) +void exp_h1_s(double& s, double step); template -void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) +void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); template -void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) +void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); template -void exp_h2_y(const double& beta0, Pos& map, double step) +void exp_h2_y(const double& beta0, Pos& map, double step); template -void exp_h2_z(const double& beta0, Pos& map, double step) +void exp_h2_z(const double& beta0, Pos& map, double step); template -void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) +void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); template -void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) +void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); template -void exp_h3_x(const double& beta0, Pos& map, double step) +void exp_h3_x(const double& beta0, Pos& map, double step); template -void exp_h3_z(const double& beta0, Pos& map, double step) +void exp_h3_z(const double& beta0, Pos& map, double step); template -void prop_h1(const double& beta0, Pos& map, double& s, double step) +void prop_h1(const double& beta0, Pos& map, double& s, double step); template -void prop_h2(const double& beta0, Pos& map, double step) +void prop_h2(const double& beta0, Pos& map, double step); template -void prop_h3(const double& beta0, Pos& map, double step) +void prop_h3(const double& beta0, Pos& map, double step); template -void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) +void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); template -void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) +void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); template -void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double& s, double step) +void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double& s, double step); #endif diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 33d984c..cb0baec 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -475,11 +475,11 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, global_2_local(pos, elem); const double brho = get_magnetic_rigidity(accelerator.energy); - const double gamma = energy / electron_rest_energy_eV; + const double gamma = accelerator.energy / electron_rest_energy_eV; const double beta0 = sqrt(1 - 1/(gamma*gamma)); double step = elem.length / float(elem.nr_steps); for (int i=0; i Date: Thu, 30 Oct 2025 15:41:38 -0300 Subject: [PATCH 07/24] change coefs in elements --- include/trackcpp/elements.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 5a00a6f..543fb72 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -65,8 +65,8 @@ class Element { double rescale_kicks = 1.0; // for kickmaps double ks = 0; // [1/m] double kx = 0; // [1/m] - double s0 = 0; // [m] - std::vector> coefs = std::vector>(5, std::vector(5, 0.0)); + double s0 = 0; // [m] + std::vector> coefs; std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; From 8b84a112ebc5ca0b6ef5fe13c99d7919f6c2edc3 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Fri, 31 Oct 2025 09:12:42 -0300 Subject: [PATCH 08/24] add field3d wrapper --- python_package/interface.cpp | 4 ++++ python_package/interface.h | 1 + 2 files changed, 5 insertions(+) diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 383c7fe..26fb6c6 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -237,6 +237,10 @@ Element kickmap_wrapper(const std::string& fam_name_, const std::string& kickta return Element::kickmap(fam_name_, kicktable_fname_, nr_steps_, rescale_length_, rescale_kicks_); } +Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_) { + return Element::field3d(fam_name_, length_, s0_, kx_, ks_, coefs_, nr_steps_); +} + Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag) { return read_flat_file(fname.data, accelerator, file_flag); } diff --git a/python_package/interface.h b/python_package/interface.h index e35e092..f0afa2c 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -91,6 +91,7 @@ Element quadrupole_wrapper(const std::string& fam_name_, const double& length_, Element sextupole_wrapper(const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); Element rfcavity_wrapper(const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag); Element kickmap_wrapper(const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length = 1.0, const double& rescale_kicks = 1.0); +Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_); Element rbend_wrapper(const std::string& fam_name_, const double& length_, const double& angle_, const double& angle_in_, const double& angle_out_, const double& gap_, const double& fint_in_, const double& fint_out_, From 8f5145c5172f84d956628da97a79d300b6a04d4f Mon Sep 17 00:00:00 2001 From: Gabriel Date: Fri, 31 Oct 2025 11:42:29 -0300 Subject: [PATCH 09/24] first compilation working --- include/trackcpp/accelerator.h | 1 - include/trackcpp/field3d.h | 4 ---- .../trackcpp/field3d.hpp | 18 ++++++++---------- include/trackcpp/passmethods.hpp | 4 +++- include/trackcpp/tracking.h | 3 +++ 5 files changed, 14 insertions(+), 16 deletions(-) rename src/field3d.cpp => include/trackcpp/field3d.hpp (95%) diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 044d425..0653507 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -18,7 +18,6 @@ #define _ACCELERATOR_H #include "kicktable.h" -#include "field3d.h" #include "elements.h" #include #include diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h index 73f8904..8fee55d 100644 --- a/include/trackcpp/field3d.h +++ b/include/trackcpp/field3d.h @@ -91,8 +91,4 @@ void prop_ix(const double& brho, const double& kx, const double& ks, const std:: template void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); -template -void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double& s, double step); - - #endif diff --git a/src/field3d.cpp b/include/trackcpp/field3d.hpp similarity index 95% rename from src/field3d.cpp rename to include/trackcpp/field3d.hpp index f6bcce8..caf598e 100644 --- a/src/field3d.cpp +++ b/include/trackcpp/field3d.hpp @@ -14,7 +14,6 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -#include #include #include #include @@ -105,19 +104,20 @@ void exp_h1_z(const double& beta0, Pos& map, double step) { map.dl += factor * step / 2.0; } -void exp_h1_s(double& s, double step) { +template +void exp_h1_s(T& s, T step) { s += step / 2.0; } template -void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { +void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { T factor = inty_day_dx(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.px += factor; } template -void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { +void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { T factor = ay(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.py += factor; } @@ -138,14 +138,14 @@ void exp_h2_z(const double& beta0, Pos& map, double step) { template -void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { +void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { T factor = ax(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.px += factor; } template -void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { +void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { T factor = intx_dax_dy(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.py += factor; } @@ -183,14 +183,14 @@ void prop_h3(const double& beta0, Pos& map, double step) { } template -void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { +void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { exp_ix_px(brho, kx, ks, coefs, map, s, sign, step); exp_ix_py(brho, kx, ks, coefs, map, s, sign, step); } template -void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign = 1, double step) { +void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { exp_iy_px(brho, kx, ks, coefs, map, s, sign, step); exp_iy_py(brho, kx, ks, coefs, map, s, sign, step); @@ -210,5 +210,3 @@ void prop_step(const double& beta0, const double& brho, const double& kx, const prop_iy(brho, kx, ks, coefs, map, s, -1, step); prop_h1(beta0, map, s, step); } - -#endif diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index cb0baec..41d5018 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -31,6 +31,7 @@ #include "auxiliary.h" #include "tpsa.h" #include "linalg.h" +#include "field3d.hpp" #include template inline T SQR(const T& X) { return X*X; } @@ -478,8 +479,9 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, const double gamma = accelerator.energy / electron_rest_energy_eV; const double beta0 = sqrt(1 - 1/(gamma*gamma)); double step = elem.length / float(elem.nr_steps); + double s0 = elem.s0; for (int i=0; i(orig_pos, el, accelerator)) != Status::success) return status; break; + case PassMethod::pm_field3d_pass: + if ((status = pm_field3d_pass(orig_pos, el, accelerator)) != Status::success) return status; + break; default: return Status::passmethod_not_defined; } From 56be44b1548ac4cb0baae0ae87fa355e3aa0f8a4 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Fri, 31 Oct 2025 13:28:01 -0300 Subject: [PATCH 10/24] fix power in field3d.hpp --- include/trackcpp/field3d.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index caf598e..85eef9c 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -132,7 +132,7 @@ void exp_h2_y(const double& beta0, Pos& map, double step) { template void exp_h2_z(const double& beta0, Pos& map, double step) { T d = calc_D(beta0, map.de); - T factor = (pow(map.py, 2)*(1/beta0+map.de)/(2*pow(d, 3))); + T factor = (map.py * map.py)*(1/beta0+map.de)/(2*d*d*d); map.dl -= factor*step/2; } @@ -160,7 +160,7 @@ void exp_h3_x(const double& beta0, Pos& map, double step) { template void exp_h3_z(const double& beta0, Pos& map, double step) { T d = calc_D(beta0, map.de); - T factor = (pow(map.px, 2)*(1/beta0+map.de)/(2*pow(d, 3))); + T factor = (map.px*map.px)*(1/beta0+map.de)/(2*d*d*d); map.dl -= factor*step; } From 866f0a4c6f9597484513a72bc2fa27e4c687cdc2 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Mon, 3 Nov 2025 07:50:22 -0300 Subject: [PATCH 11/24] calc D only once --- include/trackcpp/field3d.hpp | 43 ++++++++++++++------------------ include/trackcpp/passmethods.hpp | 3 ++- 2 files changed, 21 insertions(+), 25 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 85eef9c..fac1a72 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -98,8 +98,7 @@ T calc_D(const double& beta0, const T& delta) { } template -void exp_h1_z(const double& beta0, Pos& map, double step) { - T d = calc_D(beta0, map.de); +void exp_h1_z(const double& beta0, Pos& map, const T& d, double step) { T factor = (1.0 / beta0 - (1.0 / beta0 + map.de) / d); map.dl += factor * step / 2.0; } @@ -123,15 +122,13 @@ void exp_iy_py(const double& brho, const double& kx, const double& ks, const std } template -void exp_h2_y(const double& beta0, Pos& map, double step) { - T d = calc_D(beta0, map.de); +void exp_h2_y(const double& beta0, Pos& map, const T& d, double step) { map.ry += map.py/d*step/2.0; } template -void exp_h2_z(const double& beta0, Pos& map, double step) { - T d = calc_D(beta0, map.de); +void exp_h2_z(const double& beta0, Pos& map, const T& d, double step) { T factor = (map.py * map.py)*(1/beta0+map.de)/(2*d*d*d); map.dl -= factor*step/2; } @@ -152,34 +149,32 @@ void exp_ix_py(const double& brho, const double& kx, const double& ks, const std template -void exp_h3_x(const double& beta0, Pos& map, double step) { - T d = calc_D(beta0, map.de); +void exp_h3_x(const double& beta0, Pos& map, const T& d, double step) { map.rx += map.px/d*step; } template -void exp_h3_z(const double& beta0, Pos& map, double step) { - T d = calc_D(beta0, map.de); +void exp_h3_z(const double& beta0, Pos& map, const T& d, double step) { T factor = (map.px*map.px)*(1/beta0+map.de)/(2*d*d*d); map.dl -= factor*step; } template -void prop_h1(const double& beta0, Pos& map, double& s, double step) { - exp_h1_z(beta0, map, step); +void prop_h1(const double& beta0, Pos& map, const T& d, double& s, double step) { + exp_h1_z(beta0, map, d, step); exp_h1_s(s, step); } template -void prop_h2(const double& beta0, Pos& map, double step) { - exp_h2_y(beta0, map, step); - exp_h2_z(beta0, map, step); +void prop_h2(const double& beta0, Pos& map, const T& d, double step) { + exp_h2_y(beta0, map, d, step); + exp_h2_z(beta0, map, d, step); } template -void prop_h3(const double& beta0, Pos& map, double step) { - exp_h3_x(beta0, map, step); - exp_h3_z(beta0, map, step); +void prop_h3(const double& beta0, Pos& map, const T& d, double step) { + exp_h3_x(beta0, map, d, step); + exp_h3_z(beta0, map, d, step); } template @@ -197,16 +192,16 @@ void prop_iy(const double& brho, const double& kx, const double& ks, const std:: } template -void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double& s, double step) { - prop_h1(beta0, map, s, step); +void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, const T& d, double& s, double step) { + prop_h1(beta0, map, d, s, step); prop_iy(brho, kx, ks, coefs, map, s, +1, step); - prop_h2(beta0, map, step); + prop_h2(beta0, map, d, step); prop_iy(brho, kx, ks, coefs, map, s, -1, step); prop_ix(brho, kx, ks, coefs, map, s, +1, step); - prop_h3(beta0, map, step); + prop_h3(beta0, map, d, step); prop_ix(brho, kx, ks, coefs, map, s, -1, step); prop_iy(brho, kx, ks, coefs, map, s, +1, step); - prop_h2(beta0, map, step); + prop_h2(beta0, map, d, step); prop_iy(brho, kx, ks, coefs, map, s, -1, step); - prop_h1(beta0, map, s, step); + prop_h1(beta0, map, d, s, step); } diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 41d5018..6edf789 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -480,8 +480,9 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, const double beta0 = sqrt(1 - 1/(gamma*gamma)); double step = elem.length / float(elem.nr_steps); double s0 = elem.s0; + T d = calc_D(beta0, pos.de); for (int i=0; i Date: Tue, 4 Nov 2025 15:11:48 -0300 Subject: [PATCH 12/24] add more coefs to elements --- include/trackcpp/elements.h | 13 +++++++++++-- python_package/interface.cpp | 7 +++++-- python_package/interface.h | 5 ++++- src/elements.cpp | 15 ++++++++++++--- 4 files changed, 32 insertions(+), 8 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 543fb72..d657b30 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -67,6 +67,9 @@ class Element { double kx = 0; // [1/m] double s0 = 0; // [m] std::vector> coefs; + std::vector> coefs2; + std::vector> coefs3; + std::vector> coefs4; std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; @@ -115,7 +118,10 @@ class Element { static Element sextupole (const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); - static Element field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_ = 40); + static Element field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, + const std::vector>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs3_, const std::vector>& coefs4_, + const int nr_steps_ = 40); bool operator==(const Element& o) const; bool operator!=(const Element& o) const { return !(*this == o); }; @@ -137,6 +143,9 @@ void initialize_quadrupole(Element& element, const double& K, const int& nr_step void initialize_sextupole(Element& element, const double& S, const int& nr_steps); void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); -void initialize_field3d(Element& element, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_); +void initialize_field3d(Element& element, const double& s0_, const double& kx_, const double& ks_, + const std::vector>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs3_, const std::vector>& coefs4_, + const int nr_steps_); #endif diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 26fb6c6..47d3184 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -237,8 +237,11 @@ Element kickmap_wrapper(const std::string& fam_name_, const std::string& kickta return Element::kickmap(fam_name_, kicktable_fname_, nr_steps_, rescale_length_, rescale_kicks_); } -Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_) { - return Element::field3d(fam_name_, length_, s0_, kx_, ks_, coefs_, nr_steps_); +Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, + const std::vector>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs4_, const std::vector>& coefs3_, + const int nr_steps_) { + return Element::field3d(fam_name_, length_, s0_, kx_, ks_, coefs_, coefs2_, coefs3_, coefs4_, nr_steps_); } Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag) { diff --git a/python_package/interface.h b/python_package/interface.h index f0afa2c..4c9aba4 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -91,7 +91,10 @@ Element quadrupole_wrapper(const std::string& fam_name_, const double& length_, Element sextupole_wrapper(const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); Element rfcavity_wrapper(const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag); Element kickmap_wrapper(const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length = 1.0, const double& rescale_kicks = 1.0); -Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_); +Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, + const std::vector>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs3_, const std::vector>& coefs4_, + const int nr_steps_); Element rbend_wrapper(const std::string& fam_name_, const double& length_, const double& angle_, const double& angle_in_, const double& angle_out_, const double& gap_, const double& fint_in_, const double& fint_out_, diff --git a/src/elements.cpp b/src/elements.cpp index 1e3d748..6f44b04 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -156,9 +156,12 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt return e; } -Element Element::field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const int nr_steps_) { +Element Element::field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, + const std::vector>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs3_, const std::vector>& coefs4_, + const int nr_steps_) { Element e = Element(fam_name_, length_); - initialize_field3d(e, s0_, kx_, ks_, coefs_, nr_steps_); + initialize_field3d(e, s0_, kx_, ks_, coefs_, coefs2_, coefs3_, coefs4_,nr_steps_); return e; } @@ -315,10 +318,16 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n element.rescale_kicks = rescale_kicks; } -void initialize_field3d(Element& element, const double& s0, const double& kx, const double& ks, const std::vector>& coefs, const int nr_steps) { +void initialize_field3d(Element& element, const double& s0, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + const int nr_steps) { element.pass_method = PassMethod::pm_field3d_pass; element.nr_steps = nr_steps; element.coefs = coefs; + element.coefs2 = coefs2; + element.coefs3 = coefs3; + element.coefs4 = coefs4; element.kx = kx; element.ks = ks; element.s0 = s0; From 1d890617e46a45f218caa0b62acb40fc996f1ab1 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 4 Nov 2025 15:20:06 -0300 Subject: [PATCH 13/24] add more coefs to passmethod interface --- include/trackcpp/field3d.hpp | 5 ++++- include/trackcpp/passmethods.hpp | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index fac1a72..ab5a5b7 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -192,7 +192,10 @@ void prop_iy(const double& brho, const double& kx, const double& ks, const std:: } template -void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, const T& d, double& s, double step) { +void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + Pos& map, const T& d, double& s, double step) { prop_h1(beta0, map, d, s, step); prop_iy(brho, kx, ks, coefs, map, s, +1, step); prop_h2(beta0, map, d, step); diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 6edf789..752972c 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -482,7 +482,7 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, double s0 = elem.s0; T d = calc_D(beta0, pos.de); for (int i=0; i Date: Tue, 4 Nov 2025 15:32:20 -0300 Subject: [PATCH 14/24] add more coefs to prop py and px functions --- include/trackcpp/field3d.hpp | 50 ++++++++++++++++++++++++------------ 1 file changed, 34 insertions(+), 16 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index ab5a5b7..6787674 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -109,14 +109,20 @@ void exp_h1_s(T& s, T step) { } template -void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { +void exp_iy_px(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + Pos& map, double s, int sign, double step) { T factor = inty_day_dx(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.px += factor; } template -void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { +void exp_iy_py(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + Pos& map, double s, int sign, double step) { T factor = ay(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.py += factor; } @@ -135,14 +141,20 @@ void exp_h2_z(const double& beta0, Pos& map, const T& d, double step) { template -void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { +void exp_ix_px(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + Pos& map, double s, int sign, double step) { T factor = ax(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.px += factor; } template -void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { +void exp_ix_py(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + Pos& map, double s, int sign, double step) { T factor = intx_dax_dy(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; map.py += factor; } @@ -178,16 +190,22 @@ void prop_h3(const double& beta0, Pos& map, const T& d, double step) { } template -void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { - exp_ix_px(brho, kx, ks, coefs, map, s, sign, step); - exp_ix_py(brho, kx, ks, coefs, map, s, sign, step); +void prop_ix(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + Pos& map, double s, int sign, double step) { + exp_ix_px(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); + exp_ix_py(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); } template -void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step) { - exp_iy_px(brho, kx, ks, coefs, map, s, sign, step); - exp_iy_py(brho, kx, ks, coefs, map, s, sign, step); +void prop_iy(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + Pos& map, double s, int sign, double step) { + exp_iy_px(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); + exp_iy_py(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); } @@ -197,14 +215,14 @@ void prop_step(const double& beta0, const double& brho, const double& kx, const const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, const T& d, double& s, double step) { prop_h1(beta0, map, d, s, step); - prop_iy(brho, kx, ks, coefs, map, s, +1, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); prop_h2(beta0, map, d, step); - prop_iy(brho, kx, ks, coefs, map, s, -1, step); - prop_ix(brho, kx, ks, coefs, map, s, +1, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); prop_h3(beta0, map, d, step); - prop_ix(brho, kx, ks, coefs, map, s, -1, step); - prop_iy(brho, kx, ks, coefs, map, s, +1, step); + prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); prop_h2(beta0, map, d, step); - prop_iy(brho, kx, ks, coefs, map, s, -1, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); prop_h1(beta0, map, d, s, step); } From 912bcce8703ec828261d0ddfb657463b2636fdcc Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 4 Nov 2025 15:51:08 -0300 Subject: [PATCH 15/24] add more coeffs to potential vector --- include/trackcpp/field3d.hpp | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 6787674..3e3fd22 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -21,7 +21,10 @@ template -T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T ay(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + const T& x, const T& y, const double& s) { T ay_ = 0.0; int M = coefs.size(); @@ -39,7 +42,10 @@ T ay(const double& brho, const double& kx, const double& ks, const std::vector -T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T ax(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + const T& x, const T& y, const double& s) { T ax_ = 0.0; const int M = coefs.size(); @@ -57,7 +63,10 @@ T ax(const double& brho, const double& kx, const double& ks, const std::vector -T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T inty_day_dx(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + const T& x, const T& y, const double& s) { T day_dx = 0.0; int M = coefs.size(); @@ -75,7 +84,10 @@ T inty_day_dx(const double& brho, const double& kx, const double& ks, const std: } template -T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s) +T intx_dax_dy(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs3, const std::vector>& coefs4, + const T& x, const T& y, const double& s) { T dax_dy = 0.0; int M = coefs.size(); @@ -113,7 +125,7 @@ void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = inty_day_dx(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + T factor = inty_day_dx(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; map.px += factor; } @@ -123,7 +135,7 @@ void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = ay(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + T factor = ay(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; map.py += factor; } @@ -145,7 +157,7 @@ void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = ax(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + T factor = ax(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; map.px += factor; } @@ -155,7 +167,7 @@ void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = intx_dax_dy(brho, kx, ks, coefs, map.rx, map.ry, s) * sign * -1; + T factor = intx_dax_dy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; map.py += factor; } From 3cc34aadb893faf025b83c09c9a3958f429e1691 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 4 Nov 2025 16:10:38 -0300 Subject: [PATCH 16/24] working with all coefficients --- include/trackcpp/field3d.hpp | 48 +++++++++++++++++++++++++++--------- python_package/interface.cpp | 2 +- 2 files changed, 37 insertions(+), 13 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 3e3fd22..991f4ae 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -33,12 +33,18 @@ T ay(const double& brho, const double& kx, const double& ks, for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac = coefs[m - 1][n - 1] * (m * kx) / (n * ks * ky); - ay_ += fac * sin(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); + double fac1 = coefs[m - 1][n - 1] * (m * kx) / (n * ks * ky); + double fac2 = coefs2[m - 1][n - 1] * (m * kx) / (n * ks * ky); + double fac3 = -coefs3[m - 1][n - 1] * (m * kx) / (n * ks * ky); + double fac4 = -coefs4[m - 1][n - 1] * (m * kx) / (n * ks * ky); + ay_ += fac1 * sin(m * kx * x) * sinh(ky * y) * std::sin(n * ks * s); + ay_ += fac2 * cos(m * kx * x) * sinh(ky * y) * std::sin(n * ks * s); + ay_ += fac3 * sin(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); + ay_ += fac4 * cos(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); } } - return -1 * ay_/brho; + return 1 * ay_/brho; } template @@ -54,12 +60,18 @@ T ax(const double& brho, const double& kx, const double& ks, for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac = coefs[m - 1][n - 1] / (n * ks); - ax_ += fac * cos(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); + double fac1 = coefs[m - 1][n - 1] / (n * ks); + double fac2 = coefs2[m - 1][n - 1] / (n * ks); + double fac3 = -coefs3[m - 1][n - 1] / (n * ks); + double fac4 = -coefs4[m - 1][n - 1] / (n * ks); + ax_ += fac1 * cos(m * kx * x) * cosh(ky * y) * std::sin(n * ks * s); + ax_ += fac2 * sin(m * kx * x) * cosh(ky * y) * std::sin(n * ks * s); + ax_ += fac3 * cos(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); + ax_ += fac4 * sin(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); } } - return -1 * ax_/brho; + return 1 * ax_/brho; } template @@ -75,12 +87,18 @@ T inty_day_dx(const double& brho, const double& kx, const double& ks, for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac = coefs[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); - day_dx += fac * cos(m * kx * x)* std::cos(n * ks * s)* (cosh(ky * y) - 1.0); + double fac1 = coefs[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); + double fac2 = -coefs2[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); + double fac3 = -coefs3[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); + double fac4 = coefs4[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); + day_dx += fac1 * cos(m * kx * x) * std::sin(n * ks * s) * (cosh(ky * y) - 1.0); + day_dx += fac2 * sin(m * kx * x) * std::sin(n * ks * s) * (cosh(ky * y) - 1.0); + day_dx += fac3 * cos(m * kx * x) * std::cos(n * ks * s) * (cosh(ky * y) - 1.0); + day_dx += fac4 * sin(m * kx * x) * std::cos(n * ks * s) * (cosh(ky * y) - 1.0); } } - return -1 * day_dx/brho; + return 1 * day_dx/brho; } template @@ -96,12 +114,18 @@ T intx_dax_dy(const double& brho, const double& kx, const double& ks, for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac = coefs[m - 1][n - 1] * ky / (n * ks * m * kx); - dax_dy += fac * sin(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); + double fac1 = coefs[m - 1][n - 1] * ky / (n * ks * m * kx); + double fac2 = -coefs2[m - 1][n - 1] * ky / (n * ks * m * kx); + double fac3 = -coefs3[m - 1][n - 1] * ky / (n * ks * m * kx); + double fac4 = coefs4[m - 1][n - 1] * ky / (n * ks * m * kx); + dax_dy += fac1 * sin(m * kx * x) * std::sin(n * ks * s) * sinh(ky * y); + dax_dy += fac2 * cos(m * kx * x) * std::sin(n * ks * s) * sinh(ky * y); + dax_dy += fac3 * sin(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); + dax_dy += fac4 * cos(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); } } - return -1* dax_dy/brho; + return 1* dax_dy/brho; } template diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 47d3184..798e0e9 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -239,7 +239,7 @@ Element kickmap_wrapper(const std::string& fam_name_, const std::string& kickta Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, const std::vector>& coefs_, const std::vector>& coefs2_, - const std::vector>& coefs4_, const std::vector>& coefs3_, + const std::vector>& coefs3_, const std::vector>& coefs4_, const int nr_steps_) { return Element::field3d(fam_name_, length_, s0_, kx_, ks_, coefs_, coefs2_, coefs3_, coefs4_, nr_steps_); } From 8d4802a18ee2cc2b32ac4861d3a1ec67107ccef7 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Fri, 7 Nov 2025 13:22:54 -0300 Subject: [PATCH 17/24] change beta0 to 2 in field3d passmethod --- include/trackcpp/field3d.hpp | 3 ++- include/trackcpp/passmethods.hpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 991f4ae..8f35a89 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -130,7 +130,8 @@ T intx_dax_dy(const double& brho, const double& kx, const double& ks, template T calc_D(const double& beta0, const T& delta) { - return sqrt(1.0 + 2.0 * delta / beta0 + delta * delta); + // return sqrt(1.0 + 2.0 * delta / beta0 + delta * delta); + return 1.0 + delta; } template diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 752972c..a0374ea 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -477,7 +477,8 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, global_2_local(pos, elem); const double brho = get_magnetic_rigidity(accelerator.energy); const double gamma = accelerator.energy / electron_rest_energy_eV; - const double beta0 = sqrt(1 - 1/(gamma*gamma)); + // const double beta0 = sqrt(1 - 1/(gamma*gamma)); + const double beta0 = 1; double step = elem.length / float(elem.nr_steps); double s0 = elem.s0; T d = calc_D(beta0, pos.de); From 21bd85fd295d6538fb860c425258e10eb6d85475 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 18 Nov 2025 15:45:08 -0300 Subject: [PATCH 18/24] add id passmehotd flat_file --- include/trackcpp/auxiliary.h | 2 +- include/trackcpp/field3d.hpp | 68 +++++++++++++------------- include/trackcpp/passmethods.hpp | 19 ++++++-- src/flat_file.cpp | 83 ++++++++++++++++++++++++++++++++ 4 files changed, 134 insertions(+), 38 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 43edb5d..df2af36 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -92,7 +92,7 @@ const std::vector pm_dict = { "kicktable_pass", "matrix_pass", "drift_g2l_pass", - "pm_field3d_pass", + "field3d_pass", }; struct RadiationState { diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 8f35a89..23a6015 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -150,7 +150,7 @@ void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = inty_day_dx(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; + T factor = inty_day_dx(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } @@ -160,20 +160,21 @@ void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = ay(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; + T factor = ay(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } template -void exp_h2_y(const double& beta0, Pos& map, const T& d, double step) { - map.ry += map.py/d*step/2.0; +void exp_h2_y(const double& beta0, Pos& map, const T& pnorm, double step) { + map.ry += 0.5*step*pnorm*map.py; } template -void exp_h2_z(const double& beta0, Pos& map, const T& d, double step) { - T factor = (map.py * map.py)*(1/beta0+map.de)/(2*d*d*d); - map.dl -= factor*step/2; +void exp_h2_z(const double& beta0, Pos& map, const T& pnorm, double step) { + // T factor = (map.py * map.py)*(1/beta0+map.de)/(2*d*d*d); + // T factor = (map.py * map.py)/d/d*0.5; + map.dl += 0.25*step*pnorm*pnorm*(map.py * map.py); } @@ -182,7 +183,7 @@ void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = ax(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; + T factor = ax(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } @@ -192,38 +193,39 @@ void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = intx_dax_dy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1; + T factor = intx_dax_dy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } template -void exp_h3_x(const double& beta0, Pos& map, const T& d, double step) { - map.rx += map.px/d*step; +void exp_h3_x(const double& beta0, Pos& map, const T& pnorm, double step) { + map.rx += step*pnorm*map.px; } template -void exp_h3_z(const double& beta0, Pos& map, const T& d, double step) { - T factor = (map.px*map.px)*(1/beta0+map.de)/(2*d*d*d); - map.dl -= factor*step; +void exp_h3_z(const double& beta0, Pos& map, const T& pnorm, double step) { + // T factor = (map.px*map.px)*(1/beta0+map.de)/(2*d*d*d); + // T factor = (map.px*map.px)/d/d*0.5; + map.dl += 0.5*step*pnorm*pnorm*(map.px * map.px); } template void prop_h1(const double& beta0, Pos& map, const T& d, double& s, double step) { - exp_h1_z(beta0, map, d, step); + // exp_h1_z(beta0, map, d, step); exp_h1_s(s, step); } template -void prop_h2(const double& beta0, Pos& map, const T& d, double step) { - exp_h2_y(beta0, map, d, step); - exp_h2_z(beta0, map, d, step); +void prop_h2(const double& beta0, Pos& map, const T& pnorm, double step) { + exp_h2_y(beta0, map, pnorm, step); + exp_h2_z(beta0, map, pnorm, step); } template -void prop_h3(const double& beta0, Pos& map, const T& d, double step) { - exp_h3_x(beta0, map, d, step); - exp_h3_z(beta0, map, d, step); +void prop_h3(const double& beta0, Pos& map, const T& pnorm, double step) { + exp_h3_x(beta0, map, pnorm, step); + exp_h3_z(beta0, map, pnorm, step); } template @@ -250,16 +252,16 @@ template void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, - Pos& map, const T& d, double& s, double step) { - prop_h1(beta0, map, d, s, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); - prop_h2(beta0, map, d, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); - prop_h3(beta0, map, d, step); - prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); - prop_h2(beta0, map, d, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - prop_h1(beta0, map, d, s, step); + Pos& map, const T& pnorm, double& s, double step) { + // prop_h1(beta0, map, pnorm, s, step); + // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_h2(beta0, map, pnorm, step); + // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + // prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_h3(beta0, map, pnorm, step); + // prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_h2(beta0, map, pnorm, step); + // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + // prop_h1(beta0, map, pnorm, s, step); } diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index a0374ea..6ab4f5b 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -473,17 +473,28 @@ Status::type pm_kickmap_pass(Pos &pos, const Element &elem, template Status::type pm_field3d_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { - + global_2_local(pos, elem); const double brho = get_magnetic_rigidity(accelerator.energy); const double gamma = accelerator.energy / electron_rest_energy_eV; // const double beta0 = sqrt(1 - 1/(gamma*gamma)); - const double beta0 = 1; + const double beta0 = 1.0; double step = elem.length / float(elem.nr_steps); double s0 = elem.s0; - T d = calc_D(beta0, pos.de); + // T d = calc_D(beta0, pos.de); + T pnorm = 1 / (1 + pos.de); for (int i=0; i& p); +static bool has_coeffs(const std::vector>& matrix); static void write_6d_vector(std::ostream& fp, const std::string& label, const double* t); static void write_6d_vector(std::ostream& fp, const std::string& label, const std::vector& t); static void write_polynom(std::ostream& fp, const std::string& label, const std::vector& p); +static void write_coeffs(std::ostream& fp, const std::string& label, const std::vector>& matrix); static void synchronize_polynomials(Element& e); static void read_polynomials(std::ifstream& fp, Element& e); static void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator); @@ -127,6 +129,15 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) if (e.angle_in != 0) { fp << std::setw(pw) << "angle_in" << e.angle_in << '\n'; } if (e.angle_out != 0) { fp << std::setw(pw) << "angle_out" << e.angle_out << '\n'; } if (e.rescale_kicks != 1.0) { fp << std::setw(pw) << "rescale_kicks" << e.rescale_kicks << '\n'; } + if (e.kx != 0) { fp << std::setw(pw) << "kx" << e.kx << '\n'; } + if (e.ks != 0) { fp << std::setw(pw) << "ks" << e.ks << '\n'; } + if (e.s0 != 0) { fp << std::setw(pw) << "s0" << e.s0 << '\n'; } + if (has_coeffs(e.coefs)) write_coeffs(fp, "coefs", e.coefs); + if (has_coeffs(e.coefs2)) write_coeffs(fp, "coefs2", e.coefs2); + if (has_coeffs(e.coefs3)) write_coeffs(fp, "coefs3", e.coefs3); + if (has_coeffs(e.coefs4)) write_coeffs(fp, "coefs4", e.coefs4); + + if (e.has_t_in) write_6d_vector(fp, "t_in", e.t_in); if (e.has_t_out) write_6d_vector(fp, "t_out", e.t_out); if (e.has_r_in) { @@ -226,6 +237,9 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) if (cmd.compare("angle_in") == 0) { ss >> e.angle_in; continue; } if (cmd.compare("angle_out") == 0) { ss >> e.angle_out; continue; } if (cmd.compare("rescale_kicks") == 0) { ss >> e.rescale_kicks; continue; } + if (cmd.compare("kx") == 0) { ss >> e.kx; continue; } + if (cmd.compare("ks") == 0) { ss >> e.ks; continue; } + if (cmd.compare("s0") == 0) { ss >> e.s0; continue; } if (cmd.compare("t_in") == 0) { for(auto i=0; i<6; ++i) ss >> e.t_in[i]; e.reflag_t_in(); continue; } if (cmd.compare("t_out") == 0) { for(auto i=0; i<6; ++i) ss >> e.t_out[i]; e.reflag_t_out(); continue; } if (cmd.compare("rx|r_in") == 0) { for(auto i=0; i<6; ++i) ss >> e.r_in[0*6+i]; e.reflag_r_in(); continue; } @@ -303,6 +317,49 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) synchronize_polynomials(e); continue; } + + if (cmd.compare("coefs") == 0 || + cmd.compare("coefs2") == 0 || + cmd.compare("coefs3") == 0 || + cmd.compare("coefs4") == 0) { + std::vector> Element::* M = nullptr; + + if (cmd == "coefs") M = &Element::coefs; + if (cmd == "coefs2") M = &Element::coefs2; + if (cmd == "coefs3") M = &Element::coefs3; + if (cmd == "coefs4") M = &Element::coefs4; + + std::vector row; + std::vector col; + std::vector val; + + unsigned int max_row = 0; + unsigned int max_col = 0; + + while (!ss.eof()) { + unsigned int r, c; + double v; + ss >> r >> c >> v; + + if (ss.eof()) break; + + row.push_back(r); + col.push_back(c); + val.push_back(v); + + if (r + 1 > max_row) max_row = r + 1; + if (c + 1 > max_col) max_col = c + 1; + } + + if (max_row > 0 && max_col > 0) { + (e.*M).assign(max_row, std::vector(max_col, 0.0)); + + for (unsigned int k = 0; k < val.size(); ++k) + (e.*M)[row[k]][col[k]] = val[k]; + } + + continue; + } if (line.size()<2) continue; return Status::flat_file_error; } @@ -490,6 +547,13 @@ static bool has_matrix66(const Matrix& m) { return false; } +static bool has_coeffs(const std::vector>& matrix) { + if (!matrix.empty()) + return true; + else + return false; +} + static bool has_polynom(const std::vector& p) { for (int i=0; i>& M) +{ + fp << std::setw(pw) << label; + + for (unsigned int i = 0; i < M.size(); ++i) { + for (unsigned int j = 0; j < M[i].size(); ++j) { + + fp.unsetf(std::ios_base::showpos); + fp << i << ' ' << j << ' '; + fp.setf(std::ios_base::showpos); + fp << M[i][j] << ' '; + + } + } + fp << '\n'; +} + static void write_polynom(std::ostream& fp, const std::string& label, const std::vector& p) { fp << std::setw(pw) << label; for (int i=0; i Date: Tue, 18 Nov 2025 16:29:07 -0300 Subject: [PATCH 19/24] add prop_step to id_passmethod --- include/trackcpp/field3d.hpp | 16 ++++++++-------- include/trackcpp/passmethods.hpp | 12 ------------ 2 files changed, 8 insertions(+), 20 deletions(-) diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 23a6015..379ea6a 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -253,15 +253,15 @@ void prop_step(const double& beta0, const double& brho, const double& kx, const const std::vector>& coefs, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, const T& pnorm, double& s, double step) { - // prop_h1(beta0, map, pnorm, s, step); - // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_h1(beta0, map, pnorm, s, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); prop_h2(beta0, map, pnorm, step); - // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - // prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); prop_h3(beta0, map, pnorm, step); - // prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); prop_h2(beta0, map, pnorm, step); - // prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - // prop_h1(beta0, map, pnorm, s, step); + prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); + prop_h1(beta0, map, pnorm, s, step); } diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 6ab4f5b..9895f1c 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -477,23 +477,11 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, global_2_local(pos, elem); const double brho = get_magnetic_rigidity(accelerator.energy); const double gamma = accelerator.energy / electron_rest_energy_eV; - // const double beta0 = sqrt(1 - 1/(gamma*gamma)); const double beta0 = 1.0; double step = elem.length / float(elem.nr_steps); double s0 = elem.s0; - // T d = calc_D(beta0, pos.de); T pnorm = 1 / (1 + pos.de); for (int i=0; i Date: Mon, 8 Jun 2026 15:44:31 -0300 Subject: [PATCH 20/24] change coefs to coefs1 --- include/trackcpp/elements.h | 6 +-- include/trackcpp/field3d.h | 20 +++---- include/trackcpp/field3d.hpp | 91 ++++++++++++++------------------ include/trackcpp/passmethods.hpp | 2 +- src/elements.cpp | 8 +-- src/flat_file.cpp | 6 +-- 6 files changed, 61 insertions(+), 72 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index d657b30..61f06e4 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -66,7 +66,7 @@ class Element { double ks = 0; // [1/m] double kx = 0; // [1/m] double s0 = 0; // [m] - std::vector> coefs; + std::vector> coefs1; std::vector> coefs2; std::vector> coefs3; std::vector> coefs4; @@ -119,7 +119,7 @@ class Element { static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); static Element field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs1_, const std::vector>& coefs2_, const std::vector>& coefs3_, const std::vector>& coefs4_, const int nr_steps_ = 40); @@ -144,7 +144,7 @@ void initialize_sextupole(Element& element, const double& S, const int& nr_steps void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); void initialize_field3d(Element& element, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs1_, const std::vector>& coefs2_, const std::vector>& coefs3_, const std::vector>& coefs4_, const int nr_steps_); diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h index 8fee55d..5f83802 100644 --- a/include/trackcpp/field3d.h +++ b/include/trackcpp/field3d.h @@ -23,19 +23,19 @@ #include "pos.h" template -T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); +T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); template -T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); +T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); template -T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); +T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); template -T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, const T& x, const T& y, const double& s); +T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); template @@ -48,11 +48,11 @@ void exp_h1_z(const double& beta0, Pos& map, double step); void exp_h1_s(double& s, double step); template -void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); +void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); template -void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); +void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); template void exp_h2_y(const double& beta0, Pos& map, double step); @@ -63,11 +63,11 @@ void exp_h2_z(const double& beta0, Pos& map, double step); template -void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); +void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); template -void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); +void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); template @@ -86,9 +86,9 @@ template void prop_h3(const double& beta0, Pos& map, double step); template -void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); +void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); template -void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs, Pos& map, double s, int sign, double step); +void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); #endif diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 379ea6a..983592d 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -22,18 +22,18 @@ template T ay(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, const T& x, const T& y, const double& s) { T ay_ = 0.0; - int M = coefs.size(); - int N = coefs[0].size(); + int M = coefs1.size(); + int N = coefs1[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs[m - 1][n - 1] * (m * kx) / (n * ks * ky); + double fac1 = coefs1[m - 1][n - 1] * (m * kx) / (n * ks * ky); double fac2 = coefs2[m - 1][n - 1] * (m * kx) / (n * ks * ky); double fac3 = -coefs3[m - 1][n - 1] * (m * kx) / (n * ks * ky); double fac4 = -coefs4[m - 1][n - 1] * (m * kx) / (n * ks * ky); @@ -49,18 +49,18 @@ T ay(const double& brho, const double& kx, const double& ks, template T ax(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, const T& x, const T& y, const double& s) { T ax_ = 0.0; - const int M = coefs.size(); - const int N = coefs[0].size(); + const int M = coefs1.size(); + const int N = coefs1[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs[m - 1][n - 1] / (n * ks); + double fac1 = coefs1[m - 1][n - 1] / (n * ks); double fac2 = coefs2[m - 1][n - 1] / (n * ks); double fac3 = -coefs3[m - 1][n - 1] / (n * ks); double fac4 = -coefs4[m - 1][n - 1] / (n * ks); @@ -76,18 +76,18 @@ T ax(const double& brho, const double& kx, const double& ks, template T inty_day_dx(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, const T& x, const T& y, const double& s) { T day_dx = 0.0; - int M = coefs.size(); - int N = coefs[0].size(); + int M = coefs1.size(); + int N = coefs1[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); + double fac1 = coefs1[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); double fac2 = -coefs2[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); double fac3 = -coefs3[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); double fac4 = coefs4[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); @@ -103,18 +103,18 @@ T inty_day_dx(const double& brho, const double& kx, const double& ks, template T intx_dax_dy(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, const T& x, const T& y, const double& s) { T dax_dy = 0.0; - int M = coefs.size(); - int N = coefs[0].size(); + int M = coefs1.size(); + int N = coefs1[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs[m - 1][n - 1] * ky / (n * ks * m * kx); + double fac1 = coefs1[m - 1][n - 1] * ky / (n * ks * m * kx); double fac2 = -coefs2[m - 1][n - 1] * ky / (n * ks * m * kx); double fac3 = -coefs3[m - 1][n - 1] * ky / (n * ks * m * kx); double fac4 = coefs4[m - 1][n - 1] * ky / (n * ks * m * kx); @@ -128,12 +128,6 @@ T intx_dax_dy(const double& brho, const double& kx, const double& ks, return 1* dax_dy/brho; } -template -T calc_D(const double& beta0, const T& delta) { - // return sqrt(1.0 + 2.0 * delta / beta0 + delta * delta); - return 1.0 + delta; -} - template void exp_h1_z(const double& beta0, Pos& map, const T& d, double step) { T factor = (1.0 / beta0 - (1.0 / beta0 + map.de) / d); @@ -147,20 +141,20 @@ void exp_h1_s(T& s, T step) { template void exp_iy_px(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = inty_day_dx(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = inty_day_dx(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } template void exp_iy_py(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = ay(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = ay(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } @@ -172,28 +166,26 @@ void exp_h2_y(const double& beta0, Pos& map, const T& pnorm, double step) { template void exp_h2_z(const double& beta0, Pos& map, const T& pnorm, double step) { - // T factor = (map.py * map.py)*(1/beta0+map.de)/(2*d*d*d); - // T factor = (map.py * map.py)/d/d*0.5; map.dl += 0.25*step*pnorm*pnorm*(map.py * map.py); } template void exp_ix_px(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = ax(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = ax(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } template void exp_ix_py(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - T factor = intx_dax_dy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = intx_dax_dy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } @@ -205,14 +197,11 @@ void exp_h3_x(const double& beta0, Pos& map, const T& pnorm, double step) { template void exp_h3_z(const double& beta0, Pos& map, const T& pnorm, double step) { - // T factor = (map.px*map.px)*(1/beta0+map.de)/(2*d*d*d); - // T factor = (map.px*map.px)/d/d*0.5; map.dl += 0.5*step*pnorm*pnorm*(map.px * map.px); } template -void prop_h1(const double& beta0, Pos& map, const T& d, double& s, double step) { - // exp_h1_z(beta0, map, d, step); +void prop_h1(Pos& map, double& s, double step) { exp_h1_s(s, step); } @@ -230,38 +219,38 @@ void prop_h3(const double& beta0, Pos& map, const T& pnorm, double step) { template void prop_ix(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - exp_ix_px(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); - exp_ix_py(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); + exp_ix_px(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); + exp_ix_py(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); } template void prop_iy(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, double s, int sign, double step) { - exp_iy_px(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); - exp_iy_py(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, sign, step); + exp_iy_px(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); + exp_iy_py(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); } template void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, Pos& map, const T& pnorm, double& s, double step) { - prop_h1(beta0, map, pnorm, s, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_h1(map, s, step); + prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, +1, step); prop_h2(beta0, map, pnorm, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, -1, step); + prop_ix(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, +1, step); prop_h3(beta0, map, pnorm, step); - prop_ix(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, +1, step); + prop_ix(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, -1, step); + prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, +1, step); prop_h2(beta0, map, pnorm, step); - prop_iy(brho, kx, ks, coefs, coefs2, coefs3, coefs4, map, s, -1, step); - prop_h1(beta0, map, pnorm, s, step); + prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, -1, step); + prop_h1(map, s, step); } diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 9895f1c..095a50e 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -482,7 +482,7 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, double s0 = elem.s0; T pnorm = 1 / (1 + pos.de); for (int i=0; i>& coefs_, const std::vector>& coefs2_, + const std::vector>& coefs1_, const std::vector>& coefs2_, const std::vector>& coefs3_, const std::vector>& coefs4_, const int nr_steps_) { Element e = Element(fam_name_, length_); - initialize_field3d(e, s0_, kx_, ks_, coefs_, coefs2_, coefs3_, coefs4_,nr_steps_); + initialize_field3d(e, s0_, kx_, ks_, coefs1_, coefs2_, coefs3_, coefs4_, nr_steps_); return e; } @@ -319,12 +319,12 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n } void initialize_field3d(Element& element, const double& s0, const double& kx, const double& ks, - const std::vector>& coefs, const std::vector>& coefs2, + const std::vector>& coefs1, const std::vector>& coefs2, const std::vector>& coefs3, const std::vector>& coefs4, const int nr_steps) { element.pass_method = PassMethod::pm_field3d_pass; element.nr_steps = nr_steps; - element.coefs = coefs; + element.coefs1 = coefs1; element.coefs2 = coefs2; element.coefs3 = coefs3; element.coefs4 = coefs4; diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 81a0ef0..abe9ef3 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -132,7 +132,7 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) if (e.kx != 0) { fp << std::setw(pw) << "kx" << e.kx << '\n'; } if (e.ks != 0) { fp << std::setw(pw) << "ks" << e.ks << '\n'; } if (e.s0 != 0) { fp << std::setw(pw) << "s0" << e.s0 << '\n'; } - if (has_coeffs(e.coefs)) write_coeffs(fp, "coefs", e.coefs); + if (has_coeffs(e.coefs1)) write_coeffs(fp, "coefs1", e.coefs1); if (has_coeffs(e.coefs2)) write_coeffs(fp, "coefs2", e.coefs2); if (has_coeffs(e.coefs3)) write_coeffs(fp, "coefs3", e.coefs3); if (has_coeffs(e.coefs4)) write_coeffs(fp, "coefs4", e.coefs4); @@ -318,13 +318,13 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) continue; } - if (cmd.compare("coefs") == 0 || + if (cmd.compare("coefs1") == 0 || cmd.compare("coefs2") == 0 || cmd.compare("coefs3") == 0 || cmd.compare("coefs4") == 0) { std::vector> Element::* M = nullptr; - if (cmd == "coefs") M = &Element::coefs; + if (cmd == "coefs1") M = &Element::coefs1; if (cmd == "coefs2") M = &Element::coefs2; if (cmd == "coefs3") M = &Element::coefs3; if (cmd == "coefs4") M = &Element::coefs4; From 386fd270a01af64bec1c4030986106043e754b99 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 9 Jun 2026 09:53:46 -0300 Subject: [PATCH 21/24] change trackcpp structure to work only with 2 coefs set in symplectic integration of 3d fields --- include/trackcpp/elements.h | 10 +-- include/trackcpp/field3d.h | 2 +- include/trackcpp/field3d.hpp | 134 +++++++++++++------------------ include/trackcpp/passmethods.hpp | 3 +- python_package/interface.cpp | 6 +- python_package/interface.h | 4 +- src/elements.cpp | 12 ++- src/flat_file.cpp | 8 +- 8 files changed, 73 insertions(+), 106 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 61f06e4..892df2d 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -68,8 +68,6 @@ class Element { double s0 = 0; // [m] std::vector> coefs1; std::vector> coefs2; - std::vector> coefs3; - std::vector> coefs4; std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; @@ -119,8 +117,8 @@ class Element { static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); static Element field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs1_, const std::vector>& coefs2_, - const std::vector>& coefs3_, const std::vector>& coefs4_, + const std::vector>& coefs1_, + const std::vector>& coefs2_, const int nr_steps_ = 40); bool operator==(const Element& o) const; @@ -144,8 +142,8 @@ void initialize_sextupole(Element& element, const double& S, const int& nr_steps void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); void initialize_field3d(Element& element, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs1_, const std::vector>& coefs2_, - const std::vector>& coefs3_, const std::vector>& coefs4_, + const std::vector>& coefs1_, + const std::vector>& coefs2_, const int nr_steps_); #endif diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h index 5f83802..4556665 100644 --- a/include/trackcpp/field3d.h +++ b/include/trackcpp/field3d.h @@ -77,7 +77,7 @@ template void exp_h3_z(const double& beta0, Pos& map, double step); template -void prop_h1(const double& beta0, Pos& map, double& s, double step); +void prop_h1(Pos& map, double& s, double step); template void prop_h2(const double& beta0, Pos& map, double step); diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 983592d..c789566 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -22,8 +22,8 @@ template T ay(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, const T& x, const T& y, const double& s) { T ay_ = 0.0; @@ -34,13 +34,9 @@ T ay(const double& brho, const double& kx, const double& ks, for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); double fac1 = coefs1[m - 1][n - 1] * (m * kx) / (n * ks * ky); - double fac2 = coefs2[m - 1][n - 1] * (m * kx) / (n * ks * ky); - double fac3 = -coefs3[m - 1][n - 1] * (m * kx) / (n * ks * ky); - double fac4 = -coefs4[m - 1][n - 1] * (m * kx) / (n * ks * ky); + double fac2 = -coefs2[m - 1][n - 1] * (m * kx) / (n * ks * ky); ay_ += fac1 * sin(m * kx * x) * sinh(ky * y) * std::sin(n * ks * s); - ay_ += fac2 * cos(m * kx * x) * sinh(ky * y) * std::sin(n * ks * s); - ay_ += fac3 * sin(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); - ay_ += fac4 * cos(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); + ay_ += fac2 * sin(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); } } @@ -49,8 +45,8 @@ T ay(const double& brho, const double& kx, const double& ks, template T ax(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, const T& x, const T& y, const double& s) { T ax_ = 0.0; @@ -61,13 +57,9 @@ T ax(const double& brho, const double& kx, const double& ks, for (int n = 1; n <= N; ++n) { double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); double fac1 = coefs1[m - 1][n - 1] / (n * ks); - double fac2 = coefs2[m - 1][n - 1] / (n * ks); - double fac3 = -coefs3[m - 1][n - 1] / (n * ks); - double fac4 = -coefs4[m - 1][n - 1] / (n * ks); + double fac2 = -coefs2[m - 1][n - 1] / (n * ks); ax_ += fac1 * cos(m * kx * x) * cosh(ky * y) * std::sin(n * ks * s); - ax_ += fac2 * sin(m * kx * x) * cosh(ky * y) * std::sin(n * ks * s); - ax_ += fac3 * cos(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); - ax_ += fac4 * sin(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); + ax_ += fac2 * cos(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); } } @@ -76,8 +68,8 @@ T ax(const double& brho, const double& kx, const double& ks, template T inty_day_dx(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, const T& x, const T& y, const double& s) { T day_dx = 0.0; @@ -89,12 +81,8 @@ T inty_day_dx(const double& brho, const double& kx, const double& ks, double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); double fac1 = coefs1[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); double fac2 = -coefs2[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); - double fac3 = -coefs3[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); - double fac4 = coefs4[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); day_dx += fac1 * cos(m * kx * x) * std::sin(n * ks * s) * (cosh(ky * y) - 1.0); - day_dx += fac2 * sin(m * kx * x) * std::sin(n * ks * s) * (cosh(ky * y) - 1.0); - day_dx += fac3 * cos(m * kx * x) * std::cos(n * ks * s) * (cosh(ky * y) - 1.0); - day_dx += fac4 * sin(m * kx * x) * std::cos(n * ks * s) * (cosh(ky * y) - 1.0); + day_dx += fac2 * cos(m * kx * x) * std::cos(n * ks * s) * (cosh(ky * y) - 1.0); } } @@ -103,8 +91,8 @@ T inty_day_dx(const double& brho, const double& kx, const double& ks, template T intx_dax_dy(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, const T& x, const T& y, const double& s) { T dax_dy = 0.0; @@ -116,24 +104,14 @@ T intx_dax_dy(const double& brho, const double& kx, const double& ks, double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); double fac1 = coefs1[m - 1][n - 1] * ky / (n * ks * m * kx); double fac2 = -coefs2[m - 1][n - 1] * ky / (n * ks * m * kx); - double fac3 = -coefs3[m - 1][n - 1] * ky / (n * ks * m * kx); - double fac4 = coefs4[m - 1][n - 1] * ky / (n * ks * m * kx); dax_dy += fac1 * sin(m * kx * x) * std::sin(n * ks * s) * sinh(ky * y); - dax_dy += fac2 * cos(m * kx * x) * std::sin(n * ks * s) * sinh(ky * y); - dax_dy += fac3 * sin(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); - dax_dy += fac4 * cos(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); + dax_dy += fac2 * sin(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); } } return 1* dax_dy/brho; } -template -void exp_h1_z(const double& beta0, Pos& map, const T& d, double step) { - T factor = (1.0 / beta0 - (1.0 / beta0 + map.de) / d); - map.dl += factor * step / 2.0; -} - template void exp_h1_s(T& s, T step) { s += step / 2.0; @@ -141,62 +119,62 @@ void exp_h1_s(T& s, T step) { template void exp_iy_px(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, Pos& map, double s, int sign, double step) { - T factor = inty_day_dx(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = inty_day_dx(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } template void exp_iy_py(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, Pos& map, double s, int sign, double step) { - T factor = ay(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = ay(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } template -void exp_h2_y(const double& beta0, Pos& map, const T& pnorm, double step) { +void exp_h2_y(Pos& map, const T& pnorm, double step) { map.ry += 0.5*step*pnorm*map.py; } template -void exp_h2_z(const double& beta0, Pos& map, const T& pnorm, double step) { +void exp_h2_z(Pos& map, const T& pnorm, double step) { map.dl += 0.25*step*pnorm*pnorm*(map.py * map.py); } template void exp_ix_px(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, Pos& map, double s, int sign, double step) { - T factor = ax(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = ax(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } template void exp_ix_py(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, Pos& map, double s, int sign, double step) { - T factor = intx_dax_dy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map.rx, map.ry, s) * sign * -1.0; + T factor = intx_dax_dy(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } template -void exp_h3_x(const double& beta0, Pos& map, const T& pnorm, double step) { +void exp_h3_x(Pos& map, const T& pnorm, double step) { map.rx += step*pnorm*map.px; } template -void exp_h3_z(const double& beta0, Pos& map, const T& pnorm, double step) { +void exp_h3_z(Pos& map, const T& pnorm, double step) { map.dl += 0.5*step*pnorm*pnorm*(map.px * map.px); } @@ -206,51 +184,51 @@ void prop_h1(Pos& map, double& s, double step) { } template -void prop_h2(const double& beta0, Pos& map, const T& pnorm, double step) { - exp_h2_y(beta0, map, pnorm, step); - exp_h2_z(beta0, map, pnorm, step); +void prop_h2(Pos& map, const T& pnorm, double step) { + exp_h2_y(map, pnorm, step); + exp_h2_z(map, pnorm, step); } template -void prop_h3(const double& beta0, Pos& map, const T& pnorm, double step) { - exp_h3_x(beta0, map, pnorm, step); - exp_h3_z(beta0, map, pnorm, step); +void prop_h3(Pos& map, const T& pnorm, double step) { + exp_h3_x(map, pnorm, step); + exp_h3_z(map, pnorm, step); } template void prop_ix(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, Pos& map, double s, int sign, double step) { - exp_ix_px(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); - exp_ix_py(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); + exp_ix_px(brho, kx, ks, coefs1, coefs2, map, s, sign, step); + exp_ix_py(brho, kx, ks, coefs1, coefs2, map, s, sign, step); } template void prop_iy(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, Pos& map, double s, int sign, double step) { - exp_iy_px(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); - exp_iy_py(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, sign, step); + exp_iy_px(brho, kx, ks, coefs1, coefs2, map, s, sign, step); + exp_iy_py(brho, kx, ks, coefs1, coefs2, map, s, sign, step); } template -void prop_step(const double& beta0, const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, +void prop_step(const double& brho, const double& kx, const double& ks, + const std::vector>& coefs1, + const std::vector>& coefs2, Pos& map, const T& pnorm, double& s, double step) { prop_h1(map, s, step); - prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, +1, step); - prop_h2(beta0, map, pnorm, step); - prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, -1, step); - prop_ix(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, +1, step); - prop_h3(beta0, map, pnorm, step); - prop_ix(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, -1, step); - prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, +1, step); - prop_h2(beta0, map, pnorm, step); - prop_iy(brho, kx, ks, coefs1, coefs2, coefs3, coefs4, map, s, -1, step); + prop_iy(brho, kx, ks, coefs1, coefs2, map, s, +1, step); + prop_h2(map, pnorm, step); + prop_iy(brho, kx, ks, coefs1, coefs2, map, s, -1, step); + prop_ix(brho, kx, ks, coefs1, coefs2, map, s, +1, step); + prop_h3(map, pnorm, step); + prop_ix(brho, kx, ks, coefs1, coefs2, map, s, -1, step); + prop_iy(brho, kx, ks, coefs1, coefs2, map, s, +1, step); + prop_h2(map, pnorm, step); + prop_iy(brho, kx, ks, coefs1, coefs2, map, s, -1, step); prop_h1(map, s, step); } diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 095a50e..f41e32c 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -477,12 +477,11 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, global_2_local(pos, elem); const double brho = get_magnetic_rigidity(accelerator.energy); const double gamma = accelerator.energy / electron_rest_energy_eV; - const double beta0 = 1.0; double step = elem.length / float(elem.nr_steps); double s0 = elem.s0; T pnorm = 1 / (1 + pos.de); for (int i=0; i>& coefs_, const std::vector>& coefs2_, - const std::vector>& coefs3_, const std::vector>& coefs4_, + const std::vector>& coefs1_, + const std::vector>& coefs3_, const int nr_steps_) { - return Element::field3d(fam_name_, length_, s0_, kx_, ks_, coefs_, coefs2_, coefs3_, coefs4_, nr_steps_); + return Element::field3d(fam_name_, length_, s0_, kx_, ks_, coefs1_, coefs3_, nr_steps_); } Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag) { diff --git a/python_package/interface.h b/python_package/interface.h index 4c9aba4..7320cbe 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -92,8 +92,8 @@ Element sextupole_wrapper(const std::string& fam_name_, const double& length_, c Element rfcavity_wrapper(const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag); Element kickmap_wrapper(const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length = 1.0, const double& rescale_kicks = 1.0); Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs_, const std::vector>& coefs2_, - const std::vector>& coefs3_, const std::vector>& coefs4_, + const std::vector>& coefs1_, + const std::vector>& coefs2_, const int nr_steps_); Element rbend_wrapper(const std::string& fam_name_, const double& length_, const double& angle_, const double& angle_in_, const double& angle_out_, diff --git a/src/elements.cpp b/src/elements.cpp index eaf141f..854456f 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -157,11 +157,11 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt } Element Element::field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs1_, const std::vector>& coefs2_, - const std::vector>& coefs3_, const std::vector>& coefs4_, + const std::vector>& coefs1_, + const std::vector>& coefs2_, const int nr_steps_) { Element e = Element(fam_name_, length_); - initialize_field3d(e, s0_, kx_, ks_, coefs1_, coefs2_, coefs3_, coefs4_, nr_steps_); + initialize_field3d(e, s0_, kx_, ks_, coefs1_, coefs2_, nr_steps_); return e; } @@ -319,15 +319,13 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n } void initialize_field3d(Element& element, const double& s0, const double& kx, const double& ks, - const std::vector>& coefs1, const std::vector>& coefs2, - const std::vector>& coefs3, const std::vector>& coefs4, + const std::vector>& coefs1, + const std::vector>& coefs2, const int nr_steps) { element.pass_method = PassMethod::pm_field3d_pass; element.nr_steps = nr_steps; element.coefs1 = coefs1; element.coefs2 = coefs2; - element.coefs3 = coefs3; - element.coefs4 = coefs4; element.kx = kx; element.ks = ks; element.s0 = s0; diff --git a/src/flat_file.cpp b/src/flat_file.cpp index abe9ef3..0f71e4f 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -134,8 +134,6 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) if (e.s0 != 0) { fp << std::setw(pw) << "s0" << e.s0 << '\n'; } if (has_coeffs(e.coefs1)) write_coeffs(fp, "coefs1", e.coefs1); if (has_coeffs(e.coefs2)) write_coeffs(fp, "coefs2", e.coefs2); - if (has_coeffs(e.coefs3)) write_coeffs(fp, "coefs3", e.coefs3); - if (has_coeffs(e.coefs4)) write_coeffs(fp, "coefs4", e.coefs4); if (e.has_t_in) write_6d_vector(fp, "t_in", e.t_in); @@ -319,15 +317,11 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) } if (cmd.compare("coefs1") == 0 || - cmd.compare("coefs2") == 0 || - cmd.compare("coefs3") == 0 || - cmd.compare("coefs4") == 0) { + cmd.compare("coefs2") == 0) { std::vector> Element::* M = nullptr; if (cmd == "coefs1") M = &Element::coefs1; if (cmd == "coefs2") M = &Element::coefs2; - if (cmd == "coefs3") M = &Element::coefs3; - if (cmd == "coefs4") M = &Element::coefs4; std::vector row; std::vector col; From a479d81c1c87eeae223ae6eeac395a98d96546fa Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 23 Jun 2026 17:06:06 -0300 Subject: [PATCH 22/24] change field3d properties names --- include/trackcpp/elements.h | 22 ++--- include/trackcpp/field3d.h | 60 ++++++++++--- include/trackcpp/field3d.hpp | 150 +++++++++++++++---------------- include/trackcpp/passmethods.hpp | 4 +- python_package/interface.cpp | 9 +- python_package/interface.h | 6 +- src/elements.cpp | 24 ++--- src/flat_file.cpp | 24 ++--- 8 files changed, 170 insertions(+), 129 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 892df2d..62cd1c4 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -63,11 +63,11 @@ class Element { double phase_lag = 0; // [rad] int kicktable_idx = -1; // index of kickmap object in kicktable_list double rescale_kicks = 1.0; // for kickmaps - double ks = 0; // [1/m] - double kx = 0; // [1/m] - double s0 = 0; // [m] - std::vector> coefs1; - std::vector> coefs2; + double field3d_hori_ks = 0; // [1/m] + double field3d_hori_kx = 0; // [1/m] + double field3d_hori_s0 = 0; // [m] + std::vector> field3d_hori_coefs_cos; + std::vector> field3d_hori_coefs_sin; std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; @@ -116,9 +116,9 @@ class Element { static Element sextupole (const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); - static Element field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs1_, - const std::vector>& coefs2_, + static Element field3d (const std::string& fam_name_, const double& length_, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_, + const std::vector>& field3d_hori_coefs_cos_, + const std::vector>& field3d_hori_coefs_sin_, const int nr_steps_ = 40); bool operator==(const Element& o) const; @@ -141,9 +141,9 @@ void initialize_quadrupole(Element& element, const double& K, const int& nr_step void initialize_sextupole(Element& element, const double& S, const int& nr_steps); void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); -void initialize_field3d(Element& element, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs1_, - const std::vector>& coefs2_, +void initialize_field3d(Element& element, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_, + const std::vector>& field3d_hori_coefs_cos_, + const std::vector>& field3d_hori_coefs_sin_, const int nr_steps_); #endif diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h index 4556665..90cbc5f 100644 --- a/include/trackcpp/field3d.h +++ b/include/trackcpp/field3d.h @@ -23,19 +23,35 @@ #include "pos.h" template -T ay(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); +T ay(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + const T& x, const T& y, const double& s); template -T ax(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); +T ax(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + const T& x, const T& y, const double& s); template -T inty_day_dx(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); +T inty_day_dx(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + const T& x, const T& y, const double& s); template -T intx_dax_dy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, const T& x, const T& y, const double& s); +T intx_dax_dy(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + const T& x, const T& y, const double& s); template @@ -48,11 +64,19 @@ void exp_h1_z(const double& beta0, Pos& map, double step); void exp_h1_s(double& s, double step); template -void exp_iy_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); +void exp_iy_px(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + Pos& map, double s, int sign, double step); template -void exp_iy_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); +void exp_iy_py(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + Pos& map, double s, int sign, double step); template void exp_h2_y(const double& beta0, Pos& map, double step); @@ -63,11 +87,19 @@ void exp_h2_z(const double& beta0, Pos& map, double step); template -void exp_ix_px(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); +void exp_ix_px(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + Pos& map, double s, int sign, double step); template -void exp_ix_py(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); +void exp_ix_py(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + Pos& map, double s, int sign, double step); template @@ -86,9 +118,17 @@ template void prop_h3(const double& beta0, Pos& map, double step); template -void prop_ix(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); +void prop_ix(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + Pos& map, double s, int sign, double step); template -void prop_iy(const double& brho, const double& kx, const double& ks, const std::vector>& coefs1, Pos& map, double s, int sign, double step); +void prop_iy(const double& brho, + const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, + Pos& map, double s, int sign, double step); #endif diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index c789566..198a82a 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -21,22 +21,22 @@ template -T ay(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +T ay(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, const T& x, const T& y, const double& s) { T ay_ = 0.0; - int M = coefs1.size(); - int N = coefs1[0].size(); + int M = hori_coefs_cos.size(); + int N = hori_coefs_cos[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { - double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs1[m - 1][n - 1] * (m * kx) / (n * ks * ky); - double fac2 = -coefs2[m - 1][n - 1] * (m * kx) / (n * ks * ky); - ay_ += fac1 * sin(m * kx * x) * sinh(ky * y) * std::sin(n * ks * s); - ay_ += fac2 * sin(m * kx * x) * sinh(ky * y) * std::cos(n * ks * s); + double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); + double fac1 = hori_coefs_cos[m - 1][n - 1] * (m * hori_kx) / (n * hori_ks * hori_ky); + double fac2 = -hori_coefs_sin[m - 1][n - 1] * (m * hori_kx) / (n * hori_ks * hori_ky); + ay_ += fac1 * sin(m * hori_kx * x) * sinh(hori_ky * y) * std::sin(n * hori_ks * s); + ay_ += fac2 * sin(m * hori_kx * x) * sinh(hori_ky * y) * std::cos(n * hori_ks * s); } } @@ -44,22 +44,22 @@ T ay(const double& brho, const double& kx, const double& ks, } template -T ax(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +T ax(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, const T& x, const T& y, const double& s) { T ax_ = 0.0; - const int M = coefs1.size(); - const int N = coefs1[0].size(); + const int M = hori_coefs_cos.size(); + const int N = hori_coefs_cos[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { - double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs1[m - 1][n - 1] / (n * ks); - double fac2 = -coefs2[m - 1][n - 1] / (n * ks); - ax_ += fac1 * cos(m * kx * x) * cosh(ky * y) * std::sin(n * ks * s); - ax_ += fac2 * cos(m * kx * x) * cosh(ky * y) * std::cos(n * ks * s); + double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); + double fac1 = hori_coefs_cos[m - 1][n - 1] / (n * hori_ks); + double fac2 = -hori_coefs_sin[m - 1][n - 1] / (n * hori_ks); + ax_ += fac1 * cos(m * hori_kx * x) * cosh(hori_ky * y) * std::sin(n * hori_ks * s); + ax_ += fac2 * cos(m * hori_kx * x) * cosh(hori_ky * y) * std::cos(n * hori_ks * s); } } @@ -67,22 +67,22 @@ T ax(const double& brho, const double& kx, const double& ks, } template -T inty_day_dx(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +T inty_day_dx(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, const T& x, const T& y, const double& s) { T day_dx = 0.0; - int M = coefs1.size(); - int N = coefs1[0].size(); + int M = hori_coefs_cos.size(); + int N = hori_coefs_cos[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { - double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs1[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); - double fac2 = -coefs2[m - 1][n - 1] * std::pow(m * kx, 2) / (n * ks * std::pow(ky, 2)); - day_dx += fac1 * cos(m * kx * x) * std::sin(n * ks * s) * (cosh(ky * y) - 1.0); - day_dx += fac2 * cos(m * kx * x) * std::cos(n * ks * s) * (cosh(ky * y) - 1.0); + double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); + double fac1 = hori_coefs_cos[m - 1][n - 1] * std::pow(m * hori_kx, 2) / (n * hori_ks * std::pow(hori_ky, 2)); + double fac2 = -hori_coefs_sin[m - 1][n - 1] * std::pow(m * hori_kx, 2) / (n * hori_ks * std::pow(hori_ky, 2)); + day_dx += fac1 * cos(m * hori_kx * x) * std::sin(n * hori_ks * s) * (cosh(hori_ky * y) - 1.0); + day_dx += fac2 * cos(m * hori_kx * x) * std::cos(n * hori_ks * s) * (cosh(hori_ky * y) - 1.0); } } @@ -90,22 +90,22 @@ T inty_day_dx(const double& brho, const double& kx, const double& ks, } template -T intx_dax_dy(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +T intx_dax_dy(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, const T& x, const T& y, const double& s) { T dax_dy = 0.0; - int M = coefs1.size(); - int N = coefs1[0].size(); + int M = hori_coefs_cos.size(); + int N = hori_coefs_cos[0].size(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { - double ky = std::sqrt(std::pow(m * kx, 2) + std::pow(n * ks, 2)); - double fac1 = coefs1[m - 1][n - 1] * ky / (n * ks * m * kx); - double fac2 = -coefs2[m - 1][n - 1] * ky / (n * ks * m * kx); - dax_dy += fac1 * sin(m * kx * x) * std::sin(n * ks * s) * sinh(ky * y); - dax_dy += fac2 * sin(m * kx * x) * std::cos(n * ks * s) * sinh(ky * y); + double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); + double fac1 = hori_coefs_cos[m - 1][n - 1] * hori_ky / (n * hori_ks * m * hori_kx); + double fac2 = -hori_coefs_sin[m - 1][n - 1] * hori_ky / (n * hori_ks * m * hori_kx); + dax_dy += fac1 * sin(m * hori_kx * x) * std::sin(n * hori_ks * s) * sinh(hori_ky * y); + dax_dy += fac2 * sin(m * hori_kx * x) * std::cos(n * hori_ks * s) * sinh(hori_ky * y); } } @@ -118,21 +118,21 @@ void exp_h1_s(T& s, T step) { } template -void exp_iy_px(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void exp_iy_px(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, Pos& map, double s, int sign, double step) { - T factor = inty_day_dx(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; + T factor = inty_day_dx(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } template -void exp_iy_py(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void exp_iy_py(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, Pos& map, double s, int sign, double step) { - T factor = ay(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; + T factor = ay(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } @@ -149,21 +149,21 @@ void exp_h2_z(Pos& map, const T& pnorm, double step) { template -void exp_ix_px(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void exp_ix_px(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, Pos& map, double s, int sign, double step) { - T factor = ax(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; + T factor = ax(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.px += factor; } template -void exp_ix_py(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void exp_ix_py(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, Pos& map, double s, int sign, double step) { - T factor = intx_dax_dy(brho, kx, ks, coefs1, coefs2, map.rx, map.ry, s) * sign * -1.0; + T factor = intx_dax_dy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.py += factor; } @@ -196,39 +196,39 @@ void prop_h3(Pos& map, const T& pnorm, double step) { } template -void prop_ix(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void prop_ix(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, Pos& map, double s, int sign, double step) { - exp_ix_px(brho, kx, ks, coefs1, coefs2, map, s, sign, step); - exp_ix_py(brho, kx, ks, coefs1, coefs2, map, s, sign, step); + exp_ix_px(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); + exp_ix_py(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); } template -void prop_iy(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void prop_iy(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, Pos& map, double s, int sign, double step) { - exp_iy_px(brho, kx, ks, coefs1, coefs2, map, s, sign, step); - exp_iy_py(brho, kx, ks, coefs1, coefs2, map, s, sign, step); + exp_iy_px(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); + exp_iy_py(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); } template -void prop_step(const double& brho, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void prop_step(const double& brho, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, Pos& map, const T& pnorm, double& s, double step) { prop_h1(map, s, step); - prop_iy(brho, kx, ks, coefs1, coefs2, map, s, +1, step); + prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, +1, step); prop_h2(map, pnorm, step); - prop_iy(brho, kx, ks, coefs1, coefs2, map, s, -1, step); - prop_ix(brho, kx, ks, coefs1, coefs2, map, s, +1, step); + prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, -1, step); + prop_ix(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, +1, step); prop_h3(map, pnorm, step); - prop_ix(brho, kx, ks, coefs1, coefs2, map, s, -1, step); - prop_iy(brho, kx, ks, coefs1, coefs2, map, s, +1, step); + prop_ix(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, -1, step); + prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, +1, step); prop_h2(map, pnorm, step); - prop_iy(brho, kx, ks, coefs1, coefs2, map, s, -1, step); + prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, -1, step); prop_h1(map, s, step); } diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 09023a6..9baccdb 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -479,10 +479,10 @@ Status::type pm_field3d_pass(Pos &pos, const Element &elem, const double brho = get_magnetic_rigidity(accelerator.energy); const double gamma = accelerator.energy / electron_rest_energy_eV; double step = elem.length / float(elem.nr_steps); - double s0 = elem.s0; + double s0 = elem.field3d_hori_s0; T pnorm = 1 / (1 + pos.de); for (int i=0; i>& coefs1_, - const std::vector>& coefs3_, +Element field3d_wrapper(const std::string& fam_name_, const double& length_, + const double& hori_s0_, const double& hori_kx_, const double& hori_ks_, + const std::vector>& hori_coefs_cos_, + const std::vector>& hori_coefs_sin_, const int nr_steps_) { - return Element::field3d(fam_name_, length_, s0_, kx_, ks_, coefs1_, coefs3_, nr_steps_); + return Element::field3d(fam_name_, length_, hori_s0_, hori_kx_, hori_ks_, hori_coefs_cos_, hori_coefs_sin_, nr_steps_); } Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag) { diff --git a/python_package/interface.h b/python_package/interface.h index 7320cbe..f245f39 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -91,9 +91,9 @@ Element quadrupole_wrapper(const std::string& fam_name_, const double& length_, Element sextupole_wrapper(const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5); Element rfcavity_wrapper(const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag); Element kickmap_wrapper(const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length = 1.0, const double& rescale_kicks = 1.0); -Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs1_, - const std::vector>& coefs2_, +Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& field3d_hori_s0_, const double& field3d_hori_kx_, + const double& field3d_hori_ks_, const std::vector>& field3d_hori_coefs_cos_, + const std::vector>& field3d_hori_coefs_sin_, const int nr_steps_); Element rbend_wrapper(const std::string& fam_name_, const double& length_, const double& angle_, const double& angle_in_, const double& angle_out_, diff --git a/src/elements.cpp b/src/elements.cpp index 854456f..7c853f4 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -156,12 +156,12 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt return e; } -Element Element::field3d (const std::string& fam_name_, const double& length_, const double& s0_, const double& kx_, const double& ks_, - const std::vector>& coefs1_, - const std::vector>& coefs2_, +Element Element::field3d (const std::string& fam_name_, const double& length_, const double& hori_s0_, const double& hori_kx_, const double& hori_ks_, + const std::vector>& hori_coefs_cos_, + const std::vector>& hori_coefs_sin_, const int nr_steps_) { Element e = Element(fam_name_, length_); - initialize_field3d(e, s0_, kx_, ks_, coefs1_, coefs2_, nr_steps_); + initialize_field3d(e, hori_s0_, hori_kx_, hori_ks_, hori_coefs_cos_, hori_coefs_sin_, nr_steps_); return e; } @@ -318,15 +318,15 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n element.rescale_kicks = rescale_kicks; } -void initialize_field3d(Element& element, const double& s0, const double& kx, const double& ks, - const std::vector>& coefs1, - const std::vector>& coefs2, +void initialize_field3d(Element& element, const double& hori_s0, const double& hori_kx, const double& hori_ks, + const std::vector>& hori_coefs_cos, + const std::vector>& hori_coefs_sin, const int nr_steps) { element.pass_method = PassMethod::pm_field3d_pass; element.nr_steps = nr_steps; - element.coefs1 = coefs1; - element.coefs2 = coefs2; - element.kx = kx; - element.ks = ks; - element.s0 = s0; + element.field3d_hori_coefs_cos = hori_coefs_cos; + element.field3d_hori_coefs_sin = hori_coefs_sin; + element.field3d_hori_kx = hori_kx; + element.field3d_hori_ks = hori_ks; + element.field3d_hori_s0 = hori_s0; } diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 0f71e4f..3990154 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -129,11 +129,11 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) if (e.angle_in != 0) { fp << std::setw(pw) << "angle_in" << e.angle_in << '\n'; } if (e.angle_out != 0) { fp << std::setw(pw) << "angle_out" << e.angle_out << '\n'; } if (e.rescale_kicks != 1.0) { fp << std::setw(pw) << "rescale_kicks" << e.rescale_kicks << '\n'; } - if (e.kx != 0) { fp << std::setw(pw) << "kx" << e.kx << '\n'; } - if (e.ks != 0) { fp << std::setw(pw) << "ks" << e.ks << '\n'; } - if (e.s0 != 0) { fp << std::setw(pw) << "s0" << e.s0 << '\n'; } - if (has_coeffs(e.coefs1)) write_coeffs(fp, "coefs1", e.coefs1); - if (has_coeffs(e.coefs2)) write_coeffs(fp, "coefs2", e.coefs2); + if (e.field3d_hori_kx != 0) { fp << std::setw(pw) << "kx" << e.field3d_hori_kx << '\n'; } + if (e.field3d_hori_ks != 0) { fp << std::setw(pw) << "ks" << e.field3d_hori_ks << '\n'; } + if (e.field3d_hori_s0 != 0) { fp << std::setw(pw) << "s0" << e.field3d_hori_s0 << '\n'; } + if (has_coeffs(e.field3d_hori_coefs_cos)) write_coeffs(fp, "coefs_cos", e.field3d_hori_coefs_cos); + if (has_coeffs(e.field3d_hori_coefs_sin)) write_coeffs(fp, "coefs_sin", e.field3d_hori_coefs_sin); if (e.has_t_in) write_6d_vector(fp, "t_in", e.t_in); @@ -235,9 +235,9 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) if (cmd.compare("angle_in") == 0) { ss >> e.angle_in; continue; } if (cmd.compare("angle_out") == 0) { ss >> e.angle_out; continue; } if (cmd.compare("rescale_kicks") == 0) { ss >> e.rescale_kicks; continue; } - if (cmd.compare("kx") == 0) { ss >> e.kx; continue; } - if (cmd.compare("ks") == 0) { ss >> e.ks; continue; } - if (cmd.compare("s0") == 0) { ss >> e.s0; continue; } + if (cmd.compare("kx") == 0) { ss >> e.field3d_hori_kx; continue; } + if (cmd.compare("ks") == 0) { ss >> e.field3d_hori_ks; continue; } + if (cmd.compare("s0") == 0) { ss >> e.field3d_hori_s0; continue; } if (cmd.compare("t_in") == 0) { for(auto i=0; i<6; ++i) ss >> e.t_in[i]; e.reflag_t_in(); continue; } if (cmd.compare("t_out") == 0) { for(auto i=0; i<6; ++i) ss >> e.t_out[i]; e.reflag_t_out(); continue; } if (cmd.compare("rx|r_in") == 0) { for(auto i=0; i<6; ++i) ss >> e.r_in[0*6+i]; e.reflag_r_in(); continue; } @@ -316,12 +316,12 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) continue; } - if (cmd.compare("coefs1") == 0 || - cmd.compare("coefs2") == 0) { + if (cmd.compare("coefs_cos") == 0 || + cmd.compare("coefs_sin") == 0) { std::vector> Element::* M = nullptr; - if (cmd == "coefs1") M = &Element::coefs1; - if (cmd == "coefs2") M = &Element::coefs2; + if (cmd == "coefs_cos") M = &Element::field3d_hori_coefs_cos; + if (cmd == "coefs_sin") M = &Element::field3d_hori_coefs_sin; std::vector row; std::vector col; From bb9624b91abcc77b2ad935ca8602066c21073d66 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Fri, 17 Jul 2026 10:23:56 -0300 Subject: [PATCH 23/24] add coefmatrix.h with CoefMatrix class --- include/trackcpp/coefmatrix.h | 124 ++++++++++++++++++++++++++++++++++ include/trackcpp/elements.h | 13 ++-- include/trackcpp/field3d.h | 41 +++++------ include/trackcpp/field3d.hpp | 60 ++++++++-------- python_package/interface.cpp | 4 +- python_package/interface.h | 5 +- python_package/trackcpp.i | 2 + src/elements.cpp | 9 +-- src/flat_file.cpp | 41 +++++------ 9 files changed, 215 insertions(+), 84 deletions(-) create mode 100644 include/trackcpp/coefmatrix.h diff --git a/include/trackcpp/coefmatrix.h b/include/trackcpp/coefmatrix.h new file mode 100644 index 0000000..6886da2 --- /dev/null +++ b/include/trackcpp/coefmatrix.h @@ -0,0 +1,124 @@ +// TRACKCPP - Particle tracking code +// Copyright (C) 2015 LNLS Accelerator Physics Group +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef COEFMATRIX_H +#define COEFMATRIX_H + +#include +#include + +#ifndef SWIG +template +class MatrixRow { +public: + explicit MatrixRow(T* ptr) + : ptr_(ptr) + {} + + T& operator[](size_t j) + { + return ptr_[j]; + } + +private: + T* ptr_; +}; +#endif + + +class CoefMatrix { +public: + +#ifndef SWIG + using Row = MatrixRow; + using ConstRow = MatrixRow; +#endif + + CoefMatrix() = default; + + CoefMatrix(size_t rows, size_t cols) + : rows_(rows), + cols_(cols), + data_(rows * cols, 0.0) + {} + + void resize(size_t rows, size_t cols) + { + rows_ = rows; + cols_ = cols; + data_.resize(rows * cols); + } + + double& operator()(size_t i, size_t j) + { + return data_[i * cols_ + j]; + } + + const double& operator()(size_t i, size_t j) const + { + return data_[i * cols_ + j]; + } + +#ifndef SWIG + Row operator[](size_t i) + { + return Row(data_.data() + i * cols_); + } + + ConstRow operator[](size_t i) const + { + return ConstRow(data_.data() + i * cols_); + } +#endif + + size_t rows() const + { + return rows_; + } + + size_t cols() const + { + return cols_; + } + + size_t size() const + { + return data_.size(); + } + + bool empty() const + { + return data_.empty(); + } + + double* data() + { + return data_.data(); + } + + const double* data() const + { + return data_.data(); + } + +private: + + size_t rows_ = 0; + size_t cols_ = 0; + std::vector data_; +}; + +#endif \ No newline at end of file diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 62cd1c4..4ab1a05 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -24,6 +24,7 @@ #include #include #include +#include "coefmatrix.h" class Element { @@ -66,8 +67,8 @@ class Element { double field3d_hori_ks = 0; // [1/m] double field3d_hori_kx = 0; // [1/m] double field3d_hori_s0 = 0; // [m] - std::vector> field3d_hori_coefs_cos; - std::vector> field3d_hori_coefs_sin; + CoefMatrix field3d_hori_coefs_cos; + CoefMatrix field3d_hori_coefs_sin; std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; @@ -117,8 +118,8 @@ class Element { static Element rfcavity (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_); static Element kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0); static Element field3d (const std::string& fam_name_, const double& length_, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_, - const std::vector>& field3d_hori_coefs_cos_, - const std::vector>& field3d_hori_coefs_sin_, + const CoefMatrix& field3d_hori_coefs_cos_, + const CoefMatrix& field3d_hori_coefs_sin_, const int nr_steps_ = 40); bool operator==(const Element& o) const; @@ -142,8 +143,8 @@ void initialize_sextupole(Element& element, const double& S, const int& nr_steps void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag); void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks); void initialize_field3d(Element& element, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_, - const std::vector>& field3d_hori_coefs_cos_, - const std::vector>& field3d_hori_coefs_sin_, + const CoefMatrix& field3d_hori_coefs_cos_, + const CoefMatrix& field3d_hori_coefs_sin_, const int nr_steps_); #endif diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h index 90cbc5f..e05cd03 100644 --- a/include/trackcpp/field3d.h +++ b/include/trackcpp/field3d.h @@ -21,36 +21,37 @@ #include #include #include "pos.h" +#include "coefmatrix.h" template T ay(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s); template T ax(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s); template T inty_day_dx(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s); template T intx_dax_dy(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s); @@ -66,16 +67,16 @@ void exp_h1_s(double& s, double step); template void exp_iy_px(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step); template void exp_iy_py(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step); template @@ -89,16 +90,16 @@ void exp_h2_z(const double& beta0, Pos& map, double step); template void exp_ix_px(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step); template void exp_ix_py(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step); @@ -120,15 +121,15 @@ void prop_h3(const double& beta0, Pos& map, double step); template void prop_ix(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step); template void prop_iy(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step); #endif diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 198a82a..7cf5bcf 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -22,13 +22,13 @@ template T ay(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s) { T ay_ = 0.0; - int M = hori_coefs_cos.size(); - int N = hori_coefs_cos[0].size(); + size_t M = hori_coefs_cos.rows(); + size_t N = hori_coefs_cos.cols(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { @@ -45,13 +45,13 @@ T ay(const double& brho, const double& hori_kx, const double& hori_ks, template T ax(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s) { T ax_ = 0.0; - const int M = hori_coefs_cos.size(); - const int N = hori_coefs_cos[0].size(); + size_t M = hori_coefs_cos.rows(); + size_t N = hori_coefs_cos.cols(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { @@ -68,13 +68,13 @@ T ax(const double& brho, const double& hori_kx, const double& hori_ks, template T inty_day_dx(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s) { T day_dx = 0.0; - int M = hori_coefs_cos.size(); - int N = hori_coefs_cos[0].size(); + size_t M = hori_coefs_cos.rows(); + size_t N = hori_coefs_cos.cols(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { @@ -91,13 +91,13 @@ T inty_day_dx(const double& brho, const double& hori_kx, const double& hori_ks, template T intx_dax_dy(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const T& x, const T& y, const double& s) { T dax_dy = 0.0; - int M = hori_coefs_cos.size(); - int N = hori_coefs_cos[0].size(); + size_t M = hori_coefs_cos.rows(); + size_t N = hori_coefs_cos.cols(); for (int m = 1; m <= M; ++m) { for (int n = 1; n <= N; ++n) { @@ -119,8 +119,8 @@ void exp_h1_s(T& s, T step) { template void exp_iy_px(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step) { T factor = inty_day_dx(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.px += factor; @@ -129,8 +129,8 @@ void exp_iy_px(const double& brho, const double& hori_kx, const double& hori_ks, template void exp_iy_py(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step) { T factor = ay(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.py += factor; @@ -150,8 +150,8 @@ void exp_h2_z(Pos& map, const T& pnorm, double step) { template void exp_ix_px(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step) { T factor = ax(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.px += factor; @@ -160,8 +160,8 @@ void exp_ix_px(const double& brho, const double& hori_kx, const double& hori_ks, template void exp_ix_py(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step) { T factor = intx_dax_dy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; map.py += factor; @@ -197,8 +197,8 @@ void prop_h3(Pos& map, const T& pnorm, double step) { template void prop_ix(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step) { exp_ix_px(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); exp_ix_py(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); @@ -207,8 +207,8 @@ void prop_ix(const double& brho, const double& hori_kx, const double& hori_ks, template void prop_iy(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, double s, int sign, double step) { exp_iy_px(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); exp_iy_py(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); @@ -217,8 +217,8 @@ void prop_iy(const double& brho, const double& hori_kx, const double& hori_ks, template void prop_step(const double& brho, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, Pos& map, const T& pnorm, double& s, double step) { prop_h1(map, s, step); prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, +1, step); diff --git a/python_package/interface.cpp b/python_package/interface.cpp index aac5cf4..faf8481 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -243,8 +243,8 @@ Element kickmap_wrapper(const std::string& fam_name_, const std::string& kickta Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& hori_s0_, const double& hori_kx_, const double& hori_ks_, - const std::vector>& hori_coefs_cos_, - const std::vector>& hori_coefs_sin_, + const CoefMatrix& hori_coefs_cos_, + const CoefMatrix& hori_coefs_sin_, const int nr_steps_) { return Element::field3d(fam_name_, length_, hori_s0_, hori_kx_, hori_ks_, hori_coefs_cos_, hori_coefs_sin_, nr_steps_); } diff --git a/python_package/interface.h b/python_package/interface.h index f245f39..40f6595 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -26,6 +26,7 @@ #include #include #include +#include struct LinePassArgs { unsigned int element_offset; @@ -92,8 +93,8 @@ Element sextupole_wrapper(const std::string& fam_name_, const double& length_, c Element rfcavity_wrapper(const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag); Element kickmap_wrapper(const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length = 1.0, const double& rescale_kicks = 1.0); Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& field3d_hori_s0_, const double& field3d_hori_kx_, - const double& field3d_hori_ks_, const std::vector>& field3d_hori_coefs_cos_, - const std::vector>& field3d_hori_coefs_sin_, + const double& field3d_hori_ks_, const CoefMatrix& field3d_hori_coefs_cos_, + const CoefMatrix& field3d_hori_coefs_sin_, const int nr_steps_); Element rbend_wrapper(const std::string& fam_name_, const double& length_, const double& angle_, const double& angle_in_, const double& angle_out_, diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i index 766a189..e3b6a1e 100644 --- a/python_package/trackcpp.i +++ b/python_package/trackcpp.i @@ -29,6 +29,7 @@ #include #include #include +#include #include "interface.h" %} @@ -124,6 +125,7 @@ double get_double_max() { %include "../include/trackcpp/diffusion_matrix.h" %include "../include/trackcpp/naff.h" %include "../include/trackcpp/linalg.h" +%include "../include/trackcpp/coefmatrix.h" %include "interface.h" %template(CppDoublePos) Pos; diff --git a/src/elements.cpp b/src/elements.cpp index 7c853f4..c97b403 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -19,6 +19,7 @@ #include #include #include // necessary for memcmp +#include const std::vector Element::default_polynom = std::vector(3,0); @@ -157,8 +158,8 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt } Element Element::field3d (const std::string& fam_name_, const double& length_, const double& hori_s0_, const double& hori_kx_, const double& hori_ks_, - const std::vector>& hori_coefs_cos_, - const std::vector>& hori_coefs_sin_, + const CoefMatrix& hori_coefs_cos_, + const CoefMatrix& hori_coefs_sin_, const int nr_steps_) { Element e = Element(fam_name_, length_); initialize_field3d(e, hori_s0_, hori_kx_, hori_ks_, hori_coefs_cos_, hori_coefs_sin_, nr_steps_); @@ -319,8 +320,8 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n } void initialize_field3d(Element& element, const double& hori_s0, const double& hori_kx, const double& hori_ks, - const std::vector>& hori_coefs_cos, - const std::vector>& hori_coefs_sin, + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, const int nr_steps) { element.pass_method = PassMethod::pm_field3d_pass; element.nr_steps = nr_steps; diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 3990154..21be2cf 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -30,11 +30,11 @@ static int process_rad_property(std::istringstream& ss); static std::string get_boolean_string(bool value); static bool has_matrix66(const Matrix& r); static bool has_polynom(const std::vector& p); -static bool has_coeffs(const std::vector>& matrix); +static bool has_coeffs(const CoefMatrix& matrix); static void write_6d_vector(std::ostream& fp, const std::string& label, const double* t); static void write_6d_vector(std::ostream& fp, const std::string& label, const std::vector& t); static void write_polynom(std::ostream& fp, const std::string& label, const std::vector& p); -static void write_coeffs(std::ostream& fp, const std::string& label, const std::vector>& matrix); +static void write_coeffs(std::ostream& fp, const std::string& label, const CoefMatrix& matrix); static void synchronize_polynomials(Element& e); static void read_polynomials(std::ifstream& fp, Element& e); static void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator); @@ -318,10 +318,13 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) if (cmd.compare("coefs_cos") == 0 || cmd.compare("coefs_sin") == 0) { - std::vector> Element::* M = nullptr; + CoefMatrix Element::* M = nullptr; - if (cmd == "coefs_cos") M = &Element::field3d_hori_coefs_cos; - if (cmd == "coefs_sin") M = &Element::field3d_hori_coefs_sin; + if (cmd == "coefs_cos") + M = &Element::field3d_hori_coefs_cos; + + if (cmd == "coefs_sin") + M = &Element::field3d_hori_coefs_sin; std::vector row; std::vector col; @@ -346,7 +349,7 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) } if (max_row > 0 && max_col > 0) { - (e.*M).assign(max_row, std::vector(max_col, 0.0)); + (e.*M).resize(max_row, max_col); for (unsigned int k = 0; k < val.size(); ++k) (e.*M)[row[k]][col[k]] = val[k]; @@ -541,11 +544,9 @@ static bool has_matrix66(const Matrix& m) { return false; } -static bool has_coeffs(const std::vector>& matrix) { - if (!matrix.empty()) - return true; - else - return false; +static bool has_coeffs(const CoefMatrix& matrix) +{ + return !matrix.empty(); } static bool has_polynom(const std::vector& p) { @@ -572,20 +573,20 @@ static void write_6d_vector(std::ostream& fp, const std::string& label, const st static void write_coeffs(std::ostream& fp, const std::string& label, - const std::vector>& M) + const CoefMatrix& M) { fp << std::setw(pw) << label; - for (unsigned int i = 0; i < M.size(); ++i) { - for (unsigned int j = 0; j < M[i].size(); ++j) { - - fp.unsetf(std::ios_base::showpos); - fp << i << ' ' << j << ' '; - fp.setf(std::ios_base::showpos); - fp << M[i][j] << ' '; - + for (size_t i = 0; i < M.rows(); ++i) { + for (size_t j = 0; j < M.cols(); ++j) { + + fp.unsetf(std::ios_base::showpos); + fp << i << ' ' << j << ' '; + fp.setf(std::ios_base::showpos); + fp << M[i][j] << ' '; } } + fp << '\n'; } From 81667bd1422d188ce8f113a70b3c4ac34e3adae6 Mon Sep 17 00:00:00 2001 From: Fernando Date: Fri, 17 Jul 2026 15:39:41 -0300 Subject: [PATCH 24/24] ENH: (FIELD3D) improve performance of ID passmethod. --- include/trackcpp/field3d.h | 135 -------------------- include/trackcpp/field3d.hpp | 230 ++++++++++++++++------------------- 2 files changed, 105 insertions(+), 260 deletions(-) delete mode 100644 include/trackcpp/field3d.h diff --git a/include/trackcpp/field3d.h b/include/trackcpp/field3d.h deleted file mode 100644 index e05cd03..0000000 --- a/include/trackcpp/field3d.h +++ /dev/null @@ -1,135 +0,0 @@ -// TRACKCPP - Particle tracking code -// Copyright (C) 2015 LNLS Accelerator Physics Group -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see . - -#ifndef _FIELD3D_H -#define _FIELD3D_H - -#include "auxiliary.h" -#include -#include -#include "pos.h" -#include "coefmatrix.h" - -template -T ay(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s); - - -template -T ax(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s); - - -template -T inty_day_dx(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s); - - -template -T intx_dax_dy(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s); - - -template -T calc_D(const double& beta0, const T& delta); - -template -void exp_h1_z(const double& beta0, Pos& map, double step); - - -void exp_h1_s(double& s, double step); - -template -void exp_iy_px(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step); - - -template -void exp_iy_py(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step); - -template -void exp_h2_y(const double& beta0, Pos& map, double step); - - -template -void exp_h2_z(const double& beta0, Pos& map, double step); - - -template -void exp_ix_px(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step); - - -template -void exp_ix_py(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step); - - -template -void exp_h3_x(const double& beta0, Pos& map, double step); - -template -void exp_h3_z(const double& beta0, Pos& map, double step); - -template -void prop_h1(Pos& map, double& s, double step); - -template -void prop_h2(const double& beta0, Pos& map, double step); - -template -void prop_h3(const double& beta0, Pos& map, double step); - -template -void prop_ix(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step); - -template -void prop_iy(const double& brho, - const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step); - -#endif diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp index 7cf5bcf..b63ce5c 100644 --- a/include/trackcpp/field3d.hpp +++ b/include/trackcpp/field3d.hpp @@ -13,103 +13,108 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see . +#ifndef _FIELD3D_H +#define _FIELD3D_H -#include +#include #include #include #include - - -template -T ay(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s) +#include "auxiliary.h" +#include "pos.h" +#include "coefmatrix.h" + + +inline void calc_sdependence( + const CoefMatrix& hori_coefs_cos, + const CoefMatrix& hori_coefs_sin, + const double& hori_ks, + const double& s, + CoefMatrix& sdependence +) { - T ay_ = 0.0; size_t M = hori_coefs_cos.rows(); size_t N = hori_coefs_cos.cols(); - - for (int m = 1; m <= M; ++m) { - for (int n = 1; n <= N; ++n) { - double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); - double fac1 = hori_coefs_cos[m - 1][n - 1] * (m * hori_kx) / (n * hori_ks * hori_ky); - double fac2 = -hori_coefs_sin[m - 1][n - 1] * (m * hori_kx) / (n * hori_ks * hori_ky); - ay_ += fac1 * sin(m * hori_kx * x) * sinh(hori_ky * y) * std::sin(n * hori_ks * s); - ay_ += fac2 * sin(m * hori_kx * x) * sinh(hori_ky * y) * std::cos(n * hori_ks * s); + for (int n = 1; n <= N; ++n) { + const double nks = n * hori_ks; + const double coss = cos(nks * s); + const double sins = sin(nks * s); + for (int m = 1; m <= M; ++m) { + const double fac1 = hori_coefs_cos[m - 1][n - 1]; + const double fac2 = -hori_coefs_sin[m - 1][n - 1]; + sdependence[m-1][n-1] = fac1 * sins + fac2 * coss; } } - - return 1 * ay_/brho; } template -T ax(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s) +void ay_inty_day_dx( + const double& brho, const double& hori_kx, const double& hori_ks, + const CoefMatrix& sdependence, + const T& x, const T& y, const double& s, + T& ayn, T& idayn +) { - T ax_ = 0.0; - size_t M = hori_coefs_cos.rows(); - size_t N = hori_coefs_cos.cols(); + size_t M = sdependence.rows(); + size_t N = sdependence.cols(); - for (int m = 1; m <= M; ++m) { - for (int n = 1; n <= N; ++n) { - double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); - double fac1 = hori_coefs_cos[m - 1][n - 1] / (n * hori_ks); - double fac2 = -hori_coefs_sin[m - 1][n - 1] / (n * hori_ks); - ax_ += fac1 * cos(m * hori_kx * x) * cosh(hori_ky * y) * std::sin(n * hori_ks * s); - ax_ += fac2 * cos(m * hori_kx * x) * cosh(hori_ky * y) * std::cos(n * hori_ks * s); - } + std::vector nks(N), nks2(N); + for (int n = 1; n <= N; ++n) { + nks[n-1] = n * hori_ks; + nks2[n-1] = nks[n-1] * nks[n-1]; } - return 1 * ax_/brho; -} - -template -T inty_day_dx(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s) -{ - T day_dx = 0.0; - size_t M = hori_coefs_cos.rows(); - size_t N = hori_coefs_cos.cols(); - for (int m = 1; m <= M; ++m) { + const double mkx = m * hori_kx; + const double mkx2 = mkx * mkx; + const T sinx = sin(mkx * x); + const T cosx = cos(mkx * x); for (int n = 1; n <= N; ++n) { - double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); - double fac1 = hori_coefs_cos[m - 1][n - 1] * std::pow(m * hori_kx, 2) / (n * hori_ks * std::pow(hori_ky, 2)); - double fac2 = -hori_coefs_sin[m - 1][n - 1] * std::pow(m * hori_kx, 2) / (n * hori_ks * std::pow(hori_ky, 2)); - day_dx += fac1 * cos(m * hori_kx * x) * std::sin(n * hori_ks * s) * (cosh(hori_ky * y) - 1.0); - day_dx += fac2 * cos(m * hori_kx * x) * std::cos(n * hori_ks * s) * (cosh(hori_ky * y) - 1.0); + const double hori_ky = std::sqrt(mkx2 + nks2[n-1]); + double ratio = mkx / (nks[n-1] * hori_ky); + ayn += sinx * sinh(hori_ky * y) * ratio * sdependence[m-1][n-1]; + + ratio = mkx2 / (nks[n-1] * hori_ky * hori_ky); + idayn += cosx * (cosh(hori_ky * y) - 1.0) * ratio * sdependence[m-1][n-1]; } } - - return 1 * day_dx/brho; + ayn /= brho; + idayn /= brho; } template -T intx_dax_dy(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - const T& x, const T& y, const double& s) +void ax_intx_dax_dy( + const double& brho, const double& hori_kx, const double& hori_ks, + const CoefMatrix& sdependence, + const T& x, const T& y, const double& s, + T& axn, T& idaxn +) { - T dax_dy = 0.0; - size_t M = hori_coefs_cos.rows(); - size_t N = hori_coefs_cos.cols(); + size_t M = sdependence.rows(); + size_t N = sdependence.cols(); + + std::vector nks(N), nks2(N); + for (int n = 1; n <= N; ++n) { + nks[n-1] = n * hori_ks; + nks2[n-1] = nks[n-1] * nks[n-1]; + } for (int m = 1; m <= M; ++m) { + const double mkx = m * hori_kx; + const double mkx2 = mkx * mkx; + const T sinx = sin(mkx * x); + const T cosx = cos(mkx * x); for (int n = 1; n <= N; ++n) { - double hori_ky = std::sqrt(std::pow(m * hori_kx, 2) + std::pow(n * hori_ks, 2)); - double fac1 = hori_coefs_cos[m - 1][n - 1] * hori_ky / (n * hori_ks * m * hori_kx); - double fac2 = -hori_coefs_sin[m - 1][n - 1] * hori_ky / (n * hori_ks * m * hori_kx); - dax_dy += fac1 * sin(m * hori_kx * x) * std::sin(n * hori_ks * s) * sinh(hori_ky * y); - dax_dy += fac2 * sin(m * hori_kx * x) * std::cos(n * hori_ks * s) * sinh(hori_ky * y); + const double hori_ky = std::sqrt(mkx2 + nks2[n-1]); + double ratio = 1 / nks[n-1]; + axn += cosx * cosh(hori_ky * y) * ratio * sdependence[m-1][n-1]; + + ratio = hori_ky / (nks[n-1] * mkx); + idaxn += sinx * sinh(hori_ky * y) * ratio * sdependence[m-1][n-1]; } } - - return 1* dax_dy/brho; + axn /= brho; + idaxn /= brho; } template @@ -117,25 +122,6 @@ void exp_h1_s(T& s, T step) { s += step / 2.0; } -template -void exp_iy_px(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step) { - T factor = inty_day_dx(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; - map.px += factor; -} - - -template -void exp_iy_py(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step) { - T factor = ay(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; - map.py += factor; -} - template void exp_h2_y(Pos& map, const T& pnorm, double step) { map.ry += 0.5*step*pnorm*map.py; @@ -148,26 +134,6 @@ void exp_h2_z(Pos& map, const T& pnorm, double step) { } -template -void exp_ix_px(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step) { - T factor = ax(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; - map.px += factor; -} - - -template -void exp_ix_py(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, - Pos& map, double s, int sign, double step) { - T factor = intx_dax_dy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map.rx, map.ry, s) * sign * -1.0; - map.py += factor; -} - - template void exp_h3_x(Pos& map, const T& pnorm, double step) { map.rx += step*pnorm*map.px; @@ -196,23 +162,31 @@ void prop_h3(Pos& map, const T& pnorm, double step) { } template -void prop_ix(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, +void prop_ix(const double& brho, const double& hori_kx, const double& hori_ks, + const CoefMatrix& sdependence, Pos& map, double s, int sign, double step) { - exp_ix_px(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); - exp_ix_py(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); - + T axn = 0.0; + T idaxn = 0.0; + ax_intx_dax_dy( + brho, hori_kx, hori_ks, sdependence, + map.rx, map.ry, s, axn, idaxn + ); + map.px += axn * sign * -1.0; + map.py += idaxn * sign * -1.0; } template -void prop_iy(const double& brho, const double& hori_kx, const double& hori_ks, - const CoefMatrix& hori_coefs_cos, - const CoefMatrix& hori_coefs_sin, +void prop_iy(const double& brho, const double& hori_kx, const double& hori_ks, + const CoefMatrix& sdependence, Pos& map, double s, int sign, double step) { - exp_iy_px(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); - exp_iy_py(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, sign, step); - + T ayn = 0.0; + T idayn = 0.0; + ay_inty_day_dx( + brho, hori_kx, hori_ks, sdependence, + map.rx, map.ry, s, ayn, idayn + ); + map.px += idayn * sign * -1.0; + map.py += ayn * sign * -1.0; } template @@ -221,14 +195,20 @@ void prop_step(const double& brho, const double& hori_kx, const double& hori_ks, const CoefMatrix& hori_coefs_sin, Pos& map, const T& pnorm, double& s, double step) { prop_h1(map, s, step); - prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, +1, step); + + CoefMatrix sdependence(hori_coefs_cos.rows(), hori_coefs_cos.cols()); + calc_sdependence(hori_coefs_cos, hori_coefs_sin, hori_ks, s, sdependence); + + prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, +1, step); prop_h2(map, pnorm, step); - prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, -1, step); - prop_ix(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, +1, step); + prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, -1, step); + prop_ix(brho, hori_kx, hori_ks, sdependence, map, s, +1, step); prop_h3(map, pnorm, step); - prop_ix(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, -1, step); - prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, +1, step); + prop_ix(brho, hori_kx, hori_ks, sdependence, map, s, -1, step); + prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, +1, step); prop_h2(map, pnorm, step); - prop_iy(brho, hori_kx, hori_ks, hori_coefs_cos, hori_coefs_sin, map, s, -1, step); + prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, -1, step); prop_h1(map, s, step); } + +#endif