EnvelopeBunch.h

Go to the documentation of this file.

#ifndef _ENVELOPE_BUNCH_H
#define _ENVELOPE_BUNCH_H

#include "Algorithms/bet/profile.h"
#include "Algorithms/bet/EnvelopeSlice.h"
#include "Utilities/OpalException.h"
#include "Algorithms/PartBunch.h"
#include <Physics/Physics.h>

#include <assert.h>
#include <vector>

#include <memory>

enum EnvelopeBunchParameter {
    sp_beta,      
    sp_gamma,     
    sp_z,         
    sp_I,         
    sp_Rx,        
    sp_Ry,        
    sp_Px,        
    sp_Py,        
    sp_Pz,        
    sp_x0,        
    sp_y0,        
    sp_px0,       
    sp_py0        
};

enum SolverParameters {
    sv_fixedStep    = 0x0001,  
    sv_fieldOutside = 0x0002,  
    sv_radial       = 0x0004,  
    sv_offaxis      = 0x0008,  
    sv_lwakes       = 0x0010,  
    sv_twakes       = 0x0020,  
    sv_s_path       = 0x0100   
};

enum DataStatus {
    ds_fieldsSynchronized = 0x0001,    
    ds_slicesSynchronized = 0x0002,    
    ds_currentCalculated  = 0x0004,    
    ds_wakesZCalculated   = 0x0008,    
    ds_wakesXYCalculated  = 0x0010,    
    ds_spaceCharge        = 0x001c     
};

enum EnvelopeBunchShape {
    bsRect,
    bsGauss
};


class EnvelopeBunch : public PartBunch {

public:
    EnvelopeBunch(const PartData *ref);

    EnvelopeBunch(const std::vector<OpalParticle> &,
                  const PartData *ref);

    EnvelopeBunch(const EnvelopeBunch &) = delete;

    virtual ~EnvelopeBunch();

    void createBunch();

    void distributeSlices(int nSlice = 101);

    bool isValid_m;

    std::unique_ptr<Profile> currentProfile_m;

    IpplTimings::TimerRef calcITimer_m;
    IpplTimings::TimerRef spaceChargeTimer_m;

    using PartBunch::initialize;
    void initialize(int sli, double charge, double energy, double width, double te, double frac,
                    double current, double center, double bX, double bY, double mX, double mY, double Bz, int nbin);

    bool isRadial() { return solver_m & sv_radial; }

    bool isOffaxis() { return solver_m & sv_offaxis; }

    void calcBeamParameters();

    void computeSpaceCharge();

    void timeStep(double tStep, double zCat = 0.0);

    void derivs(double tc, double Y[], double dYdt[]);

    double Eavg();
    double zAvg();
    double zTail();
    double zHead();
    double time() { return t_m; }

    void setEx(double emi) { emtnx0_m = emi; }
    void setEy(double emi) { emtny0_m = emi; }
    void setSolverParameter(int s) { solver_m = s; }
    // set particle energy of bunch in [eV] and optional the correlated energy spread [eV/m]
    void setEnergy(double, double = 0.0);

    void setExternalFields(int i, const Vector_t &EF, const Vector_t &BF, const Vector_t &KR, const Vector_t &KT) {
        this->EF[i] = EF;
        this->BF[i] = BF;
        this->KR[i] = KR;
        this->KT[i] = KT;
    }

    void getExternalFields(int i, Vector_t &EF, Vector_t &BF, Vector_t &KR, Vector_t &KT) const {
        EF = this->EF[i];
        BF = this->BF[i];
        KR = this->KR[i];
        KT = this->KT[i];
    }

    double get_sPos();
    double AvBField();
    double AvEField();
    int getLocalNum() { return numMySlices_m; }
    int getTotalNum() { return numSlices_m; }
    double getT() { return t_m; }
    double get_meanKineticEnergy() { return Eavg(); }
    double get_dEdt() { return dEdt_m; }
    Vector_t sigmax() { return sigmax_m; }
    Vector_t sigmap() { return sigmap_m; }
    Vector_t emtn() { return emtn_m; }
    Vector_t get_norm_emit() { return norm_emtn_m; }
    Vector_t maxX() { return maxX_m; }
    Vector_t minX() { return minX_m; }
    Vector_t maxP() { return maxP_m; }
    Vector_t minP() { return minP_m; }
    double getChargePerParticle() { return Q_m / numSlices_m; }
    Vector_t get_rrms() { return sigmax_m; }
    Vector_t get_prms() { return sigmap_m; }

    size_t mySliceStartOffset() { return mySliceStartOffset_m; }
    size_t mySliceEndOffset() { return mySliceEndOffset_m; }
    size_t numMySlices() { return numMySlices_m; }

    Inform &slprint(Inform &os);

    void setCharge(double _Q) {
        sign_m = _Q < 0.0 ? -1 : 1;
        Q_m = std::abs(_Q);
    }

    double getGamma(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->computeGamma();
    }

