moment.h

#ifndef FLAME_MOMENT_H
#define FLAME_MOMENT_H

#include <ostream>
#include <limits>
#include <iomanip>
#include <math.h>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

#include "base.h"

#include "constants.h"

// Default sampling frequency [Hz].
# define SampleFreqDefault   80.5e6

struct Particle {
    double IonZ,         
           IonQ,         
           IonEs,        
           IonW,         
           gamma,        
           beta,         
           bg,           
           SampleFreq,   
           SampleLambda, 
           SampleIonK,   
           phis,         
           IonEk;        

    Particle() {
        phis = 0.0;
        // initially spoil
        IonZ = IonQ = IonEs = IonW
        = gamma = beta = bg
        = SampleFreq = SampleLambda
        = SampleIonK = IonEk
        = std::numeric_limits<double>::quiet_NaN();
    }

    void recalc() {
        IonW       = IonEs + IonEk;
        gamma      = (IonEs != 0e0)? IonW/IonEs : 1e0;
        beta       = sqrt(1e0-1e0/(gamma*gamma));
        bg         = (beta != 0e0)? beta*gamma : 1e0;
        SampleLambda = C0/SampleFreq*MtoMM;
        SampleIonK   = 2e0*M_PI/(beta*SampleLambda);
    }

    double Brho() const { return beta*IonW/(C0*IonZ); } 
};

std::ostream& operator<<(std::ostream&, const Particle&);

inline bool operator==(const Particle& lhs, const Particle& rhs)
{
    // compare only independent variables.
    return lhs.IonEk==rhs.IonEk && lhs.IonEs==rhs.IonEs
            && lhs.IonZ==rhs.IonZ && lhs.IonQ==rhs.IonQ
            && lhs.phis==rhs.phis && lhs.SampleFreq==rhs.SampleFreq;
}

inline bool operator!=(const Particle& lhs, const Particle& rhs)
{
    return !(lhs==rhs);
}

inline bool operator<=(const Particle& lhs, const Particle& rhs)
{
    // compare only independent variables.
    return lhs.IonEk==rhs.IonEk && lhs.IonEs==rhs.IonEs
            && lhs.IonZ==rhs.IonZ && lhs.IonQ==rhs.IonQ
            && lhs.SampleFreq==rhs.SampleFreq;
}

struct MomentState : public StateBase
{
    enum {maxsize=7};
    enum param_t {
        PS_X, PS_PX, PS_Y, PS_PY, PS_S, PS_PS,
        PS_QQ // ???
    };

    MomentState(const Config& c);
    virtual ~MomentState();

    typedef boost::numeric::ublas::vector<double,
                    boost::numeric::ublas::bounded_array<double, maxsize>
    > vector_t;

    typedef boost::numeric::ublas::matrix<double,
                    boost::numeric::ublas::row_major,
                    boost::numeric::ublas::bounded_array<double, maxsize*maxsize>
    > matrix_t;

    virtual void assign(const StateBase& other);

    virtual void show(std::ostream& strm, int level) const;

    Particle ref;

    // all three must have the same length, which is the # of charge states
    std::vector<Particle> real;
    std::vector<vector_t> moment0;
    std::vector<matrix_t> moment1;
    std::vector<matrix_t> transmat;

    vector_t moment0_env, moment0_rms;
    matrix_t moment1_env;

    double last_caviphi0;

    virtual bool getArray(unsigned idx, ArrayInfo& Info);

    virtual MomentState* clone() const {
        return new MomentState(*this, clone_tag());
    }

    void recalc() {
        ref.recalc();
        for(size_t i=0; i<real.size(); i++) real[i].recalc();
    }

    void calc_rms();

    inline size_t size() const { return real.size(); } 

protected:
    MomentState(const MomentState& o, clone_tag);
};

struct MomentElementBase : public ElementVoid
{
    typedef MomentState state_t;

    typedef state_t::matrix_t value_t;

    MomentElementBase(const Config& c);
    virtual ~MomentElementBase();

    void get_misalign(const state_t& ST, const Particle& real, value_t& M, value_t& IM) const;

    unsigned get_flag(const Config& c, const std::string& name, const unsigned& def_value);

    virtual void advance(StateBase& s);

    virtual bool check_cache(const state_t& S) const;

    virtual bool check_backward(const state_t& S) const;

    void resize_cache(const state_t& ST);

    virtual void recompute_matrix(state_t& ST);

    virtual void show(std::ostream& strm, int level) const;

    Particle last_ref_in, last_ref_out;
    std::vector<Particle> last_real_in, last_real_out;
    std::vector<value_t> transfer;
    std::vector<value_t> misalign, misalign_inv;

    double dx, dy, pitch, yaw, roll;

    bool skipcache;

    virtual void assign(const ElementVoid *other) =0;

protected:
    // scratch space to avoid temp. allocation in advance()
    // An Element can't be shared between multiple threads
    state_t::matrix_t scratch;
};

#endif // FLAME_MOMENT_H