FM1DProfile2.cpp

Go to the documentation of this file.

#include "Fields/FM1DProfile2.h"
#include "Fields/Fieldmap.hpp"
#include "Physics/Physics.h"
#include "Physics/Units.h"

#include <fstream>
#include <ios>
#include <cmath>

_FM1DProfile2::_FM1DProfile2(const std::string& filename)
    : _Fieldmap(filename),
      EngeCoefs_entry_m(nullptr),
      EngeCoefs_exit_m(nullptr),
      exit_slope_m(0.0),
      xExit_m(0.0),
      zExit_m(0.0),
      cosExitRotation_m(1.0),
      sinExitRotation_m(0.0) {

    int tmpInt;
    std::string tmpString;
    double tmpDouble;
    int num_gridpz = -1;

    Type = T1DProfile2;
    std::ifstream file(Filename_m.c_str());

    if (file.good()) {
        bool parsing_passed =                               \
                interpretLine<std::string, int, int, double>(file,
                        tmpString,
                        polynomialOrder_entry_m,
                        polynomialOrder_exit_m,
                        gapHeight_m);
        parsing_passed = parsing_passed &&
                         interpretLine<double, double, double, int>(file,
                                 zbegin_entry_m,
                                 polynomialOrigin_entry_m,
                                 zend_entry_m,
                                 num_gridpz);
        parsing_passed = parsing_passed &&
                         interpretLine<double, double, double, int>(file,
                                 zbegin_exit_m,
                                 polynomialOrigin_exit_m,
                                 zend_exit_m,
                                 tmpInt);
        for (int i = 0; (i < num_gridpz + 1) && parsing_passed; ++ i) {
            parsing_passed = parsing_passed &&
                             interpretLine<double>(file, tmpDouble);
        }

        parsing_passed = parsing_passed &&
                         interpreteEOF(file);

        file.close();
        lines_read_m = 0;

        if (!parsing_passed) {
            disableFieldmapWarning();
            zend_exit_m = zbegin_entry_m - 1e-3;
            zend_entry_m = zbegin_entry_m - 1e-3;
            zbegin_exit_m = zbegin_entry_m - 1e-3;
        } else {
            // conversion cm to m
            zbegin_entry_m *= Units::cm2m;
            zend_entry_m *= Units::cm2m;
            polynomialOrigin_entry_m *= Units::cm2m;
            zbegin_exit_m *= Units::cm2m;
            zend_exit_m *= Units::cm2m;
            polynomialOrigin_exit_m *= Units::cm2m;
            gapHeight_m *= Units::cm2m;
        }
        length_m = zend_exit_m - zbegin_entry_m;
    } else {
        noFieldmapWarning();
        zbegin_entry_m = 0.0;
        zend_exit_m = zbegin_entry_m - 1e-3;
        zend_entry_m = zbegin_entry_m - 1e-3;
        zbegin_exit_m = zbegin_entry_m - 1e-3;
    }
}

_FM1DProfile2::~_FM1DProfile2() {
    if (EngeCoefs_entry_m != nullptr) {
        delete[] EngeCoefs_entry_m;
        delete[] EngeCoefs_exit_m;
    }
}

FM1DProfile2 _FM1DProfile2::create(const std::string& filename)
{
    return FM1DProfile2(new _FM1DProfile2(filename));
}