    double getBeta(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_beta];
    }
    void setBeta(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_beta] = val;
    }

    void setR(int i, const Vector_t &R) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_x] = R[0];
        slices_m[i]->p[SLI_y] = R[1];
        slices_m[i]->p[SLI_z] = R[2];
    }

    Vector_t getR(int i) {
        assert(i < numMySlices_m);
        Vector_t R;
        R[0] = slices_m[i]->p[SLI_x];
        R[1] = slices_m[i]->p[SLI_y];
        R[2] = slices_m[i]->p[SLI_z];

        return R;
    }

    void setP(int i, const Vector_t &P) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_px] = P[0];
        slices_m[i]->p[SLI_py] = P[1];
        slices_m[i]->p[SLI_beta] = P[2] / sqrt(P[2]*P[2] + 1.0);
    }

    Vector_t getP(int i) {
        assert(i < numMySlices_m);
        Vector_t P;
        P[0] = slices_m[i]->p[SLI_px];
        P[1] = slices_m[i]->p[SLI_py];
        P[2] = slices_m[i]->p[SLI_beta] * slices_m[i]->computeGamma();

        return P;
    }

    void setZ(int i, double coo) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_z] = coo;
    }

    double getZ(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_z];
    }

    double getX(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_x];
    }
    void setX(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_x] = val;
    }

    double getY(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_y];
    }
    void setY(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_y] = val;
    }

    double getX0(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_x0];
    }
    void setX0(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_x0] = val;
    }

    double getY0(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_y0];
    }
    void setY0(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_y0] = val;
    }

    double getPx(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_px];
    }
    void setPx(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_px] = val;
    }

    double getPy(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_py];
    }
    void setPy(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_py] = val;
    }

    double getPz(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_beta] * Physics::m_e * slices_m[i]->computeGamma();
    }

    double getPx0(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_px0];
    }
    void setPx0(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_px0] = val;
    }

    double getPy0(int i) {
        assert(i < numMySlices_m);
        return slices_m[i]->p[SLI_py0];
    }
    void setPy0(int i, double val) {
        assert(i < numMySlices_m);
        slices_m[i]->p[SLI_py0] = val;
    }

    void get_bounds(Vector_t &min, Vector_t &max) {
        min[0] = 0.0;
        min[1] = 0.0;
        min[2] = this->zTail();
        max[0] = 0.0;
        max[1] = 0.0;
        max[2] = this->zHead();
    }


private:
    int numSlices_m;
    int numMySlices_m;
    size_t mySliceStartOffset_m;
    size_t mySliceEndOffset_m;
    std::vector<double> z_m;
    std::vector<double> b_m;
    std::vector< std::vector<int> > bins_m;
    double hbin_m;
    int nebin_m;
    int firstBinWithValue_m;
    int activeSlices_m;

    int solver_m;
    int dStat_m;
    double t_m;
    double t_offset_m;
    double emtnx0_m, emtny0_m;
    double emtbx0_m, emtby0_m;
    double dx0_m, dy0_m;
    double dfi_x_m, dfi_y_m;
    double Bz0_m;
    double Q_m;
    double I0avg_m;
    Vector_t Esl_m;
    Vector_t Bsl_m;
    Vector_t KRsl_m;
    Vector_t KTsl_m;
    std::unique_ptr<Vector_t[]> KR;
    std::unique_ptr<Vector_t[]> KT;
    std::unique_ptr<Vector_t[]> EF;
    std::unique_ptr<Vector_t[]> BF;
    std::vector< std::shared_ptr<EnvelopeSlice> > slices_m;
    unsigned int activeSlice_m;
    int sign_m;
    int currentSlice_m;
    double zCat_m;
    std::vector<double> Exw_m;
    std::vector<double> Eyw_m;
    std::vector<double> Ezw_m;
    std::vector<double> Esct_m;
    std::vector<double> G_m;

    int nValid_m;
    double z0_m;
    double emission_time_step_m;
    size_t lastEmittedBin_m;
    double E_m;
    double dEdt_m;
    double Einc_m;
    double tau_m;
    double I_m;
    double Irms_m;
    double Rx_m;
    double Ry_m;
    double RxMax_m;
    double RyMax_m;
    double RxMin_m;
    double RyMin_m;
    double Px_m;
    double Py_m;
    double x0_m;
    double y0_m;
    double x0Max_m;
    double y0Max_m;
    double x0Min_m;
    double y0Min_m;
    // double dx02_m;
    // double dy02_m;
    // double dfi_x_m;
    // double dfi_y_m;
    double Ez_m;
    double Bz_m;
    Vector_t maxX_m;
    Vector_t minX_m;
    Vector_t maxP_m;
    Vector_t minP_m;
    Vector_t sigmax_m;
    Vector_t sigmap_m;
    Vector_t emtn_m;
    Vector_t norm_emtn_m;

    // used in derivs() to remove call to zHead/Tail() containing MPI
    // collectives causing problems if slices are unevenly distributed across
    // processors.
    double curZHead_m;
    double curZTail_m;

    void synchronizeSlices();

    void runStats(EnvelopeBunchParameter sp, double *xAvg, double *xMax, double *xMin, double *rms, int *nValid);

    void calcEmittance(double *emtnx, double *emtny, double *emtx, double *emty, int *nValid);

    void calcEnergyChirp(double *g0, double *dgdt, double *gInc, int *nValid);

    void setBinnedLShape(EnvelopeBunchShape shape, double z0, double w, double frac);

    void setTShape(double enx, double eny, double rx, double ry, double b0);

    void setTOffset(double x0, double px0, double y0, double py0);

    double moveZ0(double zC);

    void backup() {
        for (auto & slice : slices_m)
            slice->backup();
    }


    double tReset(double dt = 0.0);

    void cSpaceCharge();

    void calcI();
};

inline Inform &operator<<(Inform &os, EnvelopeBunch &p) {
    return p.slprint(os);
}
#endif