08p1_RPH.f90

Go to the documentation of this file.

!----------------------------------------------------------------------------------------
!
! This file is part of EFTCAMB.
!
! Copyright (C) 2013-2016 by the EFTCAMB authors
!
! The EFTCAMB code is free software;
! You can use it, redistribute it, and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation;
! either version 3 of the License, or (at your option) any later version.
! The full text of the license can be found in the file eftcamb/LICENSE at
! the top level of the EFTCAMB distribution.
!
!----------------------------------------------------------------------------------------


!----------------------------------------------------------------------------------------


module eftcamb_reparametrized_horndeski

    use precision
    use inifile
    use amlutils
    use eftcamb_cache
    use eft_def
    use eftcamb_abstract_parametrizations_1d
    use eftcamb_neutral_parametrization_1d
    use eftcamb_constant_parametrization_1d
    use eftcamb_linear_parametrizations_1d
    use eftcamb_power_law_parametrizations_1d
    use eftcamb_exponential_parametrizations_1d
    use eftcamb_cpl_parametrizations_1d
    use eftcamb_jbp_parametrizations_1d
    use eftcamb_turning_point_parametrizations_1d
    use eftcamb_taylor_parametrizations_1d
    use eftcamb_abstract_model_designer

    implicit none

    private

    public eftcamb_rph

    !----------------------------------------------------------------------------------------
    type, extends ( eftcamb_designer_model ) :: eftcamb_rph

        ! the RPH functions model selection flags:
        integer  :: rphwde
        integer  :: rphmasspmodel
        integer  :: rphkineticitymodel
        integer  :: rphbraidingmodel
        integer  :: rphtensormodel

        ! the RPH functions:
        class( parametrized_function_1d ), allocatable :: rph_wde
        class( parametrized_function_1d ), allocatable :: rph_planckmass
        class( parametrized_function_1d ), allocatable :: rph_kineticity
        class( parametrized_function_1d ), allocatable :: rph_braiding
        class( parametrized_function_1d ), allocatable :: rph_tensor

    contains

        ! initialization of the model:
        procedure :: read_model_selection            => eftcambrphreadmodelselectionfromfile  
        procedure :: allocate_model_selection        => eftcambrphallocatemodelselection      
        procedure :: init_model_parameters           => eftcambrphinitmodelparameters         
        procedure :: init_model_parameters_from_file => eftcambrphinitmodelparametersfromfile 
        ! utility functions:
        procedure :: compute_param_number  => eftcambrphcomputeparametersnumber    
        procedure :: feedback              => eftcambrphfeedback                   
        procedure :: parameter_names       => eftcambrphparameternames             
        procedure :: parameter_names_latex => eftcambrphparameternameslatex        
        procedure :: parameter_values      => eftcambrphparametervalues            
        ! CAMB related procedures:
        procedure :: compute_background_eft_functions  => eftcambrphbackgroundeftfunctions   
        procedure :: compute_secondorder_eft_functions => eftcambrphsecondordereftfunctions  
        procedure :: compute_dtauda                    => eftcambrphcomputedtauda            
        procedure :: compute_adotoa                    => eftcambrphcomputeadotoa            
        procedure :: compute_h_derivs                  => eftcambrphcomputehubbleder         
        ! stability procedures:
        procedure :: additional_model_stability        => eftcambrphadditionalmodelstability 

    end type eftcamb_rph

    ! ---------------------------------------------------------------------------------------------