void _FM1DProfile2::readMap() {
    if (EngeCoefs_entry_m == nullptr) {
        double tolerance = 1e-8;

        std::ifstream in(Filename_m.c_str());

        int tmpInt;
        std::string tmpString;
        double tmpDouble;

        double minValue = 99999999.99, maxValue = -99999999.99;
        int num_gridpz;
        int num_gridp_fringe_entry, num_gridp_fringe_exit;
        int num_gridp_before_fringe_entry, num_gridp_before_fringe_exit;
        double *RealValues;
        double *rightHandSide;
        double *leastSquareMatrix;
        double dZ;

        interpretLine<std::string, int, int, double>(in, tmpString, tmpInt, tmpInt, tmpDouble);
        interpretLine<double, double, double, int>(in, tmpDouble, tmpDouble, tmpDouble, num_gridpz);
        interpretLine<double, double, double, int>(in, tmpDouble, tmpDouble, tmpDouble, tmpInt);

        dZ = length_m / num_gridpz;

        EngeCoefs_entry_m = new double[polynomialOrder_entry_m + 1];
        EngeCoefs_exit_m = new double[polynomialOrder_exit_m + 1];

        RealValues = new double[num_gridpz + 1];

        for (int i = 0; i < num_gridpz + 1; ++i) {
            interpretLine<double>(in, RealValues[i]);
            if (RealValues[i] > maxValue) maxValue = RealValues[i];
            else if (RealValues[i] < minValue) minValue = RealValues[i];
        }
        in.close();

        // normalise the values //
        for (int i = 0; i < num_gridpz + 1; ++i)
            RealValues[i] = (RealValues[i] - minValue) / (maxValue - minValue);

        // find begin of entry fringe field //
        int i = 0;
        while(i < num_gridpz + 1 && RealValues[i] < tolerance) ++i;
        num_gridp_before_fringe_entry = i - 1;

        // find end of entry fringe field //
        while(i < num_gridpz + 1 && RealValues[i] < 1. - tolerance) ++i;
        num_gridp_fringe_entry = i - 1 - num_gridp_before_fringe_entry;

        // find begin of exit fringe field //
        while(i < num_gridpz + 1 && RealValues[i] >= 1. - tolerance) ++i;
        num_gridp_before_fringe_exit = i - 1;

        while(i < num_gridpz + 1 && RealValues[i] > tolerance) ++i;
        num_gridp_fringe_exit = i - 1 - num_gridp_before_fringe_exit;

        // set the origin of the polynomials

        int num_gridp_fringe = std::max(num_gridp_fringe_entry, num_gridp_fringe_exit);
        int polynomialOrder = std::max(polynomialOrder_entry_m, polynomialOrder_exit_m);
        leastSquareMatrix = new double[(polynomialOrder + 1) * num_gridp_fringe];
        rightHandSide = new double[num_gridp_fringe];


        double first = polynomialOrigin_entry_m - zbegin_entry_m - dZ * (num_gridp_before_fringe_entry + 1);
        for (int i = 0; i < num_gridp_fringe_entry; ++i) {
            double powerOfZ = 1.;
            double Z = (first - dZ * i) / gapHeight_m;
            rightHandSide[i] = log(1. / RealValues[num_gridp_before_fringe_entry + i + 1] - 1.);
            for (int j = 0; j < polynomialOrder_entry_m + 1; ++j) {
                leastSquareMatrix[i * (polynomialOrder_entry_m + 1) + j] = powerOfZ;
                powerOfZ *= Z;
            }
        }

        QRDecomposition::solve(leastSquareMatrix, EngeCoefs_entry_m, rightHandSide, num_gridp_fringe_entry, polynomialOrder_entry_m + 1);

        first = polynomialOrigin_exit_m - zbegin_entry_m - dZ * (num_gridp_before_fringe_exit + 1);
        for (int i = 0; i < num_gridp_fringe_exit; ++i) {
            double powerOfZ = 1.;
            double Z = (dZ * i - first) / gapHeight_m;
            rightHandSide[i] = log(1. / RealValues[num_gridp_before_fringe_exit + i + 1] - 1.);
            for (int j = 0; j < polynomialOrder_exit_m + 1; ++j) {
                leastSquareMatrix[i * (polynomialOrder_exit_m + 1) + j] = powerOfZ;
                powerOfZ *= Z;
            }
        }
        QRDecomposition::solve(leastSquareMatrix, EngeCoefs_exit_m, rightHandSide, num_gridp_fringe_exit, polynomialOrder_exit_m + 1);

        zbegin_exit_m = zbegin_entry_m + num_gridp_before_fringe_exit * dZ;
        zend_exit_m = zbegin_exit_m + num_gridp_fringe_exit * dZ;
        zbegin_entry_m += num_gridp_before_fringe_entry * dZ;
        zend_entry_m = zbegin_entry_m + num_gridp_fringe_entry * dZ;

        delete[] RealValues;
        delete[] leastSquareMatrix;
        delete[] rightHandSide;

        INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m  + "'", "info") << "\n"
                << endl);

    }
}

