Propagator/codegen/sympy_stepper_math.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 // Note: This file is generated by generate_sympy_stepper.py
0010 //       Do not modify it manually.
0011
0012 #pragma once
0013
0014 #include <Acts/Utilities/Result.hpp>
0015
0016 #include <cmath>
0017
0018 template <typename T, typename GetB>
0019 Acts::Result<bool> rk4(const T* p, const T* d, const T t, const T h,
0020                        const T lambda, const T m, const T p_abs, GetB getB,
0021                        T* err, const T errTol, T* new_p, T* new_d, T* new_time,
0022                        T* path_derivatives, T* J) {
0023   const auto B1res = getB(p);
0024   if (!B1res.ok()) {
0025     return Acts::Result<bool>::failure(B1res.error());
0026   }
0027   const auto B1 = *B1res;
0028   const auto x5 = std::pow(h, 2);
0029   const auto x0 = B1[1] * d[2];
0030   const auto x1 = B1[0] * d[2];
0031   const auto x2 = B1[2] * d[0];
0032   const auto x3 = B1[0] * d[1];
0033   const auto x4 = h * d[0];
0034   const auto x7 = h * d[1];
0035   const auto x8 = h * d[2];
0036   const auto x6 = (1.0 / 8.0) * x5;
0037   T k1[3];
0038   k1[0] = lambda * (-x0 + B1[2] * d[1]);
0039   k1[1] = lambda * (x1 - x2);
0040   k1[2] = lambda * (-x3 + B1[1] * d[0]);
0041   T p2[3];
0042   p2[0] = (1.0 / 2.0) * x4 + x6 * k1[0] + p[0];
0043   p2[1] = x6 * k1[1] + (1.0 / 2.0) * x7 + p[1];
0044   p2[2] = x6 * k1[2] + (1.0 / 2.0) * x8 + p[2];
0045   const auto B2res = getB(p2);
0046   if (!B2res.ok()) {
0047     return Acts::Result<bool>::failure(B2res.error());
0048   }
0049   const auto B2 = *B2res;
0050   const auto x9 = (1.0 / 2.0) * h;
0051   const auto x19 = (1.0 / 2.0) * x5;
0052   const auto x11 = lambda * B2[2];
0053   const auto x13 = lambda * B2[1];
0054   const auto x15 = lambda * B2[0];
0055   const auto x20 = x4 + p[0];
0056   const auto x21 = x7 + p[1];
0057   const auto x22 = x8 + p[2];
0058   const auto x10 = x9 * k1[1] + d[1];
0059   const auto x12 = x9 * k1[2] + d[2];
0060   const auto x14 = x9 * k1[0] + d[0];
0061   T k2[3];
0062   k2[0] = x10 * x11 - x12 * x13;
0063   k2[1] = lambda * x12 * B2[0] - x11 * x14;
0064   k2[2] = -x10 * x15 + x13 * x14;
0065   const auto x16 = x9 * k2[1] + d[1];
0066   const auto x17 = x9 * k2[2] + d[2];
0067   const auto x18 = x9 * k2[0] + d[0];
0068   T k3[3];
0069   k3[0] = x11 * x16 - x13 * x17;
0070   k3[1] = lambda * x17 * B2[0] - x11 * x18;
0071   k3[2] = x13 * x18 - x15 * x16;
0072   T p3[3];
0073   p3[0] = x19 * k3[0] + x20;
0074   p3[1] = x19 * k3[1] + x21;
0075   p3[2] = x19 * k3[2] + x22;
0076   const auto B3res = getB(p3);
0077   if (!B3res.ok()) {
0078     return Acts::Result<bool>::failure(B3res.error());
0079   }
0080   const auto B3 = *B3res;
0081   const auto x24 = lambda * B3[2];
0082   const auto x26 = lambda * B3[1];
0083   const auto x28 = lambda * B3[0];
0084   const auto x23 = h * k3[1] + d[1];
0085   const auto x25 = h * k3[2] + d[2];
0086   const auto x27 = h * k3[0] + d[0];
0087   T k4[3];
0088   k4[0] = x23 * x24 - x25 * x26;
0089   k4[1] = lambda * x25 * B3[0] - x24 * x27;
0090   k4[2] = -x23 * x28 + x26 * x27;
0091   const auto x29 = k1[0] + k4[0];
0092   const auto x30 = k1[1] + k4[1];
0093   const auto x31 = k1[2] + k4[2];
0094   *err =
0095       x5 * (std::fabs(-x29 + k2[0] + k3[0]) + std::fabs(-x30 + k2[1] + k3[1]) +
0096             std::fabs(-x31 + k2[2] + k3[2]));
0097   if (*err > errTol) {
0098     return Acts::Result<bool>::success(false);
0099   }
0100   const auto x32 = (1.0 / 6.0) * x5;
0101   new_p[0] = x20 + x32 * (k1[0] + k2[0] + k3[0]);
0102   new_p[1] = x21 + x32 * (k1[1] + k2[1] + k3[1]);
0103   new_p[2] = x22 + x32 * (k1[2] + k2[2] + k3[2]);
0104   const auto x33 = (1.0 / 6.0) * h;
0105   T new_d_tmp[3];
0106   new_d_tmp[0] = x33 * (x29 + 2 * k2[0] + 2 * k3[0]) + d[0];
0107   new_d_tmp[1] = x33 * (x30 + 2 * k2[1] + 2 * k3[1]) + d[1];
0108   new_d_tmp[2] = x33 * (x31 + 2 * k2[2] + 2 * k3[2]) + d[2];
0109   const auto x34 = 1.0 / std::sqrt(std::pow(std::fabs(new_d_tmp[0]), 2) +
0110                                    std::pow(std::fabs(new_d_tmp[1]), 2) +
0111                                    std::pow(std::fabs(new_d_tmp[2]), 2));
0112   new_d[0] = x34 * new_d_tmp[0];
0113   new_d[1] = x34 * new_d_tmp[1];
0114   new_d[2] = x34 * new_d_tmp[2];
0115   const auto x35 = std::pow(m, 2);
0116   const auto dtds = std::sqrt(std::pow(p_abs, 2) + x35) / p_abs;
0117   *new_time = dtds * h + t;
0118   if (J == nullptr) {
0119     return Acts::Result<bool>::success(true);
0120   }
0121   path_derivatives[0] = new_d[0];
0122   path_derivatives[1] = new_d[1];
0123   path_derivatives[2] = new_d[2];
0124   path_derivatives[3] = dtds;
0125   path_derivatives[4] = k4[0];
0126   path_derivatives[5] = k4[1];
0127   path_derivatives[6] = k4[2];
0128   path_derivatives[7] = 0;
0129   const auto x57 = (1.0 / 3.0) * h;
0130   const auto x36 = std::pow(lambda, 2) * x9;
0131   const auto x43 = x1 - x2;
0132   const auto x47 = x11 * x9;
0133   const auto x48 = x15 * x9;
0134   const auto x49 = x13 * x9;
0135   const auto x50 = h * x26;
0136   const auto x51 = h * x24;
0137   const auto x52 = h * x28;
0138   const auto x53 = lambda * x32;
0139   const auto x58 = lambda * x33;
0140   const auto x41 = -x3 + B1[1] * d[0];
0141   const auto x46 = -x0 + B1[2] * d[1];
0142   const auto x37 = x36 * B2[1];
0143   const auto x39 = x36 * B2[2];
0144   const auto x42 = x41 * x9;
0145   const auto x44 = x36 * B2[0];
0146   const auto x54 = x53 * B1[2];
0147   const auto x55 = x53 * B1[1];
0148   const auto x56 = x53 * B1[0];
0149   const auto x59 = x58 * B1[2];
0150   const auto x60 = x58 * B1[1];
0151   const auto x61 = x58 * B1[0];
0152   const auto x38 = x37 * B1[1];
0153   const auto x40 = x39 * B1[2];
0154   const auto x45 = x44 * B1[0];
0155   T dk2dTL[12];
0156   dk2dTL[0] = -x38 - x40;
0157   dk2dTL[1] = -x11 + x44 * B1[1];
0158   dk2dTL[2] = x13 + x44 * B1[2];
0159   dk2dTL[3] = x11 + x37 * B1[0];
0160   dk2dTL[4] = -x40 - x45;
0161   dk2dTL[5] = -x15 + x37 * B1[2];
0162   dk2dTL[6] = -x13 + x39 * B1[0];
0163   dk2dTL[7] = x15 + x39 * B1[1];
0164   dk2dTL[8] = -x38 - x45;
0165   dk2dTL[9] = (1.0 / 2.0) * h * lambda * x43 * B2[2] + x10 * B2[2] -
0166               x12 * B2[1] - x13 * x42;
0167   dk2dTL[10] = x12 * B2[0] - x14 * B2[2] + x15 * x42 - x46 * x47;
0168   dk2dTL[11] = (1.0 / 2.0) * h * lambda * x46 * B2[1] - x10 * B2[0] +
0169                x14 * B2[1] - x43 * x48;
0170   T dk3dTL[12];
0171   dk3dTL[0] = (1.0 / 2.0) * h * lambda * B2[2] * dk2dTL[1] - x49 * dk2dTL[2];
0172   dk3dTL[1] =
0173       (1.0 / 2.0) * h * lambda * B2[0] * dk2dTL[2] - x11 - x47 * dk2dTL[0];
0174   dk3dTL[2] = x13 - x48 * dk2dTL[1] + x49 * dk2dTL[0];
0175   dk3dTL[3] = x11 + x47 * dk2dTL[4] - x49 * dk2dTL[5];
0176   dk3dTL[4] = -x47 * dk2dTL[3] + x48 * dk2dTL[5];
0177   dk3dTL[5] =
0178       (1.0 / 2.0) * h * lambda * B2[1] * dk2dTL[3] - x15 - x48 * dk2dTL[4];
0179   dk3dTL[6] =
0180       (1.0 / 2.0) * h * lambda * B2[2] * dk2dTL[7] - x13 - x49 * dk2dTL[8];
0181   dk3dTL[7] = x15 - x47 * dk2dTL[6] + x48 * dk2dTL[8];
0182   dk3dTL[8] = (1.0 / 2.0) * h * lambda * B2[1] * dk2dTL[6] - x48 * dk2dTL[7];
0183   dk3dTL[9] = (1.0 / 2.0) * h * lambda * B2[2] * dk2dTL[10] + x16 * B2[2] -
0184               x17 * B2[1] - x49 * dk2dTL[11];
0185   dk3dTL[10] = x17 * B2[0] - x18 * B2[2] - x47 * dk2dTL[9] + x48 * dk2dTL[11];
0186   dk3dTL[11] = (1.0 / 2.0) * h * lambda * B2[1] * dk2dTL[9] - x16 * B2[0] +
0187                x18 * B2[1] - x48 * dk2dTL[10];
0188   T dFdTL[12];
0189   dFdTL[0] = h + x32 * dk2dTL[0] + x32 * dk3dTL[0];
0190   dFdTL[1] = x32 * dk2dTL[1] + x32 * dk3dTL[1] - x54;
0191   dFdTL[2] = x32 * dk2dTL[2] + x32 * dk3dTL[2] + x55;
0192   dFdTL[3] = x32 * dk2dTL[3] + x32 * dk3dTL[3] + x54;
0193   dFdTL[4] = h + x32 * dk2dTL[4] + x32 * dk3dTL[4];
0194   dFdTL[5] = x32 * dk2dTL[5] + x32 * dk3dTL[5] - x56;
0195   dFdTL[6] = x32 * dk2dTL[6] + x32 * dk3dTL[6] - x55;
0196   dFdTL[7] = x32 * dk2dTL[7] + x32 * dk3dTL[7] + x56;
0197   dFdTL[8] = h + x32 * dk2dTL[8] + x32 * dk3dTL[8];
0198   dFdTL[9] = x32 * (x46 + dk2dTL[9] + dk3dTL[9]);
0199   dFdTL[10] = x32 * (x43 + dk2dTL[10] + dk3dTL[10]);
0200   dFdTL[11] = x32 * (x41 + dk2dTL[11] + dk3dTL[11]);
0201   T dk4dTL[12];
0202   dk4dTL[0] = h * lambda * B3[2] * dk3dTL[1] - x50 * dk3dTL[2];
0203   dk4dTL[1] = h * lambda * B3[0] * dk3dTL[2] - x24 - x51 * dk3dTL[0];
0204   dk4dTL[2] = x26 + x50 * dk3dTL[0] - x52 * dk3dTL[1];
0205   dk4dTL[3] = x24 - x50 * dk3dTL[5] + x51 * dk3dTL[4];
0206   dk4dTL[4] = -x51 * dk3dTL[3] + x52 * dk3dTL[5];
0207   dk4dTL[5] = h * lambda * B3[1] * dk3dTL[3] - x28 - x52 * dk3dTL[4];
0208   dk4dTL[6] = h * lambda * B3[2] * dk3dTL[7] - x26 - x50 * dk3dTL[8];
0209   dk4dTL[7] = x28 - x51 * dk3dTL[6] + x52 * dk3dTL[8];
0210   dk4dTL[8] = h * lambda * B3[1] * dk3dTL[6] - x52 * dk3dTL[7];
0211   dk4dTL[9] = h * lambda * B3[2] * dk3dTL[10] + x23 * B3[2] - x25 * B3[1] -
0212               x50 * dk3dTL[11];
0213   dk4dTL[10] = x25 * B3[0] - x27 * B3[2] - x51 * dk3dTL[9] + x52 * dk3dTL[11];
0214   dk4dTL[11] = h * lambda * B3[1] * dk3dTL[9] - x23 * B3[0] + x27 * B3[1] -
0215                x52 * dk3dTL[10];
0216   T dGdTL[12];
0217   dGdTL[0] = x33 * dk4dTL[0] + x57 * dk2dTL[0] + x57 * dk3dTL[0] + 1;
0218   dGdTL[1] = x33 * dk4dTL[1] + x57 * dk2dTL[1] + x57 * dk3dTL[1] - x59;
0219   dGdTL[2] = x33 * dk4dTL[2] + x57 * dk2dTL[2] + x57 * dk3dTL[2] + x60;
0220   dGdTL[3] = x33 * dk4dTL[3] + x57 * dk2dTL[3] + x57 * dk3dTL[3] + x59;
0221   dGdTL[4] = x33 * dk4dTL[4] + x57 * dk2dTL[4] + x57 * dk3dTL[4] + 1;
0222   dGdTL[5] = x33 * dk4dTL[5] + x57 * dk2dTL[5] + x57 * dk3dTL[5] - x61;
0223   dGdTL[6] = x33 * dk4dTL[6] + x57 * dk2dTL[6] + x57 * dk3dTL[6] - x60;
0224   dGdTL[7] = x33 * dk4dTL[7] + x57 * dk2dTL[7] + x57 * dk3dTL[7] + x61;
0225   dGdTL[8] = x33 * dk4dTL[8] + x57 * dk2dTL[8] + x57 * dk3dTL[8] + 1;
0226   dGdTL[9] = x33 * x46 + x33 * dk4dTL[9] + x57 * dk2dTL[9] + x57 * dk3dTL[9];
0227   dGdTL[10] =
0228       x33 * x43 + x33 * dk4dTL[10] + x57 * dk2dTL[10] + x57 * dk3dTL[10];
0229   dGdTL[11] =
0230       x33 * x41 + x33 * dk4dTL[11] + x57 * dk2dTL[11] + x57 * dk3dTL[11];
0231   T new_J[64];
0232   new_J[0] = 1;
0233   new_J[1] = 0;
0234   new_J[2] = 0;
0235   new_J[3] = 0;
0236   new_J[4] = 0;
0237   new_J[5] = 0;
0238   new_J[6] = 0;
0239   new_J[7] = 0;
0240   new_J[8] = 0;
0241   new_J[9] = 1;
0242   new_J[10] = 0;
0243   new_J[11] = 0;
0244   new_J[12] = 0;
0245   new_J[13] = 0;
0246   new_J[14] = 0;
0247   new_J[15] = 0;
0248   new_J[16] = 0;
0249   new_J[17] = 0;
0250   new_J[18] = 1;
0251   new_J[19] = 0;
0252   new_J[20] = 0;
0253   new_J[21] = 0;
0254   new_J[22] = 0;
0255   new_J[23] = 0;
0256   new_J[24] = 0;
0257   new_J[25] = 0;
0258   new_J[26] = 0;
0259   new_J[27] = 1;
0260   new_J[28] = 0;
0261   new_J[29] = 0;
0262   new_J[30] = 0;
0263   new_J[31] = 0;
0264   new_J[32] = J[32] * dGdTL[0] + J[40] * dGdTL[1] + J[48] * dGdTL[2] + dFdTL[0];
0265   new_J[33] = J[33] * dGdTL[0] + J[41] * dGdTL[1] + J[49] * dGdTL[2] + dFdTL[1];
0266   new_J[34] = J[34] * dGdTL[0] + J[42] * dGdTL[1] + J[50] * dGdTL[2] + dFdTL[2];
0267   new_J[35] = 0;
0268   new_J[36] = J[36] * dGdTL[0] + J[44] * dGdTL[1] + J[52] * dGdTL[2];
0269   new_J[37] = J[37] * dGdTL[0] + J[45] * dGdTL[1] + J[53] * dGdTL[2];
0270   new_J[38] = J[38] * dGdTL[0] + J[46] * dGdTL[1] + J[54] * dGdTL[2];
0271   new_J[39] = 0;
0272   new_J[40] = J[32] * dGdTL[3] + J[40] * dGdTL[4] + J[48] * dGdTL[5] + dFdTL[3];
0273   new_J[41] = J[33] * dGdTL[3] + J[41] * dGdTL[4] + J[49] * dGdTL[5] + dFdTL[4];
0274   new_J[42] = J[34] * dGdTL[3] + J[42] * dGdTL[4] + J[50] * dGdTL[5] + dFdTL[5];
0275   new_J[43] = 0;
0276   new_J[44] = J[36] * dGdTL[3] + J[44] * dGdTL[4] + J[52] * dGdTL[5];
0277   new_J[45] = J[37] * dGdTL[3] + J[45] * dGdTL[4] + J[53] * dGdTL[5];
0278   new_J[46] = J[38] * dGdTL[3] + J[46] * dGdTL[4] + J[54] * dGdTL[5];
0279   new_J[47] = 0;
0280   new_J[48] = J[32] * dGdTL[6] + J[40] * dGdTL[7] + J[48] * dGdTL[8] + dFdTL[6];
0281   new_J[49] = J[33] * dGdTL[6] + J[41] * dGdTL[7] + J[49] * dGdTL[8] + dFdTL[7];
0282   new_J[50] = J[34] * dGdTL[6] + J[42] * dGdTL[7] + J[50] * dGdTL[8] + dFdTL[8];
0283   new_J[51] = 0;
0284   new_J[52] = J[36] * dGdTL[6] + J[44] * dGdTL[7] + J[52] * dGdTL[8];
0285   new_J[53] = J[37] * dGdTL[6] + J[45] * dGdTL[7] + J[53] * dGdTL[8];
0286   new_J[54] = J[38] * dGdTL[6] + J[46] * dGdTL[7] + J[54] * dGdTL[8];
0287   new_J[55] = 0;
0288   new_J[56] = J[32] * dGdTL[9] + J[40] * dGdTL[10] + J[48] * dGdTL[11] + J[56] +
0289               dFdTL[9];
0290   new_J[57] = J[33] * dGdTL[9] + J[41] * dGdTL[10] + J[49] * dGdTL[11] + J[57] +
0291               dFdTL[10];
0292   new_J[58] = J[34] * dGdTL[9] + J[42] * dGdTL[10] + J[50] * dGdTL[11] + J[58] +
0293               dFdTL[11];
0294   new_J[59] = J[59] + h * lambda * x35 / dtds;
0295   new_J[60] = J[36] * dGdTL[9] + J[44] * dGdTL[10] + J[52] * dGdTL[11] + J[60];
0296   new_J[61] = J[37] * dGdTL[9] + J[45] * dGdTL[10] + J[53] * dGdTL[11] + J[61];
0297   new_J[62] = J[38] * dGdTL[9] + J[46] * dGdTL[10] + J[54] * dGdTL[11] + J[62];
0298   new_J[63] = 1;
0299   J[0] = new_J[0];
0300   J[1] = new_J[1];
0301   J[2] = new_J[2];
0302   J[3] = new_J[3];
0303   J[4] = new_J[4];
0304   J[5] = new_J[5];
0305   J[6] = new_J[6];
0306   J[7] = new_J[7];
0307   J[8] = new_J[8];
0308   J[9] = new_J[9];
0309   J[10] = new_J[10];
0310   J[11] = new_J[11];
0311   J[12] = new_J[12];
0312   J[13] = new_J[13];
0313   J[14] = new_J[14];
0314   J[15] = new_J[15];
0315   J[16] = new_J[16];
0316   J[17] = new_J[17];
0317   J[18] = new_J[18];
0318   J[19] = new_J[19];
0319   J[20] = new_J[20];
0320   J[21] = new_J[21];
0321   J[22] = new_J[22];
0322   J[23] = new_J[23];
0323   J[24] = new_J[24];
0324   J[25] = new_J[25];
0325   J[26] = new_J[26];
0326   J[27] = new_J[27];
0327   J[28] = new_J[28];
0328   J[29] = new_J[29];
0329   J[30] = new_J[30];
0330   J[31] = new_J[31];
0331   J[32] = new_J[32];
0332   J[33] = new_J[33];
0333   J[34] = new_J[34];
0334   J[35] = new_J[35];
0335   J[36] = new_J[36];
0336   J[37] = new_J[37];
0337   J[38] = new_J[38];
0338   J[39] = new_J[39];
0339   J[40] = new_J[40];
0340   J[41] = new_J[41];
0341   J[42] = new_J[42];
0342   J[43] = new_J[43];
0343   J[44] = new_J[44];
0344   J[45] = new_J[45];
0345   J[46] = new_J[46];
0346   J[47] = new_J[47];
0347   J[48] = new_J[48];
0348   J[49] = new_J[49];
0349   J[50] = new_J[50];
0350   J[51] = new_J[51];
0351   J[52] = new_J[52];
0352   J[53] = new_J[53];
0353   J[54] = new_J[54];
0354   J[55] = new_J[55];
0355   J[56] = new_J[56];
0356   J[57] = new_J[57];
0357   J[58] = new_J[58];
0358   J[59] = new_J[59];
0359   J[60] = new_J[60];
0360   J[61] = new_J[61];
0361   J[62] = new_J[62];
0362   J[63] = new_J[63];
0363   return Acts::Result<bool>::success(true);
0364 }