KerrHiggsHyp.h

#ifndef gridripper_phys_gr_fixmp_kerrhiggshyp_KerrHiggsHyp_h
#define gridripper_phys_gr_fixmp_kerrhiggshyp_KerrHiggsHyp_h


#include <gridripper/Parameters.h>
#include <gridripper/amr1d/PDE.h>
#include <gridripper/lang/IllegalArgumentException.h>
#include <gridripper/multipole/ScalarFieldMP.h>

#include <time.h> 


namespace gridripper { namespace phys { namespace gr { namespace fixmp {
namespace kerrhiggshyp {


using namespace gridripper;
using namespace gridripper::multipole;
using namespace std;




class KerrHiggsHyp: public PDE
{
public:
    static const GReal_t EPSILON=1e-8;

    /* Field component indices. */
    enum { ind_f=0, ind_fT=1, ind_fR=2 };

    const int maxOrder;
    const int arraySize;

    const int I_f;      // <- ind_f*arraySize
    const int I_fT;     // <- ind_fT*arraySize
    const int I_fR;     // <- ind_fR*arraySize
    static const int NUM_INDEP_COMPS=2;
    static const int NUM_CONSTRAINTS=1;
    static const int NUM_COMPS=NUM_INDEP_COMPS+NUM_CONSTRAINTS;

    static const string COMP_NAMES[NUM_COMPS];
    static const string CONSTRAINT_NAMES[NUM_CONSTRAINTS];
    static const string COORD_LABELS[2];

private:

    tvalarray<GReal_t> df;

    GReal_t M;          // <- mass
    GReal_t A;          // <- angular momentum
//    GReal_t e;                // <- charge
   // GReal_t K;          // koordinatarendszer K parametere XXX
    

//    GReal_t lambda;    // <- self coupling constant
//    GReal_t phi0;      // <- equilibrium field value
    // (mass=2*sqrt(lambda)*phi0)
//    bool selfInteraction;
//    bool allmassive;
    // The field is f=(phi-phi0)*r in the simulation if selfinteration 
    // is on, otherwise f=phi*r.
//    GReal_t q;

    GReal_t minR;
    GReal_t maxR;
    GReal_t minRIntq;
    GReal_t maxRIntq;

    // coefficients appearing in the diff. eq.
    GReal_t c0;
    GReal_t cP;
    GReal_t cR;
    GReal_t cT;
    GReal_t cRP;
    GReal_t cTP;
    GReal_t cTR;
    GReal_t cRR;
    
    GReal_t R2,R3,R4,R5,R6,R7,R8,R9,R10,R11;
    // GReal_t K2,K3,K4;  // XXX
    GReal_t M2,M3,M4; // XXX
    GReal_t A2,A3,A4;  // XXX
    GReal_t AA,BB,aa,bb;  // XXX
    GReal_t c00,c01,c02,c03,c04,c05,c06,c07,c08,c09;
    GReal_t cR0,cR1,cR2,cR3,cR4,cR5,cR6,cR7,cR8,cR9,cR10,cR11;
    GReal_t cT0,cT1,cT2,cT3,cT4,cT5,cT6,cT7,cT8,cT9,cT10;
   // GReal_t cTR0,cTR1,cTR2,cTR3,cTR4,cTR5,cTR6;
    GReal_t AA0,AA1,AA2,AA3,AA4,AA5,AA6,AA7; // XXX
    
    GReal_t Rprev;
    bool teszt;
    int teszt2;
    time_t time1,time2;
    
    mutable ScalarFieldMP<GComplex_t> DEN;
    mutable GComplex_t H2S2DENDEN;
    mutable GComplex_t scale;
    mutable ScalarFieldMP<GComplex_t> scaledDEN;
    
    mutable ScalarFieldMP<GComplex_t> aux;
    
    GReal_t neumannTolerance;

    /* Start time. */
    GReal_t T0;
//    GReal_t R0;

    static GReal_t horizon;
    static bool isMinkowski;
    static int transformedR;

    int exclude;

public:

    KerrHiggsHyp(const Parameters* p, const int maxorder, const int arraysize) throw(IllegalArgumentException&);

    ~KerrHiggsHyp();

    GReal_t getMinX() const;
    GReal_t getMaxX() const;
    GReal_t getPhysicalMinX() const;
    GReal_t getPhysicalMaxX() const;
    bool isXPeriodic() const;

    int eval(GReal_t* dF, GReal_t T, GReal_t R,
             const GReal_t* F, const GReal_t* F_R,
             const Grid& G, int i, Grad& d, GReal_t dT);

    void evalConstrainedComponents(Grid& G, int i, FieldWrapper& v);

    FieldWrapper* createField(int level) const;

    GReal_t energyDensity(GReal_t T, GReal_t R,
                          const tvalarray<GReal_t>& F) const;

    GReal_t energyFlow(GReal_t T, GReal_t R,
                       const tvalarray<GReal_t>& F) const;

    GReal_t angmomDensity(GReal_t T, GReal_t R,
                          const tvalarray<GReal_t>& F) const;

    GReal_t angmomFlow(GReal_t T, GReal_t R,
                       const tvalarray<GReal_t>& F) const;

    bool canBeExtended(int where) const;

    void mirrorLeft(GReal_t rh, FieldWrapper& f) const;

    void mirrorRight(GReal_t rh, FieldWrapper& f) const;

    bool hasExtendedRefinedFieldData(const Grid& G, int i) const;

    class Energy: public DensityQuantity {
    private:
        const KerrHiggsHyp& pde;
    public:
        Energy(const KerrHiggsHyp* p): DensityQuantity(p), pde(*p) { }
        GReal_t density(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                        const {
            return pde.energyDensity(t, rh, F);
        }
        GReal_t flow(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F) const {
            return pde.energyFlow(t, rh, F);
        }
        DensityQuantity* cloneDensityQuantity() const {
            return new Energy(&pde);
        }
    };
    Energy theEnergy;

    class Angmom: public DensityQuantity {
    private:
        const KerrHiggsHyp& pde;
    public:
        Angmom(const KerrHiggsHyp* p): DensityQuantity(p), pde(*p) { }
        GReal_t density(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                        const {
            return pde.angmomDensity(t, rh, F);
        }
        GReal_t flow(GReal_t t, GReal_t rh, const tvalarray<GReal_t>& F)
                     const {
            return pde.angmomFlow(t, rh, F);
        }
        DensityQuantity* cloneDensityQuantity() const {
            return new Angmom(&pde);
        }
    };
    Angmom theAngmom;

    PhysicalQuantityArray getPhysicalQuantities() const;

    class IntqBorder : public MarkerQuantity
    {
        private:
            const GReal_t position;
        public:
            IntqBorder(const PDE* eq, GReal_t positionarg)
             : MarkerQuantity(eq), position(positionarg)
            {  }
            ~IntqBorder() {  }
            std::string getQuantityName() const { return "IntqBorder"; }
            GReal_t eval(GReal_t t);
            MarkerQuantity* cloneMarkerQuantity() const
            { return new IntqBorder(&getPDE(), position); }
    };

 
private:
    static tvector<string> componentNames(int maxOrder);
    static tvector<string> constraintNames(int maxOrder);

};


} } } } }    // namespace gridripper::phys::gr::fixmp::kerrhiggshyp


#endif /* gridripper_phys_gr_fixmp_kerrhiggshyp_KerrHiggsHyp_h */