Merge branch 'release-v0.5.0'

CMakeLists.txt

@@ -4,7 +4,7 @@ if (POLICY CMP0042)
 endif ()
 PROJECT(EVOEMD CXX)
 
-SET(EVOEMD_VERSION "0.2.5")
+SET(EVOEMD_VERSION "0.5.0")
 
 FIND_PACKAGE(GSL 2.1 REQUIRED)
 FIND_PACKAGE(Eigen3 REQUIRED)
@@ -14,9 +14,10 @@ LINK_DIRECTORIES(${CUBADIR})
 SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin)
 SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
 SET(CMAKE_CXX_STANDARD 11)
+SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unknown-warning-option")
 
 ## Collecting Source Files
-FILE(GLOB lib_SOURCES src/*.cpp)
+FILE(GLOB lib_SOURCES src/EvoEMD/*.cpp)
 
 
 # For .so/.dylib library
@@ -29,4 +30,5 @@ SET_TARGET_PROPERTIES(EVOEMD_shared PROPERTIES
 TARGET_LINK_LIBRARIES(EVOEMD_shared ${GSL_LIBRARIES} cuba)
 SET(EVOEMD_LIBRARY EVOEMD_shared)
 
-add_subdirectory(test)
+add_subdirectory(src/Models/SeesawFreezeIn)
+# add_subdirectory(test)

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2020 Yongcheng Wu
+Copyright (c) 2020-2021 Yongcheng Wu
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -1,2 +1,11 @@
 # EvoEMD
-The Evolution in Early Matter Dominated Universe
+
+A framework to evaluate the evolution w/ or w/o an early matter dominated (EMD) era.
+
+## Feature
+
+- Global parameter system
+- Global particle system
+- Global process system
+- Automatically build up the Boltzmann equation according to the user's definition of particle and process (working in progress)
+- Solving the Boltzmann equation with RK to 5th order

include/EvoEMD/CollisionRate.h

@@ -0,0 +1,110 @@
+#ifndef _COLLISION_RATE_H_
+#define _COLLISION_RATE_H_
+
+#include <map>
+
+#include "EvoEMD/ParticleBase.h"
+#include "EvoEMD/RealTypes.h"
+
+namespace EvoEMD {
+struct Process_Amp {
+    // * used to store how many squared diagrams we have for corresponding amplitude calculation
+    unsigned n_diag;
+
+    // * bool determine whether the corresponding REAL is actually exactly zero
+    typedef std::vector<REAL> CTH_RES_FULL;
+    typedef CTH_RES_FULL NUMERATOR_STRUCTURE;
+    typedef int Propagator_ID;
+    typedef std::pair<Propagator_ID, CTH_RES_FULL> PROPAGATOR_STRUCTURE;
+    typedef std::pair<PROPAGATOR_STRUCTURE, PROPAGATOR_STRUCTURE> DENOMINATOR_STRUCTURE;
+
+    typedef std::vector<NUMERATOR_STRUCTURE> AMP_NUM;
+    typedef std::vector<DENOMINATOR_STRUCTURE> AMP_DEN;
+    AMP_NUM amps_numerator;
+    AMP_DEN amps_denominator;
+
+    const NUMERATOR_STRUCTURE &Get_Numerator(int diagram_id = 0) const;
+    const DENOMINATOR_STRUCTURE &Get_Denominator(int diagram_id = 0) const;
+};
+
+class Amplitude {
+    // * Amplitude used to calculate the collision rate;
+    // * As the purpose of this class is to assist the collision rate calculation and solving Boltzmann equation
+    // * The relevant dof is summed insided Get_Amp()
+    // * The dof has two parts:
+    // * 1. The spin dof: the amp is summed over spin dof for all external particles
+    // * 2. The species dof: CP conjugated component is summed as 1 -> 2 3 is added with 1bar -> 2bar 3bar etc.
+    // *                     If CP violating rate is needed, (1->23) - (1bar->2bar3bar) is provided
+    // *                     Flavor dof is also summed
+public:
+    typedef std::vector<Pseudo_Particle *> INITIAL_STATES;
+    typedef std::vector<Pseudo_Particle *> FINAL_STATES;
+
+    Amplitude();
+    virtual ~Amplitude(){};
+    virtual void Update_Amp(REAL sqrt_shat) = 0;
+    virtual const Process_Amp &Get_Amp(REAL sqrt_shat) {
+        Update_Amp(sqrt_shat);
+        return amp_res;
+    };
+    virtual REAL Get_Coeff(REAL T) = 0;
+
+    unsigned N_INITIAL;
+    INITIAL_STATES INITIAL;
+
+    unsigned N_FINAL;
+    FINAL_STATES FINAL;
+
+protected:
+    Process_Amp amp_res;
+};
+
+class Collision_Rate {
+    // For any process x -> y
+    // CP conserving rate means, gamma(x->y) + gamma(xbar -> ybar)
+    // CP violating rate means, gamma(x->y) - gamma(xbar -> ybar)
+protected:
+    // * This extra factor should be added for particular process,
+    // * if one omit the internal symmetry for simplicity during the calculation
+    // * e.g. In SM, L and Phi are SU2 doublet
+    // * In some case, we don't need to manipulate the component
+    // * but just treat them as a whole object.
+    // * In this case, extra factor should be added.
+    double Extra_Factor_for_Internal_Symmetry;
+
+    // * ptr to Amplitude class, but we don't own it.
+    Amplitude *amp;
+
+public:
+    Collision_Rate(Amplitude *amp) { this->amp = amp; }
+    virtual ~Collision_Rate(){};
+
+    virtual REAL Get_Amp_Integrate_over_Phase_Space(REAL sqrt_shat) = 0;
+    virtual REAL Get_Collision_Rate(REAL T) = 0;
+};
+
+class Decay12_Rate : public Collision_Rate {
+    // * For two body decay
+public:
+    Decay12_Rate(Amplitude *amp) : Collision_Rate(amp){};
+    ~Decay12_Rate(){};
+
+    virtual REAL Get_Amp_Integrate_over_Phase_Space(REAL sqrt_shat);
+    virtual REAL Get_Collision_Rate(REAL T);
+};
+
+class Scatter22_Rate : public Collision_Rate {
+protected:
+    REAL Get_Amp_Integrate_over_Phase_Space_Single_Channel(const Process_Amp::NUMERATOR_STRUCTURE &numerator,
+                                                           const Process_Amp::DENOMINATOR_STRUCTURE &denominator);
+
+public:
+    Scatter22_Rate(Amplitude *amp) : Collision_Rate(amp){};
+    ~Scatter22_Rate(){};
+
+    virtual REAL Get_Amp_Integrate_over_Phase_Space(REAL sqrt_shat);
+    virtual REAL Get_Collision_Rate(REAL T);
+};
+
+}  // namespace EvoEMD
+#endif  //_COLLISION_RATE_H_

include/EvoEMD/Constants.h

@@ -0,0 +1,9 @@
+#ifndef _CONSTANTS_H_
+#define _CONSTANTS_H_
+
+#define PHY_MP 2.435e18
+#define GeV 1
+#define TeV 1000
+#define eV 1e-9
+
+#endif  //_CONSTANTS_H_

include/EffDOF.h → include/EvoEMD/EffDOF.h

@@ -1,7 +1,9 @@
-#ifndef EFFDOF_H
-#define EFFDOF_H
+#ifndef _EFFDOF_H_
+#define _EFFDOF_H_
 
-#include "RealTypes.h"
+#include "EvoEMD/RealTypes.h"
+
+namespace EvoEMD {
 
 /**
  * @brief  Calculate the effective degree of freedom for energy density
@@ -19,4 +21,6 @@ REAL ge(REAL T);
  */
 REAL gs(REAL T);
 
+}  // namespace EvoEMD
+
 #endif  // EFFDOF_H

