09p1_Designer_fR.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_designer_fr

    use precision
    use inifile
    use amlutils
    use equispaced_linear_interpolation_1d
    use eft_def
    use eftcamb_rootfind
    use eftcamb_cache
    use eftcamb_abstract_parametrizations_1d
    use eftcamb_neutral_parametrization_1d
    use eftcamb_constant_parametrization_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_fr_designer

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

        ! theory parameters:
        real(dl) :: b0

        ! the pure EFT functions model selection flags:
        integer  :: eftwde

        ! the pure EFT functions:
        class( parametrized_function_1d ), allocatable :: desfrwde

        ! the interpolated EFT functions that come out of the background sover:
        type(equispaced_linear_interpolate_function_1d) :: eftomega
        type(equispaced_linear_interpolate_function_1d) :: eftlambda

        ! some designer parameters:
        integer  :: designer_num_points = 1000
        real(dl) :: x_initial           = log(10._dl**(-8._dl))           
        real(dl) :: x_final             = 0.0_dl                          

    contains

        ! initialization of the model:
        procedure :: read_model_selection            => eftcambdesignerfrreadmodelselectionfromfile   
        procedure :: allocate_model_selection        => eftcambdesignerfrallocatemodelselection       
        procedure :: init_model_parameters           => eftcambdesignerfrinitmodelparameters          
        procedure :: init_model_parameters_from_file => eftcambdesignerfrinitmodelparametersfromfile  

        ! background solver:
        procedure :: initialize_background           => eftcambdesignerfrinitbackground               
        procedure :: solve_designer_equations        => eftcambdesignerfrsolvedesignerequations       
        procedure :: find_initial_conditions         => eftcambdesignerfrfindinitialconditions        

        ! utility functions:
        procedure :: compute_param_number  => eftcambdesignerfrcomputeparametersnumber     
        procedure :: feedback              => eftcambdesignerfrfeedback                    
        procedure :: parameter_names       => eftcambdesignerfrparameternames              
        procedure :: parameter_names_latex => eftcambdesignerfrparameternameslatex         
        procedure :: parameter_values      => eftcambdesignerfrparametervalues             

        ! CAMB related procedures:
        procedure :: compute_background_eft_functions  => eftcambdesignerfrbackgroundeftfunctions   
        procedure :: compute_secondorder_eft_functions => eftcambdesignerfrsecondordereftfunctions  
        procedure :: compute_dtauda                    => eftcambdesignerfrcomputedtauda            
        procedure :: compute_adotoa                    => eftcambdesignerfrcomputeadotoa            
        procedure :: compute_h_derivs                  => eftcambdesignerfrcomputehubbleder         

        ! stability procedures:
        procedure :: additional_model_stability        => eftcambdesignerfradditionalmodelstability 

    end type eftcamb_fr_designer

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