void _FM1DProfile2::freeMap() {
    if (EngeCoefs_entry_m != nullptr) {

        delete[] EngeCoefs_entry_m;
        delete[] EngeCoefs_exit_m;
    }
}

bool _FM1DProfile2::getFieldstrength(const Vector_t &R, Vector_t &strength, Vector_t &info) const {

    info = Vector_t(0.0);

    // Find coordinates in the entrance frame.
    Vector_t REntrance(R(0), 0.0, R(2) + zbegin_entry_m);

    // Find coordinates in the exit frame.
    Vector_t RExit(0.0, R(1), 0.0);

    RExit(0) = (R(0) - xExit_m) * cosExitRotation_m - (R(2) + zbegin_entry_m - zExit_m) * sinExitRotation_m;
    RExit(2) = (R(0) - xExit_m) * sinExitRotation_m + (R(2) + zbegin_entry_m - zExit_m) * cosExitRotation_m + polynomialOrigin_exit_m;


    if (REntrance(2) >= zend_entry_m && RExit(2) <= zbegin_exit_m) {
        strength = Vector_t(1.0, 0.0, 0.0);
    } else {
        double d2Sdz2 = 0.0;
        double z;
        double *EngeCoefs;
        int polynomialOrder;
        info(0) = 1.0;
        if (REntrance(2) >= zbegin_entry_m && REntrance(2) < zend_entry_m) {
            z = -(REntrance(2) - polynomialOrigin_entry_m) / gapHeight_m;
            EngeCoefs = EngeCoefs_entry_m;
            polynomialOrder = polynomialOrder_entry_m;
        } else if (RExit(2) > zbegin_exit_m && RExit(2) <= zend_exit_m) {
            z = (RExit(2) - polynomialOrigin_exit_m) / gapHeight_m;
            EngeCoefs = EngeCoefs_exit_m;
            polynomialOrder = polynomialOrder_exit_m;
            info(1) = 1.0;
        } else {
            return true;
        }

        double S = EngeCoefs[polynomialOrder] * z;
        S += EngeCoefs[polynomialOrder - 1];
        double dSdz = polynomialOrder * EngeCoefs[polynomialOrder];

        for (int i = polynomialOrder - 2; i >= 0; i--) {
            S = S * z + EngeCoefs[i];
            dSdz = dSdz * z + (i + 1) * EngeCoefs[i + 1];
            d2Sdz2 = d2Sdz2 * z + (i + 2) * (i + 1) * EngeCoefs[i + 2];
        }
        double expS = std::exp(S);
        double f = 1.0 / (1.0 + expS);
        if (f > 1.e-30) {
            // First derivative of Enge function, f.
            double dfdz = - f * ((f * expS) * dSdz);

            // Second derivative of Enge functioin, f.
            double d2fdz2 = ((-d2Sdz2 - dSdz * dSdz * (1. - 2. * (expS * f))) * (f * expS) * f) / (gapHeight_m * gapHeight_m);

            strength(0) = f;
            strength(1) = dfdz / gapHeight_m;
            strength(2) = d2fdz2;
        } else {
            strength = Vector_t(0.0);
        }

    }
    info(2) = exit_slope_m;

    return true;

    //    info = Vector_t(0.0);
    //    const Vector_t tmpR(R(0), R(1), R(2) + zbegin_entry_m);
    //
    //    if (tmpR(2) >= zend_entry_m && tmpR(2) <= exit_slope_m * tmpR(0) + zbegin_exit_m) {
    //        strength = Vector_t(1.0, 0.0, 0.0);
    //        info(0) = 3.0;
    //    } else {
    //        double S, dSdz, d2Sdz2 = 0.0;
    //        double expS, f, dfdz, d2fdz2;
    //        double z;
    //        double *EngeCoefs;
    //        int polynomialOrder;
    //        if (tmpR(2) >= zbegin_entry_m && tmpR(2) < zend_entry_m) {
    //            z = -(tmpR(2) - polynomialOrigin_entry_m) / gapHeight_m;
    //            EngeCoefs = EngeCoefs_entry_m;
    //            polynomialOrder = polynomialOrder_entry_m;
    //            info(0) = 1.0;
    //        } else if (tmpR(2) > exit_slope_m * tmpR(0) + zbegin_exit_m && tmpR(2) <= exit_slope_m * tmpR(0) + zend_exit_m) {
    //            z = (tmpR(2) - exit_slope_m * tmpR(0) - polynomialOrigin_exit_m) / sqrt(exit_slope_m * exit_slope_m + 1) / gapHeight_m;
    //            EngeCoefs = EngeCoefs_exit_m;
    //            polynomialOrder = polynomialOrder_exit_m;
    //            info(0) = 2.0;
    //        } else {
    //            return true;
    //        }
    //
    //        S = EngeCoefs[polynomialOrder] * z;
    //        S += EngeCoefs[polynomialOrder - 1];
    //        dSdz = polynomialOrder * EngeCoefs[polynomialOrder];
    //
    //        for (int i = polynomialOrder - 2; i >= 0; i--) {
    //            S = S * z + EngeCoefs[i];
    //            dSdz = dSdz * z + (i + 1) * EngeCoefs[i+1];
    //            d2Sdz2 = d2Sdz2 * z + (i + 2) * (i + 1) * EngeCoefs[i+2];
    //        }
    //        expS = exp(S);
    //        f = 1.0 / (1.0 + expS);
    //        if (f > 1.e-30) {
    //            dfdz = - f * ((f * expS) * dSdz); // first derivative of f
    //            d2fdz2 = ((-d2Sdz2 - dSdz * dSdz * (1. - 2. * (expS * f))) * (f * expS) * f) / (gapHeight_m * gapHeight_m);  // second derivative of f
    //
    //            strength(0) = f;
    //            strength(1) = dfdz / gapHeight_m;
    //            strength(2) = d2fdz2;
    //        } else {
    //            strength = Vector_t(0.0);
    //        }
    //
    //    }
    //    info(1) = exit_slope_m;
    //    return true;

}

