KerrKGaxHyp.h

#ifndef gridripper_phys_gr_fixmp_kerrkgaxhyp_KerrKGaxHyp_h
#define gridripper_phys_gr_fixmp_kerrkgaxhyp_KerrKGaxHyp_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 kerrkgaxhyp {


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





class KerrKGaxHyp: 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;  // arraySize=2*maxOrder


    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> gradf;  // evalConstrainedComponents-hez
    
    // Result vector (time derivative). 
   // tvalarray<GReal_t> df;

   int Ym,Yl,absYm,parity;
   
   
    GReal_t M;          // <- mass
    GReal_t A;          // <- angular momentum
//    GReal_t e;                // <- charge
   // GReal_t K;          // koordinatarendszer K parametere XXX
    


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

    int maxi2,maxLevel2,gridsize,gri;
    
    // 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 R,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11;
    // GReal_t K2,K3,K4;  // XXX
    GReal_t M2,M3,M4; 
    GReal_t A2,A3,A4;  
    GReal_t AA,BB,aa,bb; 
    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;
    
    const GReal_t *f,*fT,*fR,*pf,*pfT,*pfR;
    GReal_t *df,*dfT,*dfR;
    GReal_t *p_phi_f, *p_phi_fR, *p_phi_fT, *Lapl_f;
    GReal_t *szamlalo;
    
    GReal_t *aux,*aux2;
    
    GReal_t *scale_pre,*aa_pre,*bb_pre,*nBN_pre;
    
    //GReal_t nBN;  // normBoundNeumann
    
    ScalarFieldMP<GComplex_t> DEN;
    GReal_t H2S2DENDEN;
    //mutable GComplex_t scale;
    //mutable ScalarFieldMP<GComplex_t> scaledDEN;
    
    //mutable ScalarFieldMP<GComplex_t> aux;
    
    GReal_t neumannTolerance;

    /* Start time. */
    GReal_t T0;
    
    //Start coordinate of light ray. 
//    GReal_t R0;

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

    int exclude;
    
    int j,j2; // ciklusokhoz
    
    typedef GReal_t vekt3[3]; 
    vekt3 *GC;  //Gaunt koefficienseknek
    
   
    
    void mulY20(GReal_t *v, int numofcomps, GReal_t *result);
    
    void NeumannDiv20(GReal_t *f, GReal_t a0, GReal_t a2, GReal_t *result, int numcomp, GReal_t norm, const GReal_t tolerance);
    

public:

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

    ~KerrKGaxHyp();

    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 KerrKGaxHyp& pde;
    public:
        Energy(const KerrKGaxHyp* 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 KerrKGaxHyp& pde;
    public:
        Angmom(const KerrKGaxHyp* 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::kerrkgaxhyp


#endif /* gridripper_phys_gr_fixmp_kerrkgaxhyp_KerrKGaxHyp_h */