15 EngeCoefs_entry_m(NULL),
16 EngeCoefs_exit_m(NULL),
20 cosExitRotation_m(1.0),
21 sinExitRotation_m(0.0) {
24 std::string tmpString;
32 bool parsing_passed = \
33 interpreteLine<std::string, int, int, double>(file,
38 parsing_passed = parsing_passed &&
39 interpreteLine<double, double, double, int>(file,
44 parsing_passed = parsing_passed &&
45 interpreteLine<double, double, double, int>(file,
50 for(
int i = 0; (i < num_gridpz + 1) && parsing_passed; ++ i) {
51 parsing_passed = parsing_passed &&
52 interpreteLine<double>(file, tmpDouble);
55 parsing_passed = parsing_passed &&
63 zend_exit_m = zbegin_entry_m - 1
e-3;
64 zend_entry_m = zbegin_entry_m - 1e-3;
65 zbegin_exit_m = zbegin_entry_m - 1e-3;
68 zbegin_entry_m /= 100.;
70 polynomialOrigin_entry_m /= 100.;
71 zbegin_exit_m /= 100.;
73 polynomialOrigin_exit_m /= 100.;
95 double tolerance = 1
e-8;
100 std::string tmpString;
103 double minValue = 99999999.99, maxValue = -99999999.99;
105 int num_gridp_fringe_entry, num_gridp_fringe_exit;
106 int num_gridp_before_fringe_entry, num_gridp_before_fringe_exit;
108 double *rightHandSide;
109 double *leastSquareMatrix;
112 interpreteLine<std::string, int, int, double>(in, tmpString, tmpInt, tmpInt, tmpDouble);
113 interpreteLine<double, double, double, int>(in, tmpDouble, tmpDouble, tmpDouble, num_gridpz);
114 interpreteLine<double, double, double, int>(in, tmpDouble, tmpDouble, tmpDouble, tmpInt);
121 RealValues =
new double[num_gridpz + 1];
123 for(
int i = 0; i < num_gridpz + 1; ++i) {
124 interpreteLine<double>(in, RealValues[i]);
125 if(RealValues[i] > maxValue) maxValue = RealValues[i];
126 else if(RealValues[i] < minValue) minValue = RealValues[i];
131 for(
int i = 0; i < num_gridpz + 1; ++i)
132 RealValues[i] = (RealValues[i] - minValue) / (maxValue - minValue);
136 while(i < num_gridpz + 1 && RealValues[i] < tolerance) ++i;
137 num_gridp_before_fringe_entry = i - 1;
140 while(i < num_gridpz + 1 && RealValues[i] < 1. - tolerance) ++i;
141 num_gridp_fringe_entry = i - 1 - num_gridp_before_fringe_entry;
144 while(i < num_gridpz + 1 && RealValues[i] >= 1. - tolerance) ++i;
145 num_gridp_before_fringe_exit = i - 1;
147 while(i < num_gridpz + 1 && RealValues[i] > tolerance) ++i;
148 num_gridp_fringe_exit = i - 1 - num_gridp_before_fringe_exit;
152 int num_gridp_fringe =
max(num_gridp_fringe_entry, num_gridp_fringe_exit);
154 leastSquareMatrix =
new double[(polynomialOrder + 1) * num_gridp_fringe];
155 rightHandSide =
new double[num_gridp_fringe];
159 for(
int i = 0; i < num_gridp_fringe_entry; ++i) {
160 double powerOfZ = 1.;
162 rightHandSide[i] =
log(1. / RealValues[num_gridp_before_fringe_entry + i + 1] - 1.);
164 leastSquareMatrix[i * (polynomialOrder_entry_m + 1) + j] = powerOfZ;
172 for(
int i = 0; i < num_gridp_fringe_exit; ++i) {
173 double powerOfZ = 1.;
175 rightHandSide[i] =
log(1. / RealValues[num_gridp_before_fringe_exit + i + 1] - 1.);
177 leastSquareMatrix[i * (polynomialOrder_exit_m + 1) + j] = powerOfZ;
189 delete[] leastSquareMatrix;
190 delete[] rightHandSide;
227 double S, dSdz, d2Sdz2 = 0.0;
228 double expS, f, dfdz, d2fdz2;
246 S = EngeCoefs[polynomialOrder] * z;
247 S += EngeCoefs[polynomialOrder - 1];
248 dSdz = polynomialOrder * EngeCoefs[polynomialOrder];
250 for(
int i = polynomialOrder - 2; i >= 0; i--) {
251 S = S * z + EngeCoefs[i];
252 dSdz = dSdz * z + (i + 1) * EngeCoefs[i + 1];
253 d2Sdz2 = d2Sdz2 * z + (i + 2) * (i + 1) * EngeCoefs[i + 2];
256 f = 1.0 / (1.0 + expS);
259 dfdz = - f * ((f * expS) * dSdz);
262 d2fdz2 = ((-d2Sdz2 - dSdz * dSdz * (1. - 2. * (expS * f))) * (f * expS) * f) / (
gapHeight_m *
gapHeight_m);
266 strength(2) = d2fdz2;
362 zExit_m = polynomialOrigin_entry_m + deltaZ *
cos(bendAngle / 2.0);
369 namespace QRDecomposition {
370 void solve(
double *
Matrix,
double *Solution,
double *rightHandSide,
const int &M,
const int &N) {
375 double *
R =
new double[M * N];
376 double *tempVector =
new double[M];
377 double *residuum =
new double[M];
379 for(
int i = 0; i < M; ++i) {
380 for(
int j = 0; j < N; ++j)
381 R[i * N + j] = Matrix[i * N + j];
382 tempVector[i] = rightHandSide[i];
386 for(
int i = 0; i < N; ++i) {
387 for(
int j = i + 1; j < M; ++j) {
388 len =
sqrt(R[j * N + i] * R[j * N + i] + R[i * (N + 1)] * R[i * (N + 1)]);
389 sinphi = -R[j * N + i] / len;
390 cosphi = R[i * (N + 1)] / len;
392 for(
int k = 0; k < N; ++k) {
393 tempValue = cosphi * R[ i * N + k] - sinphi * R[ j * N + k];
394 R[j * N + k] = sinphi * R[ i * N + k] + cosphi * R[ j * N + k];
395 R[i * N + k] = tempValue;
403 for(
int i = 0; i < N; ++i) {
405 for(
int j = 0; j < M; ++j) {
406 tempValue += Matrix[j * N + i] * rightHandSide[j];
408 Solution[i] = tempValue;
411 for(
int i = 0; i < N; ++i) {
413 for(
int j = 0; j < i; ++j)
414 tempValue += R[j * N + i] * residuum[j];
415 residuum[i] = (Solution[i] - tempValue) / R[i * (N + 1)];
418 for(
int i = N - 1; i >= 0; --i) {
420 for(
int j = N - 1; j > i; --j)
421 tempValue += R[i * N + j] * Solution[j];
422 Solution[i] = (residuum[i] - tempValue) / R[i * (N + 1)];
425 for(
int i = 0; i < M; ++i) {
427 for(
int j = 0; j < N; ++j)
428 tempValue += Matrix[i * N + j] * Solution[j];
429 residuum[i] = rightHandSide[i] - tempValue;
432 for(
int i = 0; i < N; ++i) {
434 for(
int j = 0; j < M; ++j)
435 tempValue += Matrix[j * N + i] * residuum[j];
436 tempVector[i] = tempValue;
439 for(
int i = 0; i < N; ++i) {
441 for(
int j = 0; j < i; ++j)
442 tempValue += R[j * N + i] * residuum[j];
443 residuum[i] = (tempVector[i] - tempValue) / R[i * (N + 1)];
446 for(
int i = N - 1; i >= 0; --i) {
448 for(
int j = N - 1; j > i; --j)
449 tempValue += R[i * N + j] * tempVector[j];
450 tempVector[i] = (residuum[i] - tempValue) / R[i * (N + 1)];
451 Solution[i] += tempVector[i];
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const
double cosExitRotation_m
Cos and sin of the exit edge rotation with respect to the local coordinates.
double polynomialOrigin_exit_m
constexpr double e
The value of .
virtual bool getFieldDerivative(const Vector_t &X, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
FM1DProfile2(std::string aFilename)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Tps< T > sin(const Tps< T > &x)
Sine.
virtual void getInfo(Inform *)
constexpr double two_pi
The value of .
Tps< T > exp(const Tps< T > &x)
Exponential.
double * EngeCoefs_exit_m
void disableFieldmapWarning()
static std::string typeset_msg(const std::string &msg, const std::string &title)
int polynomialOrder_entry_m
Tps< T > log(const Tps< T > &x)
Natural logarithm.
constexpr double mu_0
The permeability of vacuum in Vs/Am.
constexpr double c
The velocity of light in m/s.
bool interpreteEOF(std::ifstream &in)
Vektor< double, 3 > Vector_t
virtual double getFrequency() const
Tps< T > sqrt(const Tps< T > &x)
Square root.
double polynomialOrigin_entry_m
virtual void setExitFaceSlope(const double &)
Tps< T > cos(const Tps< T > &x)
Cosine.
double * EngeCoefs_entry_m
virtual void setFrequency(double freq)
virtual void setEdgeConstants(const double &bendAngle, const double &entranceAngle, const double &exitAngle)
Inform & level3(Inform &inf)
void solve(double *Matrix, double *Solution, double *rightHandSide, const int &M, const int &N)
int polynomialOrder_exit_m
virtual bool getFieldstrength(const Vector_t &X, Vector_t &strength, Vector_t &info) const
Inform & endl(Inform &inf)