bool _FM1DProfile2::getFieldDerivative(const Vector_t &/*R*/, Vector_t &/*E*/, Vector_t &/*B*/, const DiffDirection &/*dir*/) const {
    return false;
}

void _FM1DProfile2::getFieldDimensions(double &zBegin, double &zEnd) const {
    zBegin = zbegin_entry_m;
    zEnd = zend_exit_m;
}
void _FM1DProfile2::getFieldDimensions(double &/*xIni*/, double &/*xFinal*/, double &/*yIni*/, double &/*yFinal*/, double &/*zIni*/, double &/*zFinal*/) const {}

void _FM1DProfile2::swap()
{}

void _FM1DProfile2::getInfo(Inform *msg) {
    (*msg) << Filename_m << " (1D Profile type 2); zini= " << zbegin_entry_m << " m; zfinal= " << zend_exit_m << " m;" << endl;
}

double _FM1DProfile2::getFrequency() const {
    return 0.0;
}

void _FM1DProfile2::setFrequency(double /*freq*/)
{}

void _FM1DProfile2::setExitFaceSlope(const double &m) {
    exit_slope_m = m;
}

void _FM1DProfile2::setEdgeConstants(const double &bendAngle, const double &entranceAngle, const double &exitAngle) {

    double deltaZ = polynomialOrigin_exit_m - polynomialOrigin_entry_m;

    zExit_m = polynomialOrigin_entry_m + deltaZ * std::cos(bendAngle / 2.0);
    xExit_m = -deltaZ * std::sin(bendAngle / 2.0);

    cosExitRotation_m = std::cos(-bendAngle + entranceAngle + exitAngle);
    sinExitRotation_m = std::sin(-bendAngle + entranceAngle + exitAngle);
}

