SchKGHyp.h

#ifndef gridripper_phys_gr_fixmp_schkghyp_SchKGHyp_h
#define gridripper_phys_gr_fixmp_schkghyp_SchKGHyp_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 schkghyp {


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





class SchKGHyp: 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=1  

    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
  
    int Yl; //  modus l indexe
  
    GReal_t M;          // <- mass
//    GReal_t e;                // <- charge
   // GReal_t K;          // koordinatarendszer K parametere 
    
//    GReal_t q;

    GReal_t minR;
    GReal_t maxR;
    
    GReal_t dhoriz;
    
    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;  
    GReal_t M2,M3,M4; 
    GReal_t AA;  
    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; 
    
    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;
    GReal_t *Lapl_f;
    GReal_t *szamlalo;
    

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

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

    int exclude;
    
    
    int j;  // ciklusokhoz 

public:

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

    ~SchKGHyp();

    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& fw);

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


#endif /* gridripper_phys_gr_fixmp_schkghyp_SchKGHyp_h */