contains

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

        implicit none

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

        ! read model selection flags:
        self%EFTwDE             = ini_read_int_file( ini, 'EFTwDE', 0 )

    end subroutine eftcambdesignerfrreadmodelselectionfromfile

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambdesignerfrallocatemodelselection( self )

        implicit none

        class(eftcamb_fr_designer)                        :: self
        character, allocatable, dimension(:)              :: param_names
        character, allocatable, dimension(:)              :: param_names_latex

        ! allocate wDE:
        if ( allocated(self%DesfRwDE) ) deallocate(self%DesfRwDE)
        select case ( self%EFTwDE )
            case(0)
                allocate( wde_lcdm_parametrization_1d::self%DesfRwDE )
            case(1)
                allocate( constant_parametrization_1d::self%DesfRwDE )
            case(2)
                allocate( cpl_parametrization_1d::self%DesfRwDE )
                call self%DesfRwDE%set_param_names( ['EFTw0', 'EFTwa'], ['w_0', 'w_a'] )
            case(3)
                allocate( jbp_parametrization_1d::self%DesfRwDE )
                call self%DesfRwDE%set_param_names( ['EFTw0', 'EFTwa', 'EFTwn'], [ 'w_0', 'w_a', 'n  ' ] )
            case(4)
                allocate( turning_point_parametrization_1d::self%DesfRwDE )
                call self%DesfRwDE%set_param_names( ['EFTw0 ', 'EFTwa ', 'EFTwat'], ['w_0', 'w_a', 'a_t'] )
            case(5)
                allocate( taylor_parametrization_1d::self%DesfRwDE )
                call self%DesfRwDE%set_param_names( ['EFTw0', 'EFTwa', 'EFTw2', 'EFTw3'], ['w_0', 'w_a', 'w_2', 'w_3'] )
            case default
                write(*,'(a,I3)') 'No model corresponding to EFTwDE =', self%EFTwDE
                write(*,'(a)')    'Please select an appropriate model.'
        end select

        ! initialize the names:
        call self%DesfRwDE%set_name( 'EFTw', 'w' )

    end subroutine eftcambdesignerfrallocatemodelselection

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

        implicit none

        class(eftcamb_fr_designer)                             :: self
        real(dl), dimension(self%parameter_number), intent(in) :: array
        real(dl), dimension(self%parameter_number -1)          :: temp
        integer                                                :: i

        self%B0 = array(1)

        do i = 1, self%parameter_number -1
            temp(i) = array(i+1)
        end do
        call self%DesfRwDE%init_parameters(temp)

    end subroutine eftcambdesignerfrinitmodelparameters

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

        implicit none

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

        ! read B0:
        self%B0 = ini_read_double_file( ini, 'EFTB0', 0._dl )
        ! read w_DE parameters:
        call self%DesfRwDE%init_from_file( ini )

    end subroutine eftcambdesignerfrinitmodelparametersfromfile

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambdesignerfrinitbackground( self, params_cache, feedback_level, success )

        implicit none

        class(eftcamb_fr_designer)                   :: self
        type(eftcamb_parameter_cache), intent(in)    :: params_cache
        integer                      , intent(in)    :: feedback_level
        logical                      , intent(out)   :: success

        real(dl) :: a_ini, b0
        real(dl) :: tempmin, tempmax, debug_a
        integer  :: debug_maxnum, debug_n

        ! some feedback:
        if ( feedback_level>0 ) then
            write(*,'(a)') "***************************************************************"
            write(*,'(a)') ' EFTCAMB designer f(R) background solver'
            write(*,'(a)')
        end if

        ! initialize interpolating functions:
        call self%EFTOmega%initialize  ( self%designer_num_points, self%x_initial, self%x_final )
        call self%EFTLambda%initialize ( self%designer_num_points, self%x_initial, self%x_final )

        ! debug code:
        if ( debugeftcamb ) then
            ! print the function B0(A). This is used to debug the initial conditions part.
            call createtxtfile( './debug_designer_fR_B.dat', 34 )
            print*, 'EFTCAMB DEBUG ( f(R) designer ): Printing B(A) results'
            tempmin      = -10._dl
            tempmax      = +10._dl
            debug_maxnum = 1000
            do debug_n = 1, debug_maxnum
                debug_a = tempmin +REAL(debug_n-1)*(tempmax-tempmin)/REAL(debug_maxnum-1)
                call self%solve_designer_equations( params_cache, debug_a, b0, only_b0=.true., success=success )
                write(34,*) debug_a, b0
            end do
            close(34)
            ! prints f(R) quantities.
            print*, 'EFTCAMB DEBUG ( f(R) designer ): Printing F(R) results'
            call createtxtfile( './debug_designer_fR_solution.dat', 33 )
            debug_a = 1.0_dl
            call self%solve_designer_equations( params_cache, debug_a, b0, only_b0=.false., success=success )
            close(33)
        end if

        ! call boundary conditions lookup:
        call self%find_initial_conditions( params_cache, feedback_level, a_ini, success )

        ! solve the background equations and store the solution:
        call self%solve_designer_equations( params_cache, a_ini, b0, only_b0=.false., success=success )

    end subroutine eftcambdesignerfrinitbackground

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambdesignerfrsolvedesignerequations( self, params_cache, A, B0, only_B0, success )

        implicit none

        class(eftcamb_fr_designer)                   :: self
        type(eftcamb_parameter_cache), intent(in)    :: params_cache
        real(dl), intent(in)                         :: a
        real(dl), intent(out)                        :: b0
        logical , optional                           :: only_b0
        logical , intent(out)                        :: success

        integer, parameter :: num_eq = 2

        real(dl) :: omegam_eft, omegavac_eft, omegamassivenu_eft, omegagamma_eft, omeganu_eft
        real(dl) :: omegarad_eft, equivalencescale_fr, ratio_fr, initial_b_fr, initial_c_fr
        real(dl) :: pplus, yplus, coeffa_part, ystar, x

        real(dl) :: y(num_eq), ydot(num_eq)

        integer  :: itol, itask, istate, iopt, lrn, lrs, lrw, lis, lin, liw, jacobianmode, i
        real(dl) :: rtol, atol, t1, t2, b
        real(dl), allocatable :: rwork(:)
        integer,  allocatable :: iwork(:)

        ! digedt the input:
        if ( .not. present(only_b0) ) only_b0 = .false.

        ! 1) Cosmological parameters:
        omegam_eft         = params_cache%omegab + params_cache%omegac
        omegavac_eft       = params_cache%omegav
        omegamassivenu_eft = params_cache%omegan
        omegagamma_eft     = params_cache%omegag
        omeganu_eft        = params_cache%omegar

        omegarad_eft       = omegagamma_eft + omeganu_eft

        equivalencescale_fr = (omegarad_eft+ omegamassivenu_eft)/omegam_eft
        ratio_fr = equivalencescale_fr/exp( self%x_initial )

        ! 2) Growing mode solution:
        initial_b_fr = (7._dl + 8._dl*ratio_fr)/(2._dl*(1._dl + ratio_fr))
        initial_c_fr = -3._dl/(2._dl*(1._dl + ratio_fr))
        pplus = 0.5_dl*(-initial_b_fr + sqrt(initial_b_fr**2 - 4._dl*initial_c_fr))
        yplus = exp(pplus*self%x_initial)
        !    Construction of the particolar solution:
        coeffa_part = (-6._dl*initial_c_fr)/(-3._dl*exp(self%x_initial)*self%DesfRwDE%first_derivative( exp(self%x_initial) )&
            & +9._dl*self%DesfRwDE%value( exp(self%x_initial) )**2&
            & +(18._dl-3._dl*initial_b_fr)*self%DesfRwDE%value( exp(self%x_initial) )&
            & +9._dl -3._dl*initial_b_fr +initial_c_fr)
        ystar = coeffa_part*omegavac_eft*exp(-2._dl*self%x_initial)*self%DesfRwDE%integral( exp(self%x_initial) )

        ! 3) Set initial conditions:
        x    = self%x_initial
        y(1) = a*yplus + ystar
        y(2) = pplus*a*yplus - 3._dl*( 1._dl+ self%DesfRwDE%value(exp(self%x_initial)) )*ystar
        ydot = 0._dl

        ! 4) Initialize DLSODA:
        ! set-up the relative and absolute tollerances:
        itol = 1
        rtol = 1.d-10
        atol = 1.d-14
        ! initialize task to do:
        itask  = 1
        istate = 1
        iopt   = 1
        ! initialize the work space:
        lrn = 20 + 16*num_eq
        lrs = 22 + 9*num_eq + num_eq**2
        lrw = max(lrn,lrs)
        lis = 20 + num_eq
        lin = 20
        liw = max(lis,lin)
        ! allocate the arrays:
        allocate(rwork(lrw))
        allocate(iwork(liw))
        ! optional lsoda input:
        rwork(5) = 0._dl  ! the step size to be attempted on the first step. The default value is determined by the solver.
        rwork(6) = 0._dl  ! the maximum absolute step size allowed. The default value is infinite.
        rwork(7) = 0._dl  ! the minimum absolute step size allowed. The default value is 0.
        iwork(5) = 0      ! flag to generate extra printing at method switches. IXPR = 0 means no extra printing (the default). IXPR = 1 means print data on each switch.
        iwork(6) = 100    ! maximum number of (internally defined) steps allowed during one call to the solver. The default value is 500.
        iwork(7) = 0      ! maximum number of messages printed (per problem) warning that T + H = T on a step (H = step size). This must be positive to result in a non-default value.  The default value is 10.
        iwork(8) = 0      ! the maximum order to be allowed for the nonstiff (Adams) method.  the default value is 12.
        iwork(9) = 0      ! the maximum order to be allowed for the stiff (BDF) method.  The default value is 5.
        ! additional lsoda stuff:
        CALL xsetf(0) ! suppress odepack printing
        ! Jacobian mode: 1=fullJacobian, 2=not provided
        jacobianmode = 2

        ! store the first value of the EFT functions:
        t1  = self%EFTOmega%x(1)
        ! compute output EFT functions if needed:
        if ( .not. only_b0 ) then
            call output( num_eq, 1, t1, y, b0 )
        end if

        ! 5) solve the equations:
        do i=1, self%EFTOmega%num_points-1

            ! set the time step:
            t1 = self%EFTOmega%x(i)
            t2 = self%EFTOmega%x(i+1)
            ! solve the system:
            call dlsoda ( derivs, num_eq, y, t1, t2, itol, rtol, atol, itask, istate, iopt, rwork, lrw, iwork, liw, jacobian, jacobianmode)
            ! check istate for LSODA good completion:
            if ( istate < 0 ) then
                success = .false.
                return
            end if
            ! compute output EFT functions if needed:
            if ( .not. only_b0 ) then
                call output( num_eq, i+1, t2, y, b0 )
            end if

        end do

        ! compute B0:
        if ( only_b0 ) then
            call output( num_eq, self%EFTOmega%num_points, t2, y, b0 )
        end if

        return

    contains

        ! ---------------------------------------------------------------------------------------------
        subroutine derivs( num_eq, x, y, ydot )

            implicit none

            integer , intent(in)                     :: num_eq
            real(dl), intent(in)                     :: x
            real(dl), intent(in) , dimension(num_eq) :: y
            real(dl), intent(out), dimension(num_eq) :: ydot

            real(dl) :: a, eft_e_gfun, eft_e_gfunp, eft_e_gfunpp, eft_e_gfunppp
            real(dl) :: rhonu_tot, presnu_tot, presnudot_tot, presnudotdot_tot
            real(dl) :: eft_e_nu, eft_ep_nu, eft_epp_nu, eft_e3p_nu, rhonu, presnu, grhormass_t
            real(dl) :: efunction, efunprime, adotoa, hdot, presnudot, presnudotdot
            real(dl) :: efunprime2, efunprime3
            integer  :: nu_i

            ! 1) convert x in a:
            a = exp(x)

            ! 2) Compute the function g(x) and its derivatives:
            eft_e_gfun    = -(log( self%DesfRwDE%integral(a) ) -2._dl*x)/3._dl
            eft_e_gfunp   = 1._dl +self%DesfRwDE%value(a)
            eft_e_gfunpp  = exp(x)*self%DesfRwDE%first_derivative(a)
            eft_e_gfunppp = exp(x)*self%DesfRwDE%first_derivative(a) +exp(2._dl*x)*self%DesfRwDE%second_derivative(a)

            ! 3) Compute energy and its derivatives:
            ! First compute massive neutrinos contribution:
            rhonu_tot  = 0._dl
            presnu_tot = 0._dl
            eft_e_nu   = 0._dl
            eft_ep_nu  = 0._dl
            if ( params_cache%Num_Nu_Massive /= 0) then
                do nu_i = 1, params_cache%Nu_mass_eigenstates

                    rhonu  = 0._dl
                    presnu = 0._dl
                    grhormass_t= params_cache%grhormass(nu_i)/a**2
                    call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
                    rhonu_tot  = rhonu_tot + grhormass_t*rhonu
                    presnu_tot = presnu_tot + grhormass_t*presnu

                    eft_e_nu   = eft_e_nu  + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
                    eft_ep_nu  = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)

                end do
            end if

            ! Add its contribution to E and E':
            efunction = +omegarad_eft*exp(-4._dl*x)&
                & +omegam_eft*exp(-3._dl*x)&
                & +omegavac_eft*exp(-3._dl*eft_e_gfun) + eft_e_nu
            efunprime = -4._dl*omegarad_eft*exp(-4._dl*x)&
                & -3._dl*omegam_eft*exp(-3._dl*x)&
                & -3._dl*omegavac_eft*eft_e_gfunp*exp(-3._dl*eft_e_gfun) +eft_ep_nu

            ! Compute everything of massive nu again to get the time derivatives:
            rhonu_tot        = 0._dl
            presnu_tot       = 0._dl
            presnudot_tot    = 0._dl
            presnudotdot_tot = 0._dl
            eft_e_nu   = 0._dl
            eft_ep_nu  = 0._dl
            eft_epp_nu = 0._dl
            eft_e3p_nu = 0._dl
            if ( params_cache%Num_Nu_Massive /= 0 ) then
                do nu_i = 1, params_cache%Nu_mass_eigenstates

                    adotoa = +a*params_cache%h0_Mpc*sqrt(efunction)
                    hdot   = +0.5_dl*params_cache%h0_Mpc**2*a**2*efunprime +adotoa**2

                    rhonu        = 0._dl
                    presnu       = 0._dl
                    presnudot    = 0._dl
                    presnudotdot = 0._dl

                    grhormass_t = params_cache%grhormass(nu_i)/a**2

                    call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
                    presnudot = params_cache%Nu_pidot(a*params_cache%nu_masses(nu_i),adotoa,presnu)
                    presnudotdot = params_cache%Nu_pidotdot(a*params_cache%nu_masses(nu_i),adotoa,hdot,presnu,presnudot)

                    rhonu_tot  = rhonu_tot + grhormass_t*rhonu
                    presnu_tot = presnu_tot + grhormass_t*presnu
                    presnudot_tot  = presnudot_tot + grhormass_t*(presnudot -4._dl*adotoa*presnu)
                    presnudotdot_tot = presnudotdot_tot + grhormass_t*(presnudotdot &
                        & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))

                    eft_e_nu   = eft_e_nu + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
                    eft_ep_nu  = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)
                    eft_epp_nu = eft_epp_nu + 3._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
                        & -grhormass_t*(presnudot -4._dl*adotoa*presnu)/params_cache%h0_Mpc**3/sqrt(efunction)/a**3
                    eft_e3p_nu = eft_e3p_nu -9._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
                        & +(3._dl/adotoa/params_cache%h0_Mpc**2/a**2+hdot/adotoa**3/params_cache%h0_Mpc**2/a**2)&
                        &*grhormass_t*(presnudot -4._dl*adotoa*presnu)&
                        & -grhormass_t*(presnudotdot &
                        & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))/adotoa**2/params_cache%h0_Mpc**2/a**2
                end do
            end if

            efunprime2 = 16._dl*omegarad_eft*exp(-4._dl*x)&
                & +9._dl*omegam_eft*exp(-3._dl*x)&
                & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(eft_e_gfunpp -3._dl*eft_e_gfunp**2) + eft_epp_nu
            efunprime3 = -64._dl*omegarad_eft*exp(-4._dl*x)&
                & -27._dl*omegam_eft*exp(-3._dl*x)&
                & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*&
                &(eft_e_gfunppp-9._dl*eft_e_gfunp*eft_e_gfunpp+9._dl*eft_e_gfunp**3) + eft_e3p_nu

            ! 4) Get the equation of motion:
            ydot(1) = y(2)
            ydot(2) = (1._dl+0.5_dl*efunprime/efunction+(4._dl*efunprime2+efunprime3)/(4._dl*efunprime+efunprime2))*y(2) &
                & -0.5_dl*(4._dl*efunprime+efunprime2)/efunction*y(1) &
                & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(4._dl*efunprime+efunprime2)/efunction

        end subroutine

        ! ---------------------------------------------------------------------------------------------
        subroutine jacobian( num_eq, x, y, ml, mu, pd, nrowpd )

            implicit none

            integer                            :: num_eq
            integer                            :: ml
            integer                            :: mu
            integer                            :: nrowpd
            real(dl)                           :: x
            real(dl), dimension(num_eq)        :: y
            real(dl), dimension(nrowpd,num_eq) :: pd

        end subroutine jacobian

        ! ---------------------------------------------------------------------------------------------
        subroutine output( num_eq, ind, x, y, B )

            implicit none

            integer , intent(in)                     :: num_eq
            integer , intent(in)                     :: ind
            real(dl), intent(in)                     :: x
            real(dl), intent(in) , dimension(num_eq) :: y
            real(dl), intent(out)                    :: b

            real(dl) :: ydot(num_eq)
            real(dl) :: a, eft_e_gfun, eft_e_gfunp, eft_e_gfunpp, eft_e_gfunppp
            real(dl) :: rhonu_tot, presnu_tot, presnudot_tot, presnudotdot_tot
            real(dl) :: eft_e_nu, eft_ep_nu, eft_epp_nu, eft_e3p_nu, rhonu, presnu, grhormass_t
            real(dl) :: efunction, efunprime, adotoa, hdot, presnudot, presnudotdot
            real(dl) :: efunprime2, efunprime3, ricci, f_sub_r, ede, edep, edepp
            integer  :: nu_i
            logical  :: is_open

            ! 1) convert x in a:
            a = exp(x)

            ! 2) Compute the function g(x) and its derivatives:
            eft_e_gfun    = -(log( self%DesfRwDE%integral(a) ) -2._dl*x)/3._dl
            eft_e_gfunp   = 1._dl +self%DesfRwDE%value(a)
            eft_e_gfunpp  = exp(x)*self%DesfRwDE%first_derivative(a)
            eft_e_gfunppp = exp(x)*self%DesfRwDE%first_derivative(a) +exp(2._dl*x)*self%DesfRwDE%second_derivative(a)

            ! 3) Compute energy and its derivatives:
            ! First compute massive neutrinos contribution:
            rhonu_tot  = 0._dl
            presnu_tot = 0._dl
            eft_e_nu   = 0._dl
            eft_ep_nu  = 0._dl
            if ( params_cache%Num_Nu_Massive /= 0) then
                do nu_i = 1, params_cache%Nu_mass_eigenstates

                    rhonu  = 0._dl
                    presnu = 0._dl
                    grhormass_t= params_cache%grhormass(nu_i)/a**2
                    call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
                    rhonu_tot  = rhonu_tot + grhormass_t*rhonu
                    presnu_tot = presnu_tot + grhormass_t*presnu

                    eft_e_nu   = eft_e_nu  + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
                    eft_ep_nu  = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)

                end do
            end if

            ! Add its contribution to E and E':
            efunction = +omegarad_eft*exp(-4._dl*x)&
                & +omegam_eft*exp(-3._dl*x)&
                & +omegavac_eft*exp(-3._dl*eft_e_gfun) + eft_e_nu
            efunprime = -4._dl*omegarad_eft*exp(-4._dl*x)&
                & -3._dl*omegam_eft*exp(-3._dl*x)&
                & -3._dl*omegavac_eft*eft_e_gfunp*exp(-3._dl*eft_e_gfun) +eft_ep_nu

            ! Compute everything of massive nu again to get the time derivatives:
            rhonu_tot        = 0._dl
            presnu_tot       = 0._dl
            presnudot_tot    = 0._dl
            presnudotdot_tot = 0._dl
            eft_e_nu   = 0._dl
            eft_ep_nu  = 0._dl
            eft_epp_nu = 0._dl
            eft_e3p_nu = 0._dl
            if ( params_cache%Num_Nu_Massive /= 0 ) then
                do nu_i = 1, params_cache%Nu_mass_eigenstates

                    adotoa = +a*params_cache%h0_Mpc*sqrt(efunction)
                    hdot   = +0.5_dl*params_cache%h0_Mpc**2*a**2*efunprime +adotoa**2

                    rhonu        = 0._dl
                    presnu       = 0._dl
                    presnudot    = 0._dl
                    presnudotdot = 0._dl

                    grhormass_t  = params_cache%grhormass(nu_i)/a**2

                    call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
                    presnudot = params_cache%Nu_pidot(a*params_cache%nu_masses(nu_i),adotoa,presnu)
                    presnudotdot = params_cache%Nu_pidotdot(a*params_cache%nu_masses(nu_i),adotoa,hdot,presnu,presnudot)

                    rhonu_tot  = rhonu_tot + grhormass_t*rhonu
                    presnu_tot = presnu_tot + grhormass_t*presnu
                    presnudot_tot  = presnudot_tot + grhormass_t*(presnudot -4._dl*adotoa*presnu)
                    presnudotdot_tot = presnudotdot_tot + grhormass_t*(presnudotdot &
                        & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))

                    eft_e_nu   = eft_e_nu + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
                    eft_ep_nu  = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)
                    eft_epp_nu = eft_epp_nu + 3._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
                        & -grhormass_t*(presnudot -4._dl*adotoa*presnu)/params_cache%h0_Mpc**3/sqrt(efunction)/a**3
                    eft_e3p_nu = eft_e3p_nu -9._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
                        & +(3._dl/adotoa/params_cache%h0_Mpc**2/a**2+hdot/adotoa**3/params_cache%h0_Mpc**2/a**2)&
                        &*grhormass_t*(presnudot -4._dl*adotoa*presnu)&
                        & -grhormass_t*(presnudotdot &
                        & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))/adotoa**2/params_cache%h0_Mpc**2/a**2
                end do
            end if

            efunprime2 = 16._dl*omegarad_eft*exp(-4._dl*x)&
                & +9._dl*omegam_eft*exp(-3._dl*x)&
                & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(eft_e_gfunpp -3._dl*eft_e_gfunp**2) + eft_epp_nu
            efunprime3 = -64._dl*omegarad_eft*exp(-4._dl*x)&
                & -27._dl*omegam_eft*exp(-3._dl*x)&
                & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*&
                &(eft_e_gfunppp-9._dl*eft_e_gfunp*eft_e_gfunpp+9._dl*eft_e_gfunp**3) + eft_e3p_nu

            ! compute E dark energy:
            ede   = +omegavac_eft*exp(-3._dl*eft_e_gfun)
            edep  = -3._dl*omegavac_eft*eft_e_gfunp*exp(-3._dl*eft_e_gfun)
            edepp = -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(eft_e_gfunpp -3._dl*eft_e_gfunp**2)

            ! call derivs to compute ydot at a given time:
            call derivs( num_eq, x, y, ydot )

            ! Compute the Ricci:
            ricci    = 3._dl*(4._dl*efunction+efunprime)
            f_sub_r  = ydot(1)/3._dl/(4._dl*efunprime +efunprime2)

            ! compute B:
            b       = 2._dl/3._dl/(1._dl +f_sub_r)*1._dl/(4._dl*efunprime +efunprime2)*efunction/efunprime*&
                &( ydot(2) -ydot(1)*(4._dl*efunprime2+efunprime3)/(4._dl*efunprime +efunprime2) )

            ! compute the EFT functions:
            self%EFTOmega%y(ind)    = f_sub_r
            self%EFTOmega%yp(ind)   = exp(-x)/(6._dl*efunction*(efunprime2+4._dl*efunprime))*(-efunprime2*(6._dl*ede+y(1))&
                &+efunprime*(ydot(1)-4._dl*(6._dl*ede +y(1)))+2._dl*efunction*ydot(1))
            self%EFTOmega%ypp(ind)  = exp(-2._dl*x)/(12._dl*efunction**2*(efunprime2+4._dl*efunprime))*&
                &(efunprime*(efunprime*(24._dl*ede-ydot(1) +4._dl*y(1))-6._dl*efunction*(8._dl*edep +ydot(1)))&
                &+efunprime2*(efunprime*(6._dl*ede +y(1))-12._dl*efunction*edep))
            self%EFTOmega%yppp(ind) = exp(-3._dl*x)/(24._dl*efunction**3*(efunprime2+4._dl*efunprime))*&
                &(2._dl*efunction*(efunprime2**2*(6._dl*ede+y(1))&
                & +4._dl*efunprime**2*(18._dl*edep+6._dl*ede+2._dl*ydot(1) +y(1))&
                & +efunprime*efunprime2*(18._dl*edep+30._dl*ede -ydot(1)+5._dl*y(1)))&
                & +12._dl*efunction**2*(efunprime2*(-2._dl*edepp+4._dl*edep-ydot(1))&
                & +efunprime*(-8._dl*edepp+16._dl*edep +ydot(1)))&
                & +3._dl*efunprime**2*(efunprime*(ydot(1)-4._dl*(6._dl*ede +y(1)))-efunprime2*(6._dl*ede +y(1))))
            self%EFTLambda%y(ind)   = 0.5_dl*params_cache%h0_Mpc**2*( y(1) -f_sub_r*ricci )*exp(x)**2
            self%EFTLambda%yp(ind)  = -1.5_dl*params_cache%h0_Mpc**3*sqrt(efunction)*(exp(x)**4*self%EFTOmega%yp(ind)*(4._dl*efunction+efunprime))

            if ( debugeftcamb ) then
                inquire( unit=33, opened=is_open )
                if ( is_open ) then
                    write (33,'(20E15.5)') x, exp(x), ricci, y(1), &
                        & self%EFTOmega%y(ind), self%EFTOmega%yp(ind), self%EFTOmega%ypp(ind), self%EFTOmega%yppp(ind), self%EFTLambda%y(ind), self%EFTLambda%yp(ind),&
                        & b
                end if
            end if

        end subroutine

    end subroutine eftcambdesignerfrsolvedesignerequations

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambdesignerfrfindinitialconditions( self, params_cache, feedback_level, A_ini, success )

        implicit none

        class(eftcamb_fr_designer)                   :: self
        type(eftcamb_parameter_cache), intent(in)    :: params_cache
        integer                      , intent(in)    :: feedback_level
        real(dl)                     , intent(out)   :: a_ini
        logical                      , intent(out)   :: success

        real(dl) :: atemp1, atemp2, btemp1, btemp2, horizasyntb
        real(dl) :: vertasynta, atemp3, atemp4, btemp3, btemp4, realap
        integer  :: ind

        !  Find initial conditions for designer F(R).
        !  Initial conditions are found solving B0(A) = B0wanted.
        !  The next part of code is complicated by the fact that B0(A) is not continuous and we have to adopt ad-hoc strategies.

        ! 1) Find the horizontal asymptote of B0(A)
        atemp1 = 0._dl
        atemp2 = 0._dl
        do ind=10, 100, 1
            atemp1 = atemp2
            atemp2 = 10._dl**REAL(ind)
            btemp1 = desfr_bfunca(atemp1)
            btemp2 = desfr_bfunca(atemp2)
            if (abs(btemp1-btemp2)/abs(atemp1-atemp2).lt.1.d-100) exit
        end do
        horizasyntb = btemp2

        if ( feedback_level>1 ) write(*,'(a,E13.4)') '   horizontal asymptote = ', horizasyntb

        ! 2) Check that the value of B0 given is not the forbidden one: the one corresponding to the horizontal asynt.
        if ( abs(horizasyntb-self%B0)<1.d-15 ) then
            success=.false.
            return
        end if

        ! 3) Bracket the vertical asyntote
        atemp1 = -10._dl
        atemp2 = 10._dl
        call zbrac(desfr_bfunca,atemp1,atemp2,success,horizasyntb)
        if (.not.success) then
            if ( feedback_level>1 ) write(*,'(a)') '   FAILURE of vertical asymptote bracketing'
            return
        end if

        ! 4) Find the vertical asyntote by tricking the root finding algorithm.
        vertasynta = zbrent(desfr_bfunca,atemp1,atemp2,1.d-100,horizasyntb,success)
        if (.not.success) then
            if ( feedback_level>1 ) write(*,'(a)') '   FAILURE of vertical asyntote finding'
            return
        end if

        if ( feedback_level>1 ) write(*,'(a,E13.4)') '   vertical asymptote   = ', vertasynta

        ! 5) Find values for A that are on the left and on the right of the asyntote.
        do ind=-10, -1, 1
            atemp1 = vertasynta+10._dl**ind
            atemp2 = vertasynta-10._dl**ind
            btemp1 = desfr_bfunca(atemp1)
            btemp2 = desfr_bfunca(atemp2)
            if (btemp1*btemp2<0._dl) exit
        end do

        ! 6) Extablish on which side of the asyntote to look for the solution.
        atemp3 = atemp1
        atemp4 = atemp2
        do ind=1, 10, 1
            atemp3 = atemp3 + 10._dl**ind
            atemp4 = atemp4 - 10._dl**ind
            btemp3 = desfr_bfunca(atemp3)
            btemp4 = desfr_bfunca(atemp4)
            if ((btemp1-self%B0)*(btemp3-self%B0)<0.or.(btemp2-self%B0)*(btemp4-self%B0)<0) exit
        end do
        if ((btemp1-self%B0)*(btemp3-self%B0)<0.and.(btemp2-self%B0)*(btemp4-self%B0)<0) then
            if ( feedback_level>1 ) write(*,'(a)') '   FAILURE as the root seems to be on both sides'
            return
        end if

        ! 7) Solve the equation B0(A)=B0_wanted.
        if ((btemp1-self%B0)*(btemp3-self%B0)<0) then
            realap = zbrent(desfr_bfunca,atemp1,atemp3,1.d-50,self%B0,success)
            if (.not.success) then
                if ( feedback_level>1 ) write(*,'(a)') '   FAILURE right side solution not found'
                return
            end if
        else if ((btemp2-self%B0)*(btemp4-self%B0)<0) then
            realap = zbrent(desfr_bfunca,atemp2,atemp4,1.d-50,self%B0,success)
            if (.not.success) then
                if ( feedback_level>1 ) write(*,'(a)') '   FAILURE left side solution not found'
                return
            end if
        else
            if ( feedback_level>1 ) write(*,'(a)') '   FAILURE the root was not on the right side nor the left one...'
            return
        end if

        if ( feedback_level>1 ) write(*,'(a,E13.4)') '   initial condition A  = ', realap

        ! 8) Check if the result found is compatible with the requested one. This is required only for debug.
        btemp1 = desfr_bfunca(realap)
        if ( feedback_level>1 ) then
            write(*,'(a,E13.4)') '   B0 found = ', btemp1
            write(*,'(a,E13.4)') '   B0 given = ', self%B0
        end if

        if (abs(btemp1-self%B0)/abs(self%B0)>0.1_dl.and.abs(btemp1-self%B0)>1.d-8.and..false.) then
            if ( feedback_level>1 ) write(*,'(a)')  '   FAILURE designer code unable to find appropriate initial conditions'
            success = .false.
            return
        end if

        ! 9) output the value:
        a_ini = realap

        return

    contains

        ! ---------------------------------------------------------------------------------------------
        function desfr_bfunca(A)

            implicit none

            real(dl) :: a
            real(dl) :: desfr_bfunca
            real(dl) :: b0

            call self%solve_designer_equations( params_cache, a, b0, only_b0=.true., success=success )

            desfr_bfunca = b0

        end function

    end subroutine eftcambdesignerfrfindinitialconditions

    ! ---------------------------------------------------------------------------------------------
    subroutine eftcambdesignerfrcomputeparametersnumber( self )

        implicit none

        class(eftcamb_fr_designer)  :: self

        self%parameter_number = 1
        self%parameter_number = self%parameter_number +self%DesfRwDE%parameter_number

    end subroutine eftcambdesignerfrcomputeparametersnumber

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

        implicit none

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

        write(*,*)
        write(*,'(a,a)')    '   Model               =  ', self%name
        write(*,'(a,I3)')   '   Number of params    ='  , self%parameter_number
        ! print model functions informations:
        if ( self%EFTwDE /= 0 ) then
            write(*,*)
            write(*,'(a,I3)')  '   EFTwDE              =', self%EFTwDE
        end if
        write(*,*)
        write(*,'(a24,F12.6)') '   B0                  =', self%B0

        call self%DesfRwDE%feedback( print_params )

    end subroutine eftcambdesignerfrfeedback

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

        implicit none

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

        ! check validity of input:
        if ( i<=0 .or. i>self%parameter_number ) then
            write(*,'(a,I3)') 'EFTCAMB error: no parameter corresponding to number ', i
            write(*,'(a,I3)') 'Total number of parameters is ', self%parameter_number
            call mpistop('EFTCAMB error')
        ! the first parameter is B0:
        else if ( i==1 ) then
            name = trim('B0')
            return
        ! the other parameters are the w_DE parameters:
        else
            call self%DesfRwDE%parameter_names( i-1, name )
            return
        end if

    end subroutine eftcambdesignerfrparameternames

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

        implicit none

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

        ! check validity of input:
        if ( i<=0 .or. i>self%parameter_number ) then
            write(*,'(a,I3)') 'EFTCAMB error: no parameter corresponding to number ', i
            write(*,'(a,I3)') 'Total number of parameters is ', self%parameter_number
            call mpistop('EFTCAMB error')
        ! the first parameter is B0:
        else if ( i==1 ) then
            latexname = trim('B_0')
            return
        ! the other parameters are the w_DE parameters:
        else
            call self%DesfRwDE%parameter_names_latex( i-1, latexname )
            return
        end if

    end subroutine eftcambdesignerfrparameternameslatex

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

        implicit none

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

        ! check validity of input:
        if ( i<=0 .or. i>self%parameter_number ) then
            write(*,'(a,I3)') 'EFTCAMB error: no parameter corresponding to number ', i
            write(*,'(a,I3)') 'Total number of parameters is ', self%parameter_number
            call mpistop('EFTCAMB error')
        ! the first parameter is B0:
        else if ( i==1 ) then
            value = self%B0
            return
        ! the other parameters are the w_DE parameters:
        else
            call self%DesfRwDE%parameter_value( i-1, value )
            return
        end if

    end subroutine eftcambdesignerfrparametervalues

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

        implicit none

        class(eftcamb_fr_designer)                   :: 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) :: x, mu
        integer  :: ind

        x   = log(a)
        call self%EFTOmega%precompute(x, ind, mu )

        eft_cache%EFTOmegaV    = self%EFTOmega%value( x, index=ind, coeff=mu )
        eft_cache%EFTOmegaP    = self%EFTOmega%first_derivative( x, index=ind, coeff=mu )
        eft_cache%EFTOmegaPP   = self%EFTOmega%second_derivative( x, index=ind, coeff=mu )
        eft_cache%EFTOmegaPPP  = self%EFTOmega%third_derivative( x, index=ind, coeff=mu )
        eft_cache%EFTc         = 0._dl
        eft_cache%EFTLambda    = self%EFTLambda%value( x, index=ind, coeff=mu )
        eft_cache%EFTcdot      = 0._dl
        eft_cache%EFTLambdadot = self%EFTLambda%first_derivative( x, index=ind, coeff=mu )

    end subroutine eftcambdesignerfrbackgroundeftfunctions

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

        implicit none

        class(eftcamb_fr_designer)                   :: 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%EFTGamma1V  = 0._dl
        eft_cache%EFTGamma1P  = 0._dl
        eft_cache%EFTGamma2V  = 0._dl
        eft_cache%EFTGamma2P  = 0._dl
        eft_cache%EFTGamma3V  = 0._dl
        eft_cache%EFTGamma3P  = 0._dl
        eft_cache%EFTGamma4V  = 0._dl
        eft_cache%EFTGamma4P  = 0._dl
        eft_cache%EFTGamma4PP = 0._dl
        eft_cache%EFTGamma5V  = 0._dl
        eft_cache%EFTGamma5P  = 0._dl
        eft_cache%EFTGamma6V  = 0._dl
        eft_cache%EFTGamma6P  = 0._dl

    end subroutine eftcambdesignerfrsecondordereftfunctions

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

        implicit none

        class(eftcamb_fr_designer)                   :: 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) :: eftcambdesignerfrcomputedtauda

        real(dl) :: temp

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

    end function eftcambdesignerfrcomputedtauda

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

        implicit none

        class(eftcamb_fr_designer)                   :: 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%DesfRwDE%integral(a)
        eft_cache%adotoa  = sqrt( ( eft_cache%grhom_t +eft_cache%grhov_t )/3._dl )

    end subroutine eftcambdesignerfrcomputeadotoa

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

        implicit none

        class(eftcamb_fr_designer)                   :: 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%DesfRwDE%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%DesfRwDE%value(a) +1.5_dl*self%DesfRwDE%value(a)**2 -0.5_dl*a*self%DesfRwDE%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 eftcambdesignerfrcomputehubbleder

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

        implicit none

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

        logical :: eftcambdesignerfradditionalmodelstability

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

    end function eftcambdesignerfradditionalmodelstability

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

end module eftcamb_designer_fr

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