include/EvoEMD/Hubble.h

@@ -0,0 +1,104 @@
+#ifndef _HUBBLE_H_
+#define _HUBBLE_H_
+
+#include "EvoEMD/RealTypes.h"
+
+namespace EvoEMD {
+
+class Hubble_For_Single_Period {
+protected:
+    REAL T_start;
+    REAL T_end;
+    bool Isentropic;
+    double beta_T;
+    double beta_s;
+
+public:
+    /**
+     * @brief The contructor, default is for a RD or EMD period
+     * @note
+     * @param  T_start: the starting temperature of this period
+     * @param  T_end: the ending temperature of this period
+     * @param  Isentropic: whether this period is a isentropic period.
+     * @param  beta_T: T~a^{-beta_T}
+     * @param  beta_s: s~a^{-beta_s} S=s a^3
+     * @retval
+     */
+    Hubble_For_Single_Period(const REAL T_start, const REAL T_end, const bool Isentropic = true,
+                             const double beta_T = 1, const double beta_s = 3);
+    virtual ~Hubble_For_Single_Period(){};
+
+    /**
+     * @brief Calculate the Hubble parameter at temperature T, need to be overwrite by derived class
+     * @note
+     * @param  T: the temperature
+     * @retval The hubble parameter
+     */
+    virtual REAL Get_Hubble_at_T(const REAL T) = 0;
+    static REAL Get_Hubble_For_RD(const REAL T);
+    REAL Get_T_Start() const { return T_start; }
+    REAL Get_T_End() const { return T_end; }
+    double Get_beta_T() const { return beta_T; }
+    double Get_beta_s() const { return beta_s; }
+    bool Is_Isentropic() const { return Isentropic; }
+};
+
+class Hubble_RD : public Hubble_For_Single_Period {
+public:
+    Hubble_RD(const REAL T_start, const REAL T_end);
+    ~Hubble_RD(){};
+
+    virtual REAL Get_Hubble_at_T(const REAL T) override;
+};
+
+class Hubble_EMD : public Hubble_For_Single_Period {
+private:
+    REAL HRD_at_T_start;
+
+public:
+    Hubble_EMD(const REAL T_start, const REAL T_end);
+    ~Hubble_EMD(){};
+
+    virtual REAL Get_Hubble_at_T(const REAL T) override;
+};
+
+class Hubble_EP : public Hubble_For_Single_Period {
+private:
+    REAL HRD_at_T_end;
+    REAL ge_at_T_end;
+
+public:
+    Hubble_EP(const REAL T_start, const REAL T_end);
+    ~Hubble_EP(){};
+
+    virtual REAL Get_Hubble_at_T(const REAL T) override;
+};
+
+class Hubble_History {
+private:
+    std::vector<Hubble_For_Single_Period *> Periods;
+    std::vector<REAL> Temperatures;
+    REAL TRH;
+    REAL Ti;
+    REAL Te;
+    REAL Tr;
+    REAL Tf;
+
+    void Solve_Te();
+
+public:
+    Hubble_History(const REAL Ti, const REAL Tr);
+    Hubble_History(const Hubble_History &HH);
+    ~Hubble_History();
+
+    Hubble_History &operator=(const Hubble_History &HH);
+
+    int Get_N_Period() const { return Periods.size(); }
+    int Get_Period_ID_at_T(const REAL T) const;
+    REAL Get_Hubble_at_T(const REAL T);
+    Hubble_For_Single_Period *operator[](const int pid);
+};
+
+}  // namespace EvoEMD
+
+#endif  //_HUBBLE_H_