namespace QRDecomposition {
    void solve(double *Matrix, double *Solution, double *rightHandSide, const int &M, const int &N) {
        double sinphi;
        double cosphi;
        double tempValue;
        double len;
        double *R = new double[M * N];
        double *tempVector = new double[M];
        double *residuum = new double[M];

        for (int i = 0; i < M; ++i) {
            for (int j = 0; j < N; ++j)
                R[i * N + j] = Matrix[i * N + j];
            tempVector[i] = rightHandSide[i];
        }

        /* using Givens rotations */
        for (int i = 0; i < N; ++i) {
            for (int j = i + 1; j < M; ++j) {
                len = std::sqrt(R[j * N + i] * R[j * N + i] + R[i * (N + 1)] * R[i * (N + 1)]);
                sinphi = -R[j * N + i] / len;
                cosphi = R[i * (N + 1)] / len;

                for (int k = 0; k < N; ++k) {
                    tempValue = cosphi * R[ i * N + k] - sinphi * R[ j * N + k];
                    R[j * N + k] = sinphi * R[ i * N + k] + cosphi * R[ j * N + k];
                    R[i * N + k] = tempValue;
                }
            }
        }

        /* one step of iterative refinement */

        //     cout << "A^T*b" << endl;
        for (int i = 0; i < N; ++i) {     /* A^T*b */
            tempValue = 0.0;
            for (int j = 0; j < M; ++j) {
                tempValue += Matrix[j * N + i] * rightHandSide[j];
            }
            Solution[i] = tempValue;
        }
        //     cout << "R^-TA^T*b" << endl;
        for (int i = 0; i < N; ++i) {    /* R^-T*A^T*b */
            tempValue = 0.0;
            for (int j = 0; j < i; ++j)
                tempValue += R[j * N + i] * residuum[j];
            residuum[i] = (Solution[i] - tempValue) / R[i * (N + 1)];
        }
        //     cout << "R^-1R^-TA^T*b" << endl;
        for (int i = N - 1; i >= 0; --i) { /* R^-1*R^-T*A^T*b */
            tempValue = 0.0;
            for (int j = N - 1; j > i; --j)
                tempValue += R[i * N + j] * Solution[j];
            Solution[i] = (residuum[i] - tempValue) / R[i * (N + 1)];
        }
        //     cout << "b - A*x" << endl;
        for (int i = 0; i < M; ++i) {
            tempValue = 0.0;
            for (int j = 0; j < N; ++j)
                tempValue += Matrix[i * N + j] * Solution[j];
            residuum[i] = rightHandSide[i] - tempValue;
        }
        //     cout << "A^T*r" << endl;
        for (int i = 0; i < N; ++i) {
            tempValue = 0.0;
            for (int j = 0; j < M; ++j)
                tempValue += Matrix[j * N + i] * residuum[j];
            tempVector[i] = tempValue;
        }
        //     cout << "R^-TA^T*r" << endl;
        for (int i = 0; i < N; ++i) {
            tempValue = 0.0;
            for (int j = 0; j < i; ++j)
                tempValue += R[j * N + i] * residuum[j];
            residuum[i] = (tempVector[i] - tempValue) / R[i * (N + 1)];
        }
        //     cout << "R^-1R^-TA^T*r" << endl;
        for (int i = N - 1; i >= 0; --i) {
            tempValue = 0.0;
            for (int j = N - 1; j > i; --j)
                tempValue += R[i * N + j] * tempVector[j];
            tempVector[i] = (residuum[i] - tempValue) / R[i * (N + 1)];
            Solution[i] += tempVector[i];
        }

        //     for (int i = 0; i < N; ++i)
        //        cout << Solution[i] << endl;
        //     cout << endl;

        delete[] residuum;
        delete[] tempVector;
        delete[] R;

    }

}