contains

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphreadmodelselectionfromfile( self, Ini )

        implicit none

        class(eftcamb_rph) :: self
        type(tinifile)              :: ini

        ! read model selection flags:
        self%RPHwDE             = ini_read_int_file( ini, 'RPHwDE'            , 0 )
        self%RPHmassPmodel      = ini_read_int_file( ini, 'RPHmassPmodel'     , 0 )
        self%RPHkineticitymodel = ini_read_int_file( ini, 'RPHkineticitymodel', 0 )
        self%RPHbraidingmodel   = ini_read_int_file( ini, 'RPHbraidingmodel'  , 0 )
        self%RPHtensormodel     = ini_read_int_file( ini, 'RPHtensormodel'    , 0 )

    end subroutine eftcambrphreadmodelselectionfromfile

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphallocatemodelselection( self )

        implicit none

        class(eftcamb_rph)                       :: self

        ! allocate wDE:
        if ( allocated(self%RPH_wDE) ) deallocate(self%RPH_wDE)
        select case ( self%RPHwDE )
            case(0)
                allocate( wde_lcdm_parametrization_1d::self%RPH_wDE )
            case(1)
                allocate( constant_parametrization_1d::self%RPH_wDE )
            case(2)
                allocate( cpl_parametrization_1d::self%RPH_wDE )
                call self%RPH_wDE%set_param_names( ['EFTw0', 'EFTwa'], ['w_0', 'w_a'] )
            case(3)
                allocate( jbp_parametrization_1d::self%RPH_wDE )
                call self%RPH_wDE%set_param_names( ['EFTw0', 'EFTwa', 'EFTwn'], [ 'w_0', 'w_a', 'n  ' ] )
            case(4)
                allocate( turning_point_parametrization_1d::self%RPH_wDE )
                call self%RPH_wDE%set_param_names( ['EFTw0 ', 'EFTwa ', 'EFTwat'], ['w_0', 'w_a', 'a_t'] )
            case(5)
                allocate( taylor_parametrization_1d::self%RPH_wDE )
                call self%RPH_wDE%set_param_names( ['EFTw0', 'EFTwa', 'EFTw2', 'EFTw3'], ['w_0', 'w_a', 'w_2', 'w_3'] )
            case default
                write(*,'(a,I3)') 'No model corresponding to RPH_wDE =', self%RPHwDE
                write(*,'(a)')    'Please select an appropriate model.'
        end select
        ! allocate RPH_PlanckMass:
        if ( allocated(self%RPH_PlanckMass) ) deallocate(self%RPH_PlanckMass)
        select case ( self%RPHmassPmodel )
            case(0)
                allocate( zero_parametrization_1d::self%RPH_PlanckMass )
            case(1)
                allocate( constant_parametrization_1d::self%RPH_PlanckMass )
            case(2)
                allocate( linear_parametrization_1d::self%RPH_PlanckMass )
            case(3)
                allocate( power_law_parametrization_1d::self%RPH_PlanckMass )
                call self%RPH_PlanckMass%set_param_names( ['RPHmassP0  ', 'RPHmassPExp'], ['\tilde{M}_0                  ', '\tilde{M}_{\rm exp}^{\rm RPH}'] )
            case(4)
                allocate( exponential_parametrization_1d::self%RPH_PlanckMass )
                call self%RPH_PlanckMass%set_param_names( ['RPHmassP0  ', 'RPHmassPExp'], ['\tilde{M}_0                  ', '\tilde{M}_{\rm exp}^{\rm RPH}'] )
            case default
                write(*,'(a,I3)') 'No model corresponding to RPH_PlanckMass =', self%RPHmassPmodel
                write(*,'(a)')    'Please select an appropriate model.'
        end select
        ! allocate RPH_Kineticity:
        if ( allocated(self%RPH_Kineticity) ) deallocate(self%RPH_Kineticity)
        select case ( self%RPHkineticitymodel )
            case(0)
                allocate( zero_parametrization_1d::self%RPH_Kineticity )
            case(1)
                allocate( constant_parametrization_1d::self%RPH_Kineticity )
            case(2)
                allocate( linear_parametrization_1d::self%RPH_Kineticity )
            case(3)
                allocate( power_law_parametrization_1d::self%RPH_Kineticity )
                call self%RPH_Kineticity%set_param_names( ['RPHkineticity0  ', 'RPHkineticityExp'], ['\alpha_0^{\rm K}        ', '\alpha_{\rm exp}^{\rm K}'] )
            case(4)
                allocate( exponential_parametrization_1d::self%RPH_Kineticity )
                call self%RPH_Kineticity%set_param_names( ['RPHkineticity0  ', 'RPHkineticityExp'], ['\alpha_0^{\rm K}        ', '\alpha_{\rm exp}^{\rm K}'] )
            case default
                write(*,'(a,I3)') 'No model corresponding to RPH_Kineticity =', self%RPHkineticitymodel
                write(*,'(a)')    'Please select an appropriate model.'
        end select
        ! allocate RPH_Braiding:
        if ( allocated(self%RPH_Braiding) ) deallocate(self%RPH_Braiding)
        select case ( self%RPHbraidingmodel )
            case(0)
                allocate( zero_parametrization_1d::self%RPH_Braiding )
            case(1)
                allocate( constant_parametrization_1d::self%RPH_Braiding )
            case(2)
                allocate( linear_parametrization_1d::self%RPH_Braiding )
            case(3)
                allocate( power_law_parametrization_1d::self%RPH_Braiding )
                call self%RPH_Braiding%set_param_names( ['RPHbraiding0  ', 'RPHbraidingExp'], ['\alpha_0^{\rm B}        ', '\alpha_{\rm exp}^{\rm B}'] )
            case(4)
                allocate( exponential_parametrization_1d::self%RPH_Braiding )
                call self%RPH_Braiding%set_param_names( ['RPHbraiding0  ', 'RPHbraidingExp'], ['\alpha_0^{\rm B}        ', '\alpha_{\rm exp}^{\rm B}'] )
            case default
                write(*,'(a,I3)') 'No model corresponding to RPH_Braiding =', self%RPHbraidingmodel
                write(*,'(a)')    'Please select an appropriate model.'
        end select
        ! allocate RPH_Braiding:
        if ( allocated(self%RPH_Tensor) ) deallocate(self%RPH_Tensor)
        select case ( self%RPHtensormodel )
            case(0)
                allocate( zero_parametrization_1d::self%RPH_Tensor )
            case(1)
                allocate( constant_parametrization_1d::self%RPH_Tensor )
            case(2)
                allocate( linear_parametrization_1d::self%RPH_Tensor )
            case(3)
                allocate( power_law_parametrization_1d::self%RPH_Tensor )
                call self%RPH_Tensor%set_param_names( ['RPHtensor0  ', 'RPHtensorExp'], ['\alpha_0^{\rm T}        ', '\alpha_{\rm exp}^{\rm T}'] )
            case(4)
                allocate( exponential_parametrization_1d::self%RPH_Tensor )
                call self%RPH_Tensor%set_param_names( ['RPHtensor0  ', 'RPHtensorExp'], ['\alpha_0^{\rm T}        ', '\alpha_{\rm exp}^{\rm T}'] )
            case default
                write(*,'(a,I3)') 'No model corresponding to RPH_Tensor =', self%RPHtensormodel
                write(*,'(a)')    'Please select an appropriate model.'
        end select

        ! initialize the names:
        call self%RPH_wDE%set_name        ( 'RPHw'         , 'w'              )
        call self%RPH_PlanckMass%set_name ( 'RPHmassP'     , '\tilde{M}'      )
        call self%RPH_Kineticity%set_name ( 'RPHkineticity', '\alpha^{\rm K}' )
        call self%RPH_Braiding%set_name   ( 'RPHbraiding'  , '\alpha^{\rm B}' )
        call self%RPH_Tensor%set_name     ( 'RPHtensor'    , '\alpha^{\rm T}' )

    end subroutine eftcambrphallocatemodelselection

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphinitmodelparameters( self, array )

        implicit none

        class(eftcamb_rph)                            :: self
        real(dl), dimension(self%parameter_number), intent(in) :: array

        real(dl), allocatable, dimension(:)                    :: temp
        integer :: num_params_function, num_params_temp, i

        num_params_temp     = 1

        ! first elements are w_DE parameters:
        num_params_function = self%RPH_wDE%parameter_number
        allocate( temp(num_params_function) )
        do i = 1, num_params_function
            temp(i)         = array(num_params_temp)
            num_params_temp = num_params_temp +1
        end do
        call self%RPH_wDE%init_parameters(temp)
        deallocate( temp )
        ! then RPH_PlanckMass parameters:
        num_params_function = self%RPH_PlanckMass%parameter_number
        allocate( temp(num_params_function) )
        do i = 1, num_params_function
            temp(i)         = array(num_params_temp)
            num_params_temp = num_params_temp +1
        end do
        call self%RPH_PlanckMass%init_parameters(temp)
        deallocate(temp)
        ! then RPH_Kineticity parameters:
        num_params_function = self%RPH_Kineticity%parameter_number
        allocate( temp(num_params_function) )
        do i = 1, num_params_function
            temp(i)         = array(num_params_temp)
            num_params_temp = num_params_temp +1
        end do
        call self%RPH_Kineticity%init_parameters(temp)
        deallocate(temp)
        ! then RPH_Braiding parameters:
        num_params_function = self%RPH_Braiding%parameter_number
        allocate( temp(num_params_function) )
        do i = 1, num_params_function
            temp(i)         = array(num_params_temp)
            num_params_temp = num_params_temp +1
        end do
        call self%RPH_Braiding%init_parameters(temp)
        deallocate(temp)
        ! then RPH_Tensor parameters:
        num_params_function = self%RPH_Tensor%parameter_number
        allocate( temp(num_params_function) )
        do i = 1, num_params_function
            temp(i)         = array(num_params_temp)
            num_params_temp = num_params_temp +1
        end do
        call self%RPH_Tensor%init_parameters(temp)
        deallocate(temp)

        ! now check the length of the parameters:
        if ( num_params_temp-1 /= self%parameter_number ) then
            write(*,*) 'In EFTCAMBRPHInitModelParameters:'
            write(*,*) 'Length of num_params_temp and self%parameter_number do not coincide.'
            write(*,*) 'num_params_temp:', num_params_temp-1
            write(*,*) 'self%parameter_number:', self%parameter_number
            call mpistop('EFTCAMB error')
        end if

    end subroutine eftcambrphinitmodelparameters

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphinitmodelparametersfromfile( self, Ini )

        implicit none

        class(eftcamb_rph)          :: self
        type(tinifile)              :: ini

        call self%RPH_wDE%init_from_file( ini )
        call self%RPH_PlanckMass%init_from_file( ini )
        call self%RPH_Kineticity%init_from_file( ini )
        call self%RPH_Braiding%init_from_file( ini )
        call self%RPH_Tensor%init_from_file( ini )

    end subroutine eftcambrphinitmodelparametersfromfile

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphcomputeparametersnumber( self )

        implicit none

        class(eftcamb_rph)  :: self

        self%parameter_number = 0
        self%parameter_number = self%parameter_number +self%RPH_wDE%parameter_number
        self%parameter_number = self%parameter_number +self%RPH_PlanckMass%parameter_number
        self%parameter_number = self%parameter_number +self%RPH_Kineticity%parameter_number
        self%parameter_number = self%parameter_number +self%RPH_Braiding%parameter_number
        self%parameter_number = self%parameter_number +self%RPH_Tensor%parameter_number

    end subroutine eftcambrphcomputeparametersnumber

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphfeedback( self, print_params )

        implicit none

        class(eftcamb_rph)  :: self
        logical, optional   :: print_params

        ! print general model informations:
        write(*,*)
        write(*,'(a,a)')    '   Model               =  ', self%name
        write(*,'(a,I3)')   '   Number of params    ='  , self%parameter_number

        ! print model functions informations:
        write(*,*)
        if ( self%RPHwDE             /= 0 ) write(*,'(a,I3)') '   RPHwDE             =', self%RPHwDE
        if ( self%RPHmassPmodel      /= 0 ) write(*,'(a,I3)') '   RPHmassPmodel      =', self%RPHmassPmodel
        if ( self%RPHkineticitymodel /= 0 ) write(*,'(a,I3)') '   RPHkineticitymodel =', self%RPHkineticitymodel
        if ( self%RPHbraidingmodel   /= 0 ) write(*,'(a,I3)') '   RPHbraidingmodel   =', self%RPHbraidingmodel
        if ( self%RPHtensormodel     /= 0 ) write(*,'(a,I3)') '   RPHtensormodel     =', self%RPHtensormodel

        write(*,*)
        ! print functions informations:
        call self%RPH_wDE%feedback        ( print_params )
        call self%RPH_PlanckMass%feedback ( print_params )
        call self%RPH_Kineticity%feedback ( print_params )
        call self%RPH_Braiding%feedback   ( print_params )
        call self%RPH_Tensor%feedback     ( print_params )

    end subroutine eftcambrphfeedback

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphparameternames( self, i, name )

        implicit none

        class(eftcamb_rph) :: self
        integer     , intent(in)    :: i
        character(*), intent(out)   :: name

        integer  :: nw, nm, nk, nb, nt
        integer  :: j

        ! compute the incremental number of parameters:
        nw = self%RPH_wDE%parameter_number
        nm = nw + self%RPH_PlanckMass%parameter_number
        nk = nm + self%RPH_Kineticity%parameter_number
        nb = nk + self%RPH_Braiding%parameter_number
        nt = nb + self%RPH_Tensor%parameter_number

        ! check validity of input:
        if ( i > self%parameter_number .or. i <= 0 ) then
            write(*,'(a,I3)') 'No parameter corresponding to: ', i
            write(*,'(a,I3)') 'Total number of parameters is: ', self%parameter_number
            call mpistop('EFTCAMB error')

        ! parameter from wDE function
        else if ( i <= nw ) then
            do j = 1, self%RPH_wDE%parameter_number
                if ( i == j ) call self%RPH_wDE%parameter_names( j, name )
            end do
            return

        ! parameter from RPH_PlanckMass function
        else if ( i <= nm ) then
            do j = 1, self%RPH_PlanckMass%parameter_number
                if ( i-nw == j ) call self%RPH_PlanckMass%parameter_names( j, name )
            end do
            return

        ! parameter from RPH_Kineticity function
        else if ( i <= nk) then
            do j = 1, self%RPH_Kineticity%parameter_number
                if ( i-nm == j ) call self%RPH_Kineticity%parameter_names( j, name )
            end do
            return

        !parameter from RPH_Braiding function
        else if ( i <= nb ) then
            do j = 1, self%RPH_Braiding%parameter_number
                if ( i-nk == j ) call self%RPH_Braiding%parameter_names( j, name )
            end do
            return

        !parameter from RPH_Tensor function
        else if ( i <= nt ) then

            do j = 1, self%RPH_Tensor%parameter_number
                if ( i-nb == j ) call self%RPH_Tensor%parameter_names( j, name )
            end do
            return

        end if

    end subroutine eftcambrphparameternames

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphparameternameslatex( self, i, latexname )

        implicit none

        class(eftcamb_rph) :: self
        integer     , intent(in)    :: i
        character(*), intent(out)   :: latexname

        integer  :: nw, nm, nk, nb, nt
        integer  :: j

        ! compute the incremental number of parameters:
        nw = self%RPH_wDE%parameter_number
        nm = nw + self%RPH_PlanckMass%parameter_number
        nk = nm + self%RPH_Kineticity%parameter_number
        nb = nk + self%RPH_Braiding%parameter_number
        nt = nb + self%RPH_Tensor%parameter_number

        ! check validity of input:
        if ( i > self%parameter_number .or. i <= 0 ) then
            write(*,'(a,I3)') 'No parameter corresponding to: ', i
            write(*,'(a,I3)') 'Total number of parameters is: ', self%parameter_number
            call mpistop('EFTCAMB error')

        ! parameter from wDE function
        else if ( i <= nw ) then
            do j = 1, self%RPH_wDE%parameter_number
                if ( i == j ) call self%RPH_wDE%parameter_names_latex( j, latexname )
            end do
            return

        ! parameter from RPH_PlanckMass function
        else if ( i <= nm ) then
            do j = 1, self%RPH_PlanckMass%parameter_number
                if ( i-nw == j ) call self%RPH_PlanckMass%parameter_names_latex( j, latexname )
            end do
            return

        ! parameter from RPH_Kineticity function
        else if ( i <= nk) then
            do j = 1, self%RPH_Kineticity%parameter_number
                if ( i-nm == j ) call self%RPH_Kineticity%parameter_names_latex( j, latexname )
            end do
            return

        !parameter from RPH_Braiding function
        else if ( i <= nb ) then
            do j = 1, self%RPH_Braiding%parameter_number
                if ( i-nk == j ) call self%RPH_Braiding%parameter_names_latex( j, latexname )
            end do
            return

        !parameter from RPH_Tensor function
        else if ( i <= nt ) then

            do j = 1, self%RPH_Tensor%parameter_number
                if ( i-nb == j ) call self%RPH_Tensor%parameter_names_latex( j, latexname )
            end do
            return

        end if

    end subroutine eftcambrphparameternameslatex

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphparametervalues( self, i, value )

        implicit none

        class(eftcamb_rph) :: self
        integer , intent(in)        :: i
        real(dl), intent(out)       :: value

        integer  :: nw, nm, nk, nb, nt
        integer  :: j

        ! compute the incremental number of parameters:
        nw = self%RPH_wDE%parameter_number
        nm = nw + self%RPH_PlanckMass%parameter_number
        nk = nm + self%RPH_Kineticity%parameter_number
        nb = nk + self%RPH_Braiding%parameter_number
        nt = nb + self%RPH_Tensor%parameter_number

        ! check validity of input:
        if ( i > self%parameter_number .or. i <= 0 ) then
            write(*,'(a,I3)') 'No parameter corresponding to: ', i
            write(*,'(a,I3)') 'Total number of parameters is: ', self%parameter_number
            call mpistop('EFTCAMB error')

        ! parameter from wDE function
        else if ( i <= nw ) then
            do j = 1, self%RPH_wDE%parameter_number
                if ( i == j ) call self%RPH_wDE%parameter_value( j, value )
            end do
            return

        ! parameter from RPH_PlanckMass function
        else if ( i <= nm ) then
            do j = 1, self%RPH_PlanckMass%parameter_number
                if ( i-nw == j ) call self%RPH_PlanckMass%parameter_value( j, value )
            end do
            return

        ! parameter from RPH_Kineticity function
        else if ( i <= nk) then
            do j = 1, self%RPH_Kineticity%parameter_number
                if ( i-nm == j ) call self%RPH_Kineticity%parameter_value( j, value )
            end do
            return

        !parameter from RPH_Braiding function
        else if ( i <= nb ) then
            do j = 1, self%RPH_Braiding%parameter_number
                if ( i-nk == j ) call self%RPH_Braiding%parameter_value( j, value )
            end do
            return

        !parameter from RPH_Tensor function
        else if ( i <= nt ) then

            do j = 1, self%RPH_Tensor%parameter_number
                if ( i-nb == j ) call self%RPH_Tensor%parameter_value( j, value )
            end do
            return

        end if

    end subroutine eftcambrphparametervalues

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphbackgroundeftfunctions( self, a, eft_par_cache, eft_cache )

        implicit none

        class(eftcamb_rph)                           :: self
        real(dl), intent(in)                         :: a
        type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
        type(eftcamb_timestep_cache ), intent(inout) :: eft_cache

        real(dl) :: rph_pm_v, rph_at_v, rph_pm_p, rph_at_p, rph_pm_pp, rph_at_pp, rph_pm_ppp, rph_at_ppp

        ! precompute some functions:
        rph_pm_v               = self%RPH_PlanckMass%value(a)
        rph_at_v               = self%RPH_Tensor%value(a)
        rph_pm_p               = self%RPH_PlanckMass%first_derivative(a)
        rph_at_p               = self%RPH_Tensor%first_derivative(a)
        rph_pm_pp              = self%RPH_PlanckMass%second_derivative(a)
        rph_at_pp              = self%RPH_Tensor%second_derivative(a)
        rph_pm_ppp             = self%RPH_PlanckMass%third_derivative(a)
        rph_at_ppp             = self%RPH_Tensor%third_derivative(a)
        ! compute the EFT functions:
        eft_cache%EFTOmegaV    = rph_pm_v +rph_at_v*(1._dl +rph_pm_v)
        eft_cache%EFTOmegaP    = rph_pm_v*rph_at_p +(1._dl +rph_at_v)*rph_pm_p
        eft_cache%EFTOmegaPP   = rph_pm_v*rph_at_pp +2._dl*rph_pm_p*rph_at_p +(1._dl +rph_at_v)*rph_pm_pp
        eft_cache%EFTOmegaPPP  = rph_pm_v*rph_at_ppp +3._dl*rph_pm_pp*rph_at_p +3._dl*rph_pm_p*rph_at_pp &
            & +(1._dl +rph_at_v)*rph_pm_ppp
        eft_cache%EFTc         = ( eft_cache%adotoa**2 - eft_cache%Hdot )*( eft_cache%EFTOmegaV + 0.5_dl*a*eft_cache%EFTOmegaP ) &
            & -0.5_dl*( a*eft_cache%adotoa )**2*eft_cache%EFTOmegaPP&
            & +0.5_dl*eft_cache%grhov_t*( 1._dl+ self%RPH_wDE%value(a) )
        eft_cache%EFTLambda    = -eft_cache%EFTOmegaV*( 2._dl*eft_cache%Hdot +eft_cache%adotoa**2 ) &
            & -a*eft_cache%EFTOmegaP*( 2._dl*eft_cache%adotoa**2 + eft_cache%Hdot ) &
            & -( a*eft_cache%adotoa )**2*eft_cache%EFTOmegaPP &
            & +self%RPH_wDE%value(a)*eft_cache%grhov_t
        eft_cache%EFTcdot      = +0.5_dl*eft_cache%adotoa*eft_cache%grhov_t*( -3._dl*(1._dl +self%RPH_wDE%value(a))**2 + a*self%RPH_wDE%first_derivative(a) ) &
            & -eft_cache%EFTOmegaV*( eft_cache%Hdotdot -4._dl*eft_cache%adotoa*eft_cache%Hdot +2._dl*eft_cache%adotoa**3 ) &
            & +0.5_dl*a*eft_cache%EFTOmegaP*( -eft_cache%Hdotdot +eft_cache%adotoa*eft_cache%Hdot +eft_cache%adotoa**3) &
            & +0.5_dl*a**2*eft_cache%adotoa*eft_cache%EFTOmegaPP*( eft_cache%adotoa**2 -3._dl*eft_cache%Hdot ) &
            & -0.5_dl*(a*eft_cache%adotoa)**3*eft_cache%EFTOmegaPPP
        eft_cache%EFTLambdadot = -2._dl*eft_cache%EFTOmegaV*( eft_cache%Hdotdot -eft_cache%adotoa*eft_cache%Hdot -eft_cache%adotoa**3 ) &
            & -a*eft_cache%EFTOmegaP*( +eft_cache%Hdotdot +5._dl*eft_cache%adotoa*eft_cache%Hdot -eft_cache%adotoa**3  ) &
            & -a**2*eft_cache%EFTOmegaPP*eft_cache%adotoa*( +2._dl*eft_cache%adotoa**2 +3._dl*eft_cache%Hdot )&
            & -(a*eft_cache%adotoa)**3*eft_cache%EFTOmegaPPP &
            & +eft_cache%grhov_t*eft_cache%adotoa*( a*self%RPH_wDE%first_derivative(a) -3._dl*self%RPH_wDE%value(a)*(1._dl +self%RPH_wDE%value(a) ))

    end subroutine eftcambrphbackgroundeftfunctions

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphsecondordereftfunctions( self, a, eft_par_cache, eft_cache )

        implicit none

        class(eftcamb_rph)                           :: self
        real(dl), intent(in)                         :: a
        type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
        type(eftcamb_timestep_cache ), intent(inout) :: eft_cache

        real(dl) :: rph_pm_v, rph_at_v, rph_ak_v, rph_ab_v, rph_pm_p, rph_at_p, rph_ak_p, rph_ab_p
        real(dl) :: rph_pm_pp, rph_at_pp

        ! precompute some functions:
        rph_pm_v               = self%RPH_PlanckMass%value(a)
        rph_at_v               = self%RPH_Tensor%value(a)
        rph_ak_v               = self%RPH_Kineticity%value(a)
        rph_ab_v               = self%RPH_Braiding%value(a)
        rph_pm_p               = self%RPH_PlanckMass%first_derivative(a)
        rph_at_p               = self%RPH_Tensor%first_derivative(a)
        rph_ak_p               = self%RPH_Kineticity%first_derivative(a)
        rph_ab_p               = self%RPH_Braiding%first_derivative(a)
        rph_pm_pp              = self%RPH_PlanckMass%second_derivative(a)
        rph_at_pp              = self%RPH_Tensor%second_derivative(a)
        ! compute the EFT functions:
        eft_cache%EFTGamma1V  = 0.25_dl*( rph_ak_v*(1._dl +rph_pm_v)*eft_cache%adotoa**2 &
            & -2._dl*eft_cache%EFTc )/(eft_par_cache%h0_Mpc**2*a**2)
        eft_cache%EFTGamma1P  = - 0.5_dl*( rph_ak_v*(1._dl +rph_pm_v)*eft_cache%adotoa**2 &
            & -2._dl*eft_cache%EFTc )/(eft_par_cache%h0_Mpc**2*a**3) &
            & +0.25_dl*( rph_ak_p*(1._dl +rph_pm_v)*eft_cache%adotoa**2 &
            & +rph_ak_v*rph_pm_p*eft_cache%adotoa**2 &
            & +2._dl*rph_ak_v*(1._dl +rph_pm_v)*eft_cache%Hdot/a &
            & -4._dl*eft_cache%EFTc -2._dl*eft_cache%EFTcdot/a/eft_cache%adotoa )/(eft_par_cache%h0_Mpc**2*a**2)
        eft_cache%EFTGamma2V  = ( +2._dl*rph_ab_v*(1._dl +rph_pm_v) &
            & -a*eft_cache%EFTOmegaP )*eft_cache%adotoa/(eft_par_cache%h0_Mpc*a)
        eft_cache%EFTGamma2P  = -0.5_dl*( +2._dl*rph_ab_v*(1._dl +rph_pm_v) &
            & -a*eft_cache%EFTOmegaP )*eft_cache%adotoa/(eft_par_cache%h0_Mpc*a) &
            & -( -2._dl*(1._dl +rph_pm_v)*( rph_ab_p*eft_cache%adotoa**2 &
            & + rph_ab_v*eft_cache%Hdot/a) &
            & - 2._dl*rph_ab_v*eft_cache%adotoa**2*rph_pm_p &
            & + eft_cache%EFTOmegaP*( eft_cache%adotoa**2 +eft_cache%Hdot ) &
            & + a*eft_cache%adotoa**2*eft_cache%EFTOmegaPP )/(eft_par_cache%h0_Mpc*a*eft_cache%adotoa)
        eft_cache%EFTGamma3V  = -rph_at_v*(1._dl +rph_pm_v)
        eft_cache%EFTGamma3P  = -rph_pm_p*rph_at_v -(1._dl +rph_pm_v)*rph_at_p
        eft_cache%EFTGamma4V  = -eft_cache%EFTGamma3V
        eft_cache%EFTGamma4P  = -eft_cache%EFTGamma3P
        eft_cache%EFTGamma4PP = +(1._dl +rph_pm_v)*rph_at_pp +rph_pm_pp*rph_at_v &
            & +2._dl*rph_pm_p*rph_at_p
        eft_cache%EFTGamma5V  = +0.5_dl*eft_cache%EFTGamma3V
        eft_cache%EFTGamma5P  = +0.5_dl*eft_cache%EFTGamma3P
        eft_cache%EFTGamma6V  = 0._dl
        eft_cache%EFTGamma6P  = 0._dl

    end subroutine eftcambrphsecondordereftfunctions

    ! ---------------------------------------------------------------------------------------------
    function eftcambrphcomputedtauda( self, a, eft_par_cache, eft_cache )

        implicit none

        class(eftcamb_rph)                           :: self
        real(dl), intent(in)                         :: a
        type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
        type(eftcamb_timestep_cache ), intent(inout) :: eft_cache

        real(dl) :: eftcambrphcomputedtauda

        real(dl) :: temp

        temp = eft_cache%grhoa2 +eft_par_cache%grhov*a*a*self%RPH_wDE%integral(a)
        eftcambrphcomputedtauda = sqrt(3/temp)

    end function eftcambrphcomputedtauda

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphcomputeadotoa( self, a, eft_par_cache, eft_cache )

        implicit none

        class(eftcamb_rph)                           :: self
        real(dl), intent(in)                         :: a
        type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
        type(eftcamb_timestep_cache ), intent(inout) :: eft_cache

        eft_cache%grhov_t = eft_par_cache%grhov*self%RPH_wDE%integral(a)
        eft_cache%adotoa  = sqrt( ( eft_cache%grhom_t +eft_cache%grhov_t )/3._dl )

    end subroutine eftcambrphcomputeadotoa

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambrphcomputehubbleder( self, a, eft_par_cache, eft_cache )

        implicit none

        class(eftcamb_rph)                           :: self
        real(dl), intent(in)                         :: a
        type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
        type(eftcamb_timestep_cache ), intent(inout) :: eft_cache

        eft_cache%gpiv_t  = self%RPH_wDE%value(a)*eft_cache%grhov_t
        eft_cache%Hdot    = -0.5_dl*( eft_cache%adotoa**2 +eft_cache%gpresm_t +eft_cache%gpiv_t )
        eft_cache%Hdotdot = eft_cache%adotoa*( ( eft_cache%grhob_t +eft_cache%grhoc_t)/6._dl +2._dl*( eft_cache%grhor_t +eft_cache%grhog_t)/3._dl ) &
            & +eft_cache%adotoa*eft_cache%grhov_t*( 1._dl/6._dl +self%RPH_wDE%value(a) +1.5_dl*self%RPH_wDE%value(a)**2 -0.5_dl*a*self%RPH_wDE%first_derivative(a) ) &
            & +eft_cache%adotoa*eft_cache%grhonu_tot/6._dl -0.5_dl*eft_cache%adotoa*eft_cache%gpinu_tot -0.5_dl*eft_cache%gpinudot_tot

    end subroutine eftcambrphcomputehubbleder

    ! ---------------------------------------------------------------------------------------------
    function eftcambrphadditionalmodelstability( self, a, eft_par_cache, eft_cache )

        implicit none

        class(eftcamb_rph)                           :: self
        real(dl), intent(in)                         :: a
        type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
        type(eftcamb_timestep_cache ), intent(inout) :: eft_cache

        logical :: eftcambrphadditionalmodelstability

        eftcambrphadditionalmodelstability = .true.
        if ( self%RPH_wDE%value(a) > -1._dl/3._dl ) eftcambrphadditionalmodelstability = .false.

    end function eftcambrphadditionalmodelstability

    ! ---------------------------------------------------------------------------------------------

end module eftcamb_reparametrized_horndeski