OutputPlane.h

Go to the documentation of this file.

// Copyright (c) 2023, Chris Rogers
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL 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.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//

#ifndef CLASSIC_OutputPlane_HH
#define CLASSIC_OutputPlane_HH

#include "AbsBeamline/PluginElement.h"
#include "BeamlineGeometry/StraightGeometry.h"
#include "Fields/NullField.h"
#include "Steppers/Steppers.h"
#include "Steppers/RK4.h"

#include <string>

class Component;

class OutputPlane: public PluginElement {

public:
    enum class algorithm {RK4STEP, INTERPOLATION};

    explicit OutputPlane(const std::string& name);

    OutputPlane();

    OutputPlane(const OutputPlane&);

    virtual ~OutputPlane();

    ElementBase* clone() const override;
    
    virtual void accept(BeamlineVisitor&) const override;

    inline Component* getGlobalFieldMap() const;
    inline void setGlobalFieldMap(Component* field);

    inline Vector_t getNormal() const;
    inline void setNormal(Vector_t normal);

    inline Vector_t getCentre() const;
    inline void setCentre(Vector_t centre);

    inline double getTolerance() const;
    inline void setTolerance(double tolerance);

    inline size_t getMaxIterations() const;
    inline void setMaxIterations(size_t max);

    inline int getRecentre() const;
    inline void setRecentre(int willRecentre);

    inline double getRadialExtent() const;
    inline void setRadialExtent(double r);

    inline double getVerticalExtent() const;
    inline void setVerticalExtent(double z);

    inline double getHorizontalExtent() const;
    inline void setHorizontalExtent(double width);

    inline algorithm getAlgorithm() const;
    inline void setAlgorithm(algorithm alg);

    inline int getVerboseLevel() const;
    inline void setVerboseLevel(int verbose);

    bool checkOne(const int index, const double tstep, double chargeToMass,
                  double& t, Vector_t& R, Vector_t& P);

    NullField& getField() override {return nullfield_m;}
    const NullField& getField() const override {return nullfield_m;}

    StraightGeometry& getGeometry() override {return geom_m;}
    const StraightGeometry& getGeometry() const override {return geom_m;}
    void RK4Step(const double& tstep,
                 const double& chargeToMass,
                 const double& t,
                 Vector_t& R,
                 Vector_t& P) const;

    void operator=(const OutputPlane&) = delete;
    
    //ElementBase::ElementType getType() const;
    ElementType getType() const override;
    void recentre(Vector_t R, Vector_t P);


private:
    virtual void doInitialise(PartBunchBase<double, 3>* /*bunch*/) override;

    inline bool doPreCheck(PartBunchBase<double, 3> *bunch) override;

    virtual bool doCheck(PartBunchBase<double, 3> *bunch,
                         const int turnnumber, 
                         const double t, 
                         const double tstep) override;
    virtual void doGoOffline() override;
    void getDerivatives(const Vector_t& R,
                        const Vector_t& P,
                        const double& t,
                        const double& chargeToMass,
                        double* yp) const;
    double distanceToPlane(Vector_t point) const;
    void rk4Test(double tstep, double chargeToMass,
                 double& t, Vector_t& R, Vector_t& P);
    void interpolation(double& t, Vector_t& R, Vector_t& P);
    bool getFieldsAtPoint(const Vector_t& R, const Vector_t& P, const double& t, Vector_t& Efield, Vector_t& Bfield);
    
    Component* field_m = NULL; // field map - this is a borrowed pointer
    Vector_t normal_m; // normal to the output plane
    Vector_t centre_m; // centre of the output plane
    double radialExtent_m = -1.0; // maximum radial extent of the plane (circular)
    double verticalExtent_m = -1.0; // maximum vertical extent of the plane (rectangular)
    double horizontalExtent_m = -1.0; // maximum horizontal extent of the plane (rectangular)
    size_t maxIterations_m = 10; // maximum number of iterations when finding intercept
    double tolerance_m = 1e-9; // tolerance on distance from plane when finding intercept
    NullField nullfield_m; // dummy variable for inheritance
    StraightGeometry geom_m; // dummy variable for inheritance
    algorithm algorithm_m; // sets interpolation or RK4
    int recentre_m = -1; // particle index for recentring
    int nHits_m = 0; // counter for number of hits on the plane
    int verbose_m = 0; // verbosity

    typedef std::function<bool(const double&,
                               const size_t&,
                               Vector_t&,
                               Vector_t&)> function_t;
    std::unique_ptr<Stepper<function_t> > stepper_m;
};

Component* OutputPlane::getGlobalFieldMap() const {
    return field_m;
}

void OutputPlane::setGlobalFieldMap(Component* field) {
    field_m = field;
}

Vector_t OutputPlane::getNormal() const {
    return normal_m;
}

void OutputPlane::setNormal(Vector_t normal) {
    normal_m = normal;
    normal_m /= euclidean_norm(normal);
}

Vector_t OutputPlane::getCentre() const {
    return centre_m;
}

void OutputPlane::setCentre(Vector_t centre) {
    centre_m = centre;
}

double OutputPlane::getTolerance() const {
    return tolerance_m;
}

void OutputPlane::setTolerance(double tolerance) {
    tolerance_m = tolerance;
}

size_t OutputPlane::getMaxIterations() const {
    return maxIterations_m;
}

void OutputPlane::setMaxIterations(size_t max) {
    maxIterations_m = max;
}

double OutputPlane::getHorizontalExtent() const {
    return horizontalExtent_m;
}

void OutputPlane::setHorizontalExtent(double width) {
    horizontalExtent_m = width;
}

double OutputPlane::getVerticalExtent() const {
    return verticalExtent_m;
}

void OutputPlane::setVerticalExtent(double width) {
    verticalExtent_m = width;
}

double OutputPlane::getRadialExtent() const {
    return radialExtent_m;
}

void OutputPlane::setRadialExtent(double radius) {
    radialExtent_m = radius;
}

int OutputPlane::getRecentre() const {
    return recentre_m;
}

void OutputPlane::setRecentre(int recentre) {
    recentre_m = recentre;
}

OutputPlane::algorithm OutputPlane::getAlgorithm() const {
    return algorithm_m;
}

void OutputPlane::setAlgorithm(OutputPlane::algorithm alg) {
    algorithm_m = alg;
}

bool OutputPlane::doPreCheck(PartBunchBase<double, 3>* /*bunch*/) {
    return true;
}

int OutputPlane::getVerboseLevel() const {
    return verbose_m;
}

void OutputPlane::setVerboseLevel(int verbose) {
    verbose_m = verbose;
}


#endif // CLASSIC_OutputPlane_HH