include/EvoEMD/ParameterBase.h

@@ -0,0 +1,97 @@
+#ifndef _PARAMETER_BASE_H_
+#define _PARAMETER_BASE_H_
+
+#include <map>
+#include <set>
+#include <string>
+
+#include "EvoEMD/RealTypes.h"
+
+namespace EvoEMD {
+
+class Parameter_Base {
+protected:
+    std::string name;
+    REAL value;
+    bool updated;
+    std::set<Parameter_Base *> underlying_parameters;  // Parameter in this set depends on current parameter
+
+public:
+    Parameter_Base(std::string par_name) : name(par_name), value(0), updated(true){};
+    virtual ~Parameter_Base(){};
+
+    bool Is_Independent() { return (underlying_parameters.size() == 0); }
+    void Claim_Dependence(Parameter_Base *par) { underlying_parameters.insert(par); }
+    void Notify() {
+        updated = false;
+        for (auto &&bp : underlying_parameters) {
+            bp->Notify();
+        }
+    }
+
+    // * Any parameter derived from this Base_Parameter should re-implement this function
+    virtual void Update_Value(REAL input) = 0;
+
+    void Set_Value(REAL input = 0) {
+        Notify();
+        Update_Value(input);
+        updated = true;
+    }
+    REAL Get_Value() {
+        if (updated) {
+            return value;
+        } else {
+            Set_Value();
+            return value;
+        }
+    }
+    std::string Get_Name() const { return name; }
+};
+
+class Free_Parameter : public Parameter_Base {
+public:
+    Free_Parameter(std::string name, REAL default_value = 0) : Parameter_Base(name) {
+        value = default_value;
+        updated = true;
+    }
+    ~Free_Parameter(){};
+
+    virtual void Update_Value(REAL input) override { value = input; }
+};
+
+class Parameter_Factory {
+private:
+    Parameter_Factory();
+    ~Parameter_Factory();
+    typedef std::map<std::string, Parameter_Base *> Parameter_List;
+    Parameter_List PL;
+
+public:
+    static Parameter_Factory &Get_Parameter_Factory();
+
+    void Register_Parameter(Parameter_Base *par);
+    bool Set_Parameter_Value(std::string name, REAL value);
+    REAL Get_Parameter_Value(std::string name, REAL default_value = 0);
+    Parameter_Base *Get_Parameter(std::string name);
+    void List_Parameters();
+};
+}  // namespace EvoEMD
+
+class Register_Parameter {
+public:
+    Register_Parameter(EvoEMD::Parameter_Base *par) {
+        EvoEMD::Parameter_Factory::Get_Parameter_Factory().Register_Parameter(par);
+    };
+};
+
+#define REGISTER_PARAMETER(paramName) Register_Parameter g_register_parameter_##paramName(new paramName)
+#define DECLARE_FREE_PARAMETER(paramName, value)                   \
+    class param_##paramName : public Free_Parameter {              \
+    public:                                                        \
+        param_##paramName() : Free_Parameter(#paramName, value){}; \
+    };                                                             \
+    REGISTER_PARAMETER(param_##paramName)
+
+#define RETRIVE_PARAMETER(paramName) EvoEMD::Parameter_Factory::Get_Parameter_Factory().Get_Parameter(#paramName)
+
+#endif  //_PARAMETER_BASE_H_