EFTCAMB  Reference documentation for version 3.0
09p1_Designer_fR.f90
Go to the documentation of this file.
1 !----------------------------------------------------------------------------------------
2 !
3 ! This file is part of EFTCAMB.
4 !
5 ! Copyright (C) 2013-2016 by the EFTCAMB authors
6 !
7 ! The EFTCAMB code is free software;
8 ! You can use it, redistribute it, and/or modify it under the terms
9 ! of the GNU General Public License as published by the Free Software Foundation;
10 ! either version 3 of the License, or (at your option) any later version.
11 ! The full text of the license can be found in the file eftcamb/LICENSE at
12 ! the top level of the EFTCAMB distribution.
13 !
14 !----------------------------------------------------------------------------------------
15 
18 
19 
20 !----------------------------------------------------------------------------------------
22 
24 
26 
27  use precision
28  use inifile
29  use amlutils
31  use eft_def
33  use eftcamb_cache
42 
43  implicit none
44 
45  private
46 
47  public eftcamb_fr_designer
48 
49  !----------------------------------------------------------------------------------------
52  type, extends ( eftcamb_designer_model ) :: eftcamb_fr_designer
53 
54  ! theory parameters:
55  real(dl) :: b0
56 
57  ! the pure EFT functions model selection flags:
58  integer :: eftwde
59 
60  ! the pure EFT functions:
61  class( parametrized_function_1d ), allocatable :: desfrwde
62 
63  ! the interpolated EFT functions that come out of the background sover:
64  type(equispaced_linear_interpolate_function_1d) :: eftomega
65  type(equispaced_linear_interpolate_function_1d) :: eftlambda
66 
67  ! some designer parameters:
68  integer :: designer_num_points = 1000
69  real(dl) :: x_initial = log(10._dl**(-8._dl))
70  real(dl) :: x_final = 0.0_dl
71 
72  contains
73 
74  ! initialization of the model:
75  procedure :: read_model_selection => eftcambdesignerfrreadmodelselectionfromfile
76  procedure :: allocate_model_selection => eftcambdesignerfrallocatemodelselection
77  procedure :: init_model_parameters => eftcambdesignerfrinitmodelparameters
78  procedure :: init_model_parameters_from_file => eftcambdesignerfrinitmodelparametersfromfile
79 
80  ! background solver:
81  procedure :: initialize_background => eftcambdesignerfrinitbackground
82  procedure :: solve_designer_equations => eftcambdesignerfrsolvedesignerequations
83  procedure :: find_initial_conditions => eftcambdesignerfrfindinitialconditions
84 
85  ! utility functions:
86  procedure :: compute_param_number => eftcambdesignerfrcomputeparametersnumber
87  procedure :: feedback => eftcambdesignerfrfeedback
88  procedure :: parameter_names => eftcambdesignerfrparameternames
89  procedure :: parameter_names_latex => eftcambdesignerfrparameternameslatex
90  procedure :: parameter_values => eftcambdesignerfrparametervalues
91 
92  ! CAMB related procedures:
93  procedure :: compute_background_eft_functions => eftcambdesignerfrbackgroundeftfunctions
94  procedure :: compute_secondorder_eft_functions => eftcambdesignerfrsecondordereftfunctions
95  procedure :: compute_dtauda => eftcambdesignerfrcomputedtauda
96  procedure :: compute_adotoa => eftcambdesignerfrcomputeadotoa
97  procedure :: compute_h_derivs => eftcambdesignerfrcomputehubbleder
98 
99  ! stability procedures:
100  procedure :: additional_model_stability => eftcambdesignerfradditionalmodelstability
101 
102  end type eftcamb_fr_designer
103 
104  ! ---------------------------------------------------------------------------------------------
105 
106 contains
107 
108  ! ---------------------------------------------------------------------------------------------
110  subroutine eftcambdesignerfrreadmodelselectionfromfile( self, Ini )
111 
112  implicit none
113 
114  class(eftcamb_fr_designer) :: self
115  type(tinifile) :: ini
116 
117  ! read model selection flags:
118  self%EFTwDE = ini_read_int_file( ini, 'EFTwDE', 0 )
119 
120  end subroutine eftcambdesignerfrreadmodelselectionfromfile
121 
122  ! ---------------------------------------------------------------------------------------------
124  subroutine eftcambdesignerfrallocatemodelselection( self )
125 
126  implicit none
127 
128  class(eftcamb_fr_designer) :: self
129  character, allocatable, dimension(:) :: param_names
130  character, allocatable, dimension(:) :: param_names_latex
131 
132  ! allocate wDE:
133  if ( allocated(self%DesfRwDE) ) deallocate(self%DesfRwDE)
134  select case ( self%EFTwDE )
135  case(0)
136  allocate( wde_lcdm_parametrization_1d::self%DesfRwDE )
137  case(1)
138  allocate( constant_parametrization_1d::self%DesfRwDE )
139  case(2)
140  allocate( cpl_parametrization_1d::self%DesfRwDE )
141  call self%DesfRwDE%set_param_names( ['EFTw0', 'EFTwa'], ['w_0', 'w_a'] )
142  case(3)
143  allocate( jbp_parametrization_1d::self%DesfRwDE )
144  call self%DesfRwDE%set_param_names( ['EFTw0', 'EFTwa', 'EFTwn'], [ 'w_0', 'w_a', 'n ' ] )
145  case(4)
146  allocate( turning_point_parametrization_1d::self%DesfRwDE )
147  call self%DesfRwDE%set_param_names( ['EFTw0 ', 'EFTwa ', 'EFTwat'], ['w_0', 'w_a', 'a_t'] )
148  case(5)
149  allocate( taylor_parametrization_1d::self%DesfRwDE )
150  call self%DesfRwDE%set_param_names( ['EFTw0', 'EFTwa', 'EFTw2', 'EFTw3'], ['w_0', 'w_a', 'w_2', 'w_3'] )
151  case default
152  write(*,'(a,I3)') 'No model corresponding to EFTwDE =', self%EFTwDE
153  write(*,'(a)') 'Please select an appropriate model.'
154  end select
155 
156  ! initialize the names:
157  call self%DesfRwDE%set_name( 'EFTw', 'w' )
158 
159  end subroutine eftcambdesignerfrallocatemodelselection
160 
161  ! ---------------------------------------------------------------------------------------------
163  subroutine eftcambdesignerfrinitmodelparameters( self, array )
164 
165  implicit none
166 
167  class(eftcamb_fr_designer) :: self
168  real(dl), dimension(self%parameter_number), intent(in) :: array
169  real(dl), dimension(self%parameter_number -1) :: temp
170  integer :: i
171 
172  self%B0 = array(1)
173 
174  do i = 1, self%parameter_number -1
175  temp(i) = array(i+1)
176  end do
177  call self%DesfRwDE%init_parameters(temp)
178 
179  end subroutine eftcambdesignerfrinitmodelparameters
180 
181  ! ---------------------------------------------------------------------------------------------
183  subroutine eftcambdesignerfrinitmodelparametersfromfile( self, Ini )
184 
185  implicit none
186 
187  class(eftcamb_fr_designer) :: self
188  type(tinifile) :: ini
189 
190  ! read B0:
191  self%B0 = ini_read_double_file( ini, 'EFTB0', 0._dl )
192  ! read w_DE parameters:
193  call self%DesfRwDE%init_from_file( ini )
194 
195  end subroutine eftcambdesignerfrinitmodelparametersfromfile
196 
197  ! ---------------------------------------------------------------------------------------------
199  subroutine eftcambdesignerfrinitbackground( self, params_cache, feedback_level, success )
200 
201  implicit none
202 
203  class(eftcamb_fr_designer) :: self
204  type(eftcamb_parameter_cache), intent(in) :: params_cache
205  integer , intent(in) :: feedback_level
206  logical , intent(out) :: success
207 
208  real(dl) :: a_ini, b0
209  real(dl) :: tempmin, tempmax, debug_a
210  integer :: debug_maxnum, debug_n
211 
212  ! some feedback:
213  if ( feedback_level>0 ) then
214  write(*,'(a)') "***************************************************************"
215  write(*,'(a)') ' EFTCAMB designer f(R) background solver'
216  write(*,'(a)')
217  end if
218 
219  ! initialize interpolating functions:
220  call self%EFTOmega%initialize ( self%designer_num_points, self%x_initial, self%x_final )
221  call self%EFTLambda%initialize ( self%designer_num_points, self%x_initial, self%x_final )
222 
223  ! debug code:
224  if ( debugeftcamb ) then
225  ! print the function B0(A). This is used to debug the initial conditions part.
226  call createtxtfile( './debug_designer_fR_B.dat', 34 )
227  print*, 'EFTCAMB DEBUG ( f(R) designer ): Printing B(A) results'
228  tempmin = -10._dl
229  tempmax = +10._dl
230  debug_maxnum = 1000
231  do debug_n = 1, debug_maxnum
232  debug_a = tempmin +REAL(debug_n-1)*(tempmax-tempmin)/REAL(debug_maxnum-1)
233  call self%solve_designer_equations( params_cache, debug_a, b0, only_b0=.true., success=success )
234  write(34,*) debug_a, b0
235  end do
236  close(34)
237  ! prints f(R) quantities.
238  print*, 'EFTCAMB DEBUG ( f(R) designer ): Printing F(R) results'
239  call createtxtfile( './debug_designer_fR_solution.dat', 33 )
240  debug_a = 1.0_dl
241  call self%solve_designer_equations( params_cache, debug_a, b0, only_b0=.false., success=success )
242  close(33)
243  end if
244 
245  ! call boundary conditions lookup:
246  call self%find_initial_conditions( params_cache, feedback_level, a_ini, success )
247 
248  ! solve the background equations and store the solution:
249  call self%solve_designer_equations( params_cache, a_ini, b0, only_b0=.false., success=success )
250 
251  end subroutine eftcambdesignerfrinitbackground
252 
253  ! ---------------------------------------------------------------------------------------------
255  subroutine eftcambdesignerfrsolvedesignerequations( self, params_cache, A, B0, only_B0, success )
256 
257  implicit none
258 
259  class(eftcamb_fr_designer) :: self
260  type(eftcamb_parameter_cache), intent(in) :: params_cache
261  real(dl), intent(in) :: a
262  real(dl), intent(out) :: b0
263  logical , optional :: only_b0
264  logical , intent(out) :: success
265 
266  integer, parameter :: num_eq = 2
267 
268  real(dl) :: omegam_eft, omegavac_eft, omegamassivenu_eft, omegagamma_eft, omeganu_eft
269  real(dl) :: omegarad_eft, equivalencescale_fr, ratio_fr, initial_b_fr, initial_c_fr
270  real(dl) :: pplus, yplus, coeffa_part, ystar, x
271 
272  real(dl) :: y(num_eq), ydot(num_eq)
273 
274  integer :: itol, itask, istate, iopt, lrn, lrs, lrw, lis, lin, liw, jacobianmode, i
275  real(dl) :: rtol, atol, t1, t2, b
276  real(dl), allocatable :: rwork(:)
277  integer, allocatable :: iwork(:)
278 
279  ! digedt the input:
280  if ( .not. present(only_b0) ) only_b0 = .false.
281 
282  ! 1) Cosmological parameters:
283  omegam_eft = params_cache%omegab + params_cache%omegac
284  omegavac_eft = params_cache%omegav
285  omegamassivenu_eft = params_cache%omegan
286  omegagamma_eft = params_cache%omegag
287  omeganu_eft = params_cache%omegar
288 
289  omegarad_eft = omegagamma_eft + omeganu_eft
290 
291  equivalencescale_fr = (omegarad_eft+ omegamassivenu_eft)/omegam_eft
292  ratio_fr = equivalencescale_fr/exp( self%x_initial )
293 
294  ! 2) Growing mode solution:
295  initial_b_fr = (7._dl + 8._dl*ratio_fr)/(2._dl*(1._dl + ratio_fr))
296  initial_c_fr = -3._dl/(2._dl*(1._dl + ratio_fr))
297  pplus = 0.5_dl*(-initial_b_fr + sqrt(initial_b_fr**2 - 4._dl*initial_c_fr))
298  yplus = exp(pplus*self%x_initial)
299  ! Construction of the particolar solution:
300  coeffa_part = (-6._dl*initial_c_fr)/(-3._dl*exp(self%x_initial)*self%DesfRwDE%first_derivative( exp(self%x_initial) )&
301  & +9._dl*self%DesfRwDE%value( exp(self%x_initial) )**2&
302  & +(18._dl-3._dl*initial_b_fr)*self%DesfRwDE%value( exp(self%x_initial) )&
303  & +9._dl -3._dl*initial_b_fr +initial_c_fr)
304  ystar = coeffa_part*omegavac_eft*exp(-2._dl*self%x_initial)*self%DesfRwDE%integral( exp(self%x_initial) )
305 
306  ! 3) Set initial conditions:
307  x = self%x_initial
308  y(1) = a*yplus + ystar
309  y(2) = pplus*a*yplus - 3._dl*( 1._dl+ self%DesfRwDE%value(exp(self%x_initial)) )*ystar
310  ydot = 0._dl
311 
312  ! 4) Initialize DLSODA:
313  ! set-up the relative and absolute tollerances:
314  itol = 1
315  rtol = 1.d-10
316  atol = 1.d-14
317  ! initialize task to do:
318  itask = 1
319  istate = 1
320  iopt = 1
321  ! initialize the work space:
322  lrn = 20 + 16*num_eq
323  lrs = 22 + 9*num_eq + num_eq**2
324  lrw = max(lrn,lrs)
325  lis = 20 + num_eq
326  lin = 20
327  liw = max(lis,lin)
328  ! allocate the arrays:
329  allocate(rwork(lrw))
330  allocate(iwork(liw))
331  ! optional lsoda input:
332  rwork(5) = 0._dl ! the step size to be attempted on the first step. The default value is determined by the solver.
333  rwork(6) = 0._dl ! the maximum absolute step size allowed. The default value is infinite.
334  rwork(7) = 0._dl ! the minimum absolute step size allowed. The default value is 0.
335  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.
336  iwork(6) = 100 ! maximum number of (internally defined) steps allowed during one call to the solver. The default value is 500.
337  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.
338  iwork(8) = 0 ! the maximum order to be allowed for the nonstiff (Adams) method. the default value is 12.
339  iwork(9) = 0 ! the maximum order to be allowed for the stiff (BDF) method. The default value is 5.
340  ! additional lsoda stuff:
341  CALL xsetf(0) ! suppress odepack printing
342  ! Jacobian mode: 1=fullJacobian, 2=not provided
343  jacobianmode = 2
344 
345  ! store the first value of the EFT functions:
346  t1 = self%EFTOmega%x(1)
347  ! compute output EFT functions if needed:
348  if ( .not. only_b0 ) then
349  call output( num_eq, 1, t1, y, b0 )
350  end if
351 
352  ! 5) solve the equations:
353  do i=1, self%EFTOmega%num_points-1
354 
355  ! set the time step:
356  t1 = self%EFTOmega%x(i)
357  t2 = self%EFTOmega%x(i+1)
358  ! solve the system:
359  call dlsoda ( derivs, num_eq, y, t1, t2, itol, rtol, atol, itask, istate, iopt, rwork, lrw, iwork, liw, jacobian, jacobianmode)
360  ! check istate for LSODA good completion:
361  if ( istate < 0 ) then
362  success = .false.
363  return
364  end if
365  ! compute output EFT functions if needed:
366  if ( .not. only_b0 ) then
367  call output( num_eq, i+1, t2, y, b0 )
368  end if
369 
370  end do
371 
372  ! compute B0:
373  if ( only_b0 ) then
374  call output( num_eq, self%EFTOmega%num_points, t2, y, b0 )
375  end if
376 
377  return
378 
379  contains
380 
381  ! ---------------------------------------------------------------------------------------------
383  subroutine derivs( num_eq, x, y, ydot )
384 
385  implicit none
386 
387  integer , intent(in) :: num_eq
388  real(dl), intent(in) :: x
389  real(dl), intent(in) , dimension(num_eq) :: y
390  real(dl), intent(out), dimension(num_eq) :: ydot
391 
392  real(dl) :: a, eft_e_gfun, eft_e_gfunp, eft_e_gfunpp, eft_e_gfunppp
393  real(dl) :: rhonu_tot, presnu_tot, presnudot_tot, presnudotdot_tot
394  real(dl) :: eft_e_nu, eft_ep_nu, eft_epp_nu, eft_e3p_nu, rhonu, presnu, grhormass_t
395  real(dl) :: efunction, efunprime, adotoa, hdot, presnudot, presnudotdot
396  real(dl) :: efunprime2, efunprime3
397  integer :: nu_i
398 
399  ! 1) convert x in a:
400  a = exp(x)
401 
402  ! 2) Compute the function g(x) and its derivatives:
403  eft_e_gfun = -(log( self%DesfRwDE%integral(a) ) -2._dl*x)/3._dl
404  eft_e_gfunp = 1._dl +self%DesfRwDE%value(a)
405  eft_e_gfunpp = exp(x)*self%DesfRwDE%first_derivative(a)
406  eft_e_gfunppp = exp(x)*self%DesfRwDE%first_derivative(a) +exp(2._dl*x)*self%DesfRwDE%second_derivative(a)
407 
408  ! 3) Compute energy and its derivatives:
409  ! First compute massive neutrinos contribution:
410  rhonu_tot = 0._dl
411  presnu_tot = 0._dl
412  eft_e_nu = 0._dl
413  eft_ep_nu = 0._dl
414  if ( params_cache%Num_Nu_Massive /= 0) then
415  do nu_i = 1, params_cache%Nu_mass_eigenstates
416 
417  rhonu = 0._dl
418  presnu = 0._dl
419  grhormass_t= params_cache%grhormass(nu_i)/a**2
420  call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
421  rhonu_tot = rhonu_tot + grhormass_t*rhonu
422  presnu_tot = presnu_tot + grhormass_t*presnu
423 
424  eft_e_nu = eft_e_nu + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
425  eft_ep_nu = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)
426 
427  end do
428  end if
429 
430  ! Add its contribution to E and E':
431  efunction = +omegarad_eft*exp(-4._dl*x)&
432  & +omegam_eft*exp(-3._dl*x)&
433  & +omegavac_eft*exp(-3._dl*eft_e_gfun) + eft_e_nu
434  efunprime = -4._dl*omegarad_eft*exp(-4._dl*x)&
435  & -3._dl*omegam_eft*exp(-3._dl*x)&
436  & -3._dl*omegavac_eft*eft_e_gfunp*exp(-3._dl*eft_e_gfun) +eft_ep_nu
437 
438  ! Compute everything of massive nu again to get the time derivatives:
439  rhonu_tot = 0._dl
440  presnu_tot = 0._dl
441  presnudot_tot = 0._dl
442  presnudotdot_tot = 0._dl
443  eft_e_nu = 0._dl
444  eft_ep_nu = 0._dl
445  eft_epp_nu = 0._dl
446  eft_e3p_nu = 0._dl
447  if ( params_cache%Num_Nu_Massive /= 0 ) then
448  do nu_i = 1, params_cache%Nu_mass_eigenstates
449 
450  adotoa = +a*params_cache%h0_Mpc*sqrt(efunction)
451  hdot = +0.5_dl*params_cache%h0_Mpc**2*a**2*efunprime +adotoa**2
452 
453  rhonu = 0._dl
454  presnu = 0._dl
455  presnudot = 0._dl
456  presnudotdot = 0._dl
457 
458  grhormass_t = params_cache%grhormass(nu_i)/a**2
459 
460  call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
461  presnudot = params_cache%Nu_pidot(a*params_cache%nu_masses(nu_i),adotoa,presnu)
462  presnudotdot = params_cache%Nu_pidotdot(a*params_cache%nu_masses(nu_i),adotoa,hdot,presnu,presnudot)
463 
464  rhonu_tot = rhonu_tot + grhormass_t*rhonu
465  presnu_tot = presnu_tot + grhormass_t*presnu
466  presnudot_tot = presnudot_tot + grhormass_t*(presnudot -4._dl*adotoa*presnu)
467  presnudotdot_tot = presnudotdot_tot + grhormass_t*(presnudotdot &
468  & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))
469 
470  eft_e_nu = eft_e_nu + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
471  eft_ep_nu = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)
472  eft_epp_nu = eft_epp_nu + 3._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
473  & -grhormass_t*(presnudot -4._dl*adotoa*presnu)/params_cache%h0_Mpc**3/sqrt(efunction)/a**3
474  eft_e3p_nu = eft_e3p_nu -9._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
475  & +(3._dl/adotoa/params_cache%h0_Mpc**2/a**2+hdot/adotoa**3/params_cache%h0_Mpc**2/a**2)&
476  &*grhormass_t*(presnudot -4._dl*adotoa*presnu)&
477  & -grhormass_t*(presnudotdot &
478  & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))/adotoa**2/params_cache%h0_Mpc**2/a**2
479  end do
480  end if
481 
482  efunprime2 = 16._dl*omegarad_eft*exp(-4._dl*x)&
483  & +9._dl*omegam_eft*exp(-3._dl*x)&
484  & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(eft_e_gfunpp -3._dl*eft_e_gfunp**2) + eft_epp_nu
485  efunprime3 = -64._dl*omegarad_eft*exp(-4._dl*x)&
486  & -27._dl*omegam_eft*exp(-3._dl*x)&
487  & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*&
488  &(eft_e_gfunppp-9._dl*eft_e_gfunp*eft_e_gfunpp+9._dl*eft_e_gfunp**3) + eft_e3p_nu
489 
490  ! 4) Get the equation of motion:
491  ydot(1) = y(2)
492  ydot(2) = (1._dl+0.5_dl*efunprime/efunction+(4._dl*efunprime2+efunprime3)/(4._dl*efunprime+efunprime2))*y(2) &
493  & -0.5_dl*(4._dl*efunprime+efunprime2)/efunction*y(1) &
494  & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(4._dl*efunprime+efunprime2)/efunction
495 
496  end subroutine
497 
498  ! ---------------------------------------------------------------------------------------------
501  subroutine jacobian( num_eq, x, y, ml, mu, pd, nrowpd )
502 
503  implicit none
504 
505  integer :: num_eq
506  integer :: ml
507  integer :: mu
508  integer :: nrowpd
509  real(dl) :: x
510  real(dl), dimension(num_eq) :: y
511  real(dl), dimension(nrowpd,num_eq) :: pd
512 
513  end subroutine jacobian
514 
515  ! ---------------------------------------------------------------------------------------------
518  subroutine output( num_eq, ind, x, y, B )
519 
520  implicit none
521 
522  integer , intent(in) :: num_eq
523  integer , intent(in) :: ind
524  real(dl), intent(in) :: x
525  real(dl), intent(in) , dimension(num_eq) :: y
526  real(dl), intent(out) :: b
527 
528  real(dl) :: ydot(num_eq)
529  real(dl) :: a, eft_e_gfun, eft_e_gfunp, eft_e_gfunpp, eft_e_gfunppp
530  real(dl) :: rhonu_tot, presnu_tot, presnudot_tot, presnudotdot_tot
531  real(dl) :: eft_e_nu, eft_ep_nu, eft_epp_nu, eft_e3p_nu, rhonu, presnu, grhormass_t
532  real(dl) :: efunction, efunprime, adotoa, hdot, presnudot, presnudotdot
533  real(dl) :: efunprime2, efunprime3, ricci, f_sub_r, ede, edep, edepp
534  integer :: nu_i
535  logical :: is_open
536 
537  ! 1) convert x in a:
538  a = exp(x)
539 
540  ! 2) Compute the function g(x) and its derivatives:
541  eft_e_gfun = -(log( self%DesfRwDE%integral(a) ) -2._dl*x)/3._dl
542  eft_e_gfunp = 1._dl +self%DesfRwDE%value(a)
543  eft_e_gfunpp = exp(x)*self%DesfRwDE%first_derivative(a)
544  eft_e_gfunppp = exp(x)*self%DesfRwDE%first_derivative(a) +exp(2._dl*x)*self%DesfRwDE%second_derivative(a)
545 
546  ! 3) Compute energy and its derivatives:
547  ! First compute massive neutrinos contribution:
548  rhonu_tot = 0._dl
549  presnu_tot = 0._dl
550  eft_e_nu = 0._dl
551  eft_ep_nu = 0._dl
552  if ( params_cache%Num_Nu_Massive /= 0) then
553  do nu_i = 1, params_cache%Nu_mass_eigenstates
554 
555  rhonu = 0._dl
556  presnu = 0._dl
557  grhormass_t= params_cache%grhormass(nu_i)/a**2
558  call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
559  rhonu_tot = rhonu_tot + grhormass_t*rhonu
560  presnu_tot = presnu_tot + grhormass_t*presnu
561 
562  eft_e_nu = eft_e_nu + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
563  eft_ep_nu = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)
564 
565  end do
566  end if
567 
568  ! Add its contribution to E and E':
569  efunction = +omegarad_eft*exp(-4._dl*x)&
570  & +omegam_eft*exp(-3._dl*x)&
571  & +omegavac_eft*exp(-3._dl*eft_e_gfun) + eft_e_nu
572  efunprime = -4._dl*omegarad_eft*exp(-4._dl*x)&
573  & -3._dl*omegam_eft*exp(-3._dl*x)&
574  & -3._dl*omegavac_eft*eft_e_gfunp*exp(-3._dl*eft_e_gfun) +eft_ep_nu
575 
576  ! Compute everything of massive nu again to get the time derivatives:
577  rhonu_tot = 0._dl
578  presnu_tot = 0._dl
579  presnudot_tot = 0._dl
580  presnudotdot_tot = 0._dl
581  eft_e_nu = 0._dl
582  eft_ep_nu = 0._dl
583  eft_epp_nu = 0._dl
584  eft_e3p_nu = 0._dl
585  if ( params_cache%Num_Nu_Massive /= 0 ) then
586  do nu_i = 1, params_cache%Nu_mass_eigenstates
587 
588  adotoa = +a*params_cache%h0_Mpc*sqrt(efunction)
589  hdot = +0.5_dl*params_cache%h0_Mpc**2*a**2*efunprime +adotoa**2
590 
591  rhonu = 0._dl
592  presnu = 0._dl
593  presnudot = 0._dl
594  presnudotdot = 0._dl
595 
596  grhormass_t = params_cache%grhormass(nu_i)/a**2
597 
598  call params_cache%Nu_background(a*params_cache%nu_masses(nu_i),rhonu,presnu)
599  presnudot = params_cache%Nu_pidot(a*params_cache%nu_masses(nu_i),adotoa,presnu)
600  presnudotdot = params_cache%Nu_pidotdot(a*params_cache%nu_masses(nu_i),adotoa,hdot,presnu,presnudot)
601 
602  rhonu_tot = rhonu_tot + grhormass_t*rhonu
603  presnu_tot = presnu_tot + grhormass_t*presnu
604  presnudot_tot = presnudot_tot + grhormass_t*(presnudot -4._dl*adotoa*presnu)
605  presnudotdot_tot = presnudotdot_tot + grhormass_t*(presnudotdot &
606  & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))
607 
608  eft_e_nu = eft_e_nu + params_cache%grhormass(nu_i)/3._dl/a**4/params_cache%h0_Mpc**2*rhonu
609  eft_ep_nu = eft_ep_nu - params_cache%grhormass(nu_i)/params_cache%h0_Mpc**2/a**4*(rhonu +presnu)
610  eft_epp_nu = eft_epp_nu + 3._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
611  & -grhormass_t*(presnudot -4._dl*adotoa*presnu)/params_cache%h0_Mpc**3/sqrt(efunction)/a**3
612  eft_e3p_nu = eft_e3p_nu -9._dl/params_cache%h0_Mpc**2*params_cache%grhormass(nu_i)/a**4*(rhonu +presnu)&
613  & +(3._dl/adotoa/params_cache%h0_Mpc**2/a**2+hdot/adotoa**3/params_cache%h0_Mpc**2/a**2)&
614  &*grhormass_t*(presnudot -4._dl*adotoa*presnu)&
615  & -grhormass_t*(presnudotdot &
616  & -8._dl*adotoa*presnudot +4._dl*presnu*(+4._dl*adotoa**2-hdot))/adotoa**2/params_cache%h0_Mpc**2/a**2
617  end do
618  end if
619 
620  efunprime2 = 16._dl*omegarad_eft*exp(-4._dl*x)&
621  & +9._dl*omegam_eft*exp(-3._dl*x)&
622  & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(eft_e_gfunpp -3._dl*eft_e_gfunp**2) + eft_epp_nu
623  efunprime3 = -64._dl*omegarad_eft*exp(-4._dl*x)&
624  & -27._dl*omegam_eft*exp(-3._dl*x)&
625  & -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*&
626  &(eft_e_gfunppp-9._dl*eft_e_gfunp*eft_e_gfunpp+9._dl*eft_e_gfunp**3) + eft_e3p_nu
627 
628  ! compute E dark energy:
629  ede = +omegavac_eft*exp(-3._dl*eft_e_gfun)
630  edep = -3._dl*omegavac_eft*eft_e_gfunp*exp(-3._dl*eft_e_gfun)
631  edepp = -3._dl*omegavac_eft*exp(-3._dl*eft_e_gfun)*(eft_e_gfunpp -3._dl*eft_e_gfunp**2)
632 
633  ! call derivs to compute ydot at a given time:
634  call derivs( num_eq, x, y, ydot )
635 
636  ! Compute the Ricci:
637  ricci = 3._dl*(4._dl*efunction+efunprime)
638  f_sub_r = ydot(1)/3._dl/(4._dl*efunprime +efunprime2)
639 
640  ! compute B:
641  b = 2._dl/3._dl/(1._dl +f_sub_r)*1._dl/(4._dl*efunprime +efunprime2)*efunction/efunprime*&
642  &( ydot(2) -ydot(1)*(4._dl*efunprime2+efunprime3)/(4._dl*efunprime +efunprime2) )
643 
644  ! compute the EFT functions:
645  self%EFTOmega%y(ind) = f_sub_r
646  self%EFTOmega%yp(ind) = exp(-x)/(6._dl*efunction*(efunprime2+4._dl*efunprime))*(-efunprime2*(6._dl*ede+y(1))&
647  &+efunprime*(ydot(1)-4._dl*(6._dl*ede +y(1)))+2._dl*efunction*ydot(1))
648  self%EFTOmega%ypp(ind) = exp(-2._dl*x)/(12._dl*efunction**2*(efunprime2+4._dl*efunprime))*&
649  &(efunprime*(efunprime*(24._dl*ede-ydot(1) +4._dl*y(1))-6._dl*efunction*(8._dl*edep +ydot(1)))&
650  &+efunprime2*(efunprime*(6._dl*ede +y(1))-12._dl*efunction*edep))
651  self%EFTOmega%yppp(ind) = exp(-3._dl*x)/(24._dl*efunction**3*(efunprime2+4._dl*efunprime))*&
652  &(2._dl*efunction*(efunprime2**2*(6._dl*ede+y(1))&
653  & +4._dl*efunprime**2*(18._dl*edep+6._dl*ede+2._dl*ydot(1) +y(1))&
654  & +efunprime*efunprime2*(18._dl*edep+30._dl*ede -ydot(1)+5._dl*y(1)))&
655  & +12._dl*efunction**2*(efunprime2*(-2._dl*edepp+4._dl*edep-ydot(1))&
656  & +efunprime*(-8._dl*edepp+16._dl*edep +ydot(1)))&
657  & +3._dl*efunprime**2*(efunprime*(ydot(1)-4._dl*(6._dl*ede +y(1)))-efunprime2*(6._dl*ede +y(1))))
658  self%EFTLambda%y(ind) = 0.5_dl*params_cache%h0_Mpc**2*( y(1) -f_sub_r*ricci )*exp(x)**2
659  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))
660 
661  if ( debugeftcamb ) then
662  inquire( unit=33, opened=is_open )
663  if ( is_open ) then
664  write (33,'(20E15.5)') x, exp(x), ricci, y(1), &
665  & self%EFTOmega%y(ind), self%EFTOmega%yp(ind), self%EFTOmega%ypp(ind), self%EFTOmega%yppp(ind), self%EFTLambda%y(ind), self%EFTLambda%yp(ind),&
666  & b
667  end if
668  end if
669 
670  end subroutine
671 
672  end subroutine eftcambdesignerfrsolvedesignerequations
673 
674  ! ---------------------------------------------------------------------------------------------
677  subroutine eftcambdesignerfrfindinitialconditions( self, params_cache, feedback_level, A_ini, success )
678 
679  implicit none
680 
681  class(eftcamb_fr_designer) :: self
682  type(eftcamb_parameter_cache), intent(in) :: params_cache
683  integer , intent(in) :: feedback_level
684  real(dl) , intent(out) :: a_ini
685  logical , intent(out) :: success
686 
687  real(dl) :: atemp1, atemp2, btemp1, btemp2, horizasyntb
688  real(dl) :: vertasynta, atemp3, atemp4, btemp3, btemp4, realap
689  integer :: ind
690 
691  ! Find initial conditions for designer F(R).
692  ! Initial conditions are found solving B0(A) = B0wanted.
693  ! The next part of code is complicated by the fact that B0(A) is not continuous and we have to adopt ad-hoc strategies.
694 
695  ! 1) Find the horizontal asymptote of B0(A)
696  atemp1 = 0._dl
697  atemp2 = 0._dl
698  do ind=10, 100, 1
699  atemp1 = atemp2
700  atemp2 = 10._dl**REAL(ind)
701  btemp1 = desfr_bfunca(atemp1)
702  btemp2 = desfr_bfunca(atemp2)
703  if (abs(btemp1-btemp2)/abs(atemp1-atemp2).lt.1.d-100) exit
704  end do
705  horizasyntb = btemp2
706 
707  if ( feedback_level>1 ) write(*,'(a,E13.4)') ' horizontal asymptote = ', horizasyntb
708 
709  ! 2) Check that the value of B0 given is not the forbidden one: the one corresponding to the horizontal asynt.
710  if ( abs(horizasyntb-self%B0)<1.d-15 ) then
711  success=.false.
712  return
713  end if
714 
715  ! 3) Bracket the vertical asyntote
716  atemp1 = -10._dl
717  atemp2 = 10._dl
718  call zbrac(desfr_bfunca,atemp1,atemp2,success,horizasyntb)
719  if (.not.success) then
720  if ( feedback_level>1 ) write(*,'(a)') ' FAILURE of vertical asymptote bracketing'
721  return
722  end if
723 
724  ! 4) Find the vertical asyntote by tricking the root finding algorithm.
725  vertasynta = zbrent(desfr_bfunca,atemp1,atemp2,1.d-100,horizasyntb,success)
726  if (.not.success) then
727  if ( feedback_level>1 ) write(*,'(a)') ' FAILURE of vertical asyntote finding'
728  return
729  end if
730 
731  if ( feedback_level>1 ) write(*,'(a,E13.4)') ' vertical asymptote = ', vertasynta
732 
733  ! 5) Find values for A that are on the left and on the right of the asyntote.
734  do ind=-10, -1, 1
735  atemp1 = vertasynta+10._dl**ind
736  atemp2 = vertasynta-10._dl**ind
737  btemp1 = desfr_bfunca(atemp1)
738  btemp2 = desfr_bfunca(atemp2)
739  if (btemp1*btemp2<0._dl) exit
740  end do
741 
742  ! 6) Extablish on which side of the asyntote to look for the solution.
743  atemp3 = atemp1
744  atemp4 = atemp2
745  do ind=1, 10, 1
746  atemp3 = atemp3 + 10._dl**ind
747  atemp4 = atemp4 - 10._dl**ind
748  btemp3 = desfr_bfunca(atemp3)
749  btemp4 = desfr_bfunca(atemp4)
750  if ((btemp1-self%B0)*(btemp3-self%B0)<0.or.(btemp2-self%B0)*(btemp4-self%B0)<0) exit
751  end do
752  if ((btemp1-self%B0)*(btemp3-self%B0)<0.and.(btemp2-self%B0)*(btemp4-self%B0)<0) then
753  if ( feedback_level>1 ) write(*,'(a)') ' FAILURE as the root seems to be on both sides'
754  return
755  end if
756 
757  ! 7) Solve the equation B0(A)=B0_wanted.
758  if ((btemp1-self%B0)*(btemp3-self%B0)<0) then
759  realap = zbrent(desfr_bfunca,atemp1,atemp3,1.d-50,self%B0,success)
760  if (.not.success) then
761  if ( feedback_level>1 ) write(*,'(a)') ' FAILURE right side solution not found'
762  return
763  end if
764  else if ((btemp2-self%B0)*(btemp4-self%B0)<0) then
765  realap = zbrent(desfr_bfunca,atemp2,atemp4,1.d-50,self%B0,success)
766  if (.not.success) then
767  if ( feedback_level>1 ) write(*,'(a)') ' FAILURE left side solution not found'
768  return
769  end if
770  else
771  if ( feedback_level>1 ) write(*,'(a)') ' FAILURE the root was not on the right side nor the left one...'
772  return
773  end if
774 
775  if ( feedback_level>1 ) write(*,'(a,E13.4)') ' initial condition A = ', realap
776 
777  ! 8) Check if the result found is compatible with the requested one. This is required only for debug.
778  btemp1 = desfr_bfunca(realap)
779  if ( feedback_level>1 ) then
780  write(*,'(a,E13.4)') ' B0 found = ', btemp1
781  write(*,'(a,E13.4)') ' B0 given = ', self%B0
782  end if
783 
784  if (abs(btemp1-self%B0)/abs(self%B0)>0.1_dl.and.abs(btemp1-self%B0)>1.d-8.and..false.) then
785  if ( feedback_level>1 ) write(*,'(a)') ' FAILURE designer code unable to find appropriate initial conditions'
786  success = .false.
787  return
788  end if
789 
790  ! 9) output the value:
791  a_ini = realap
792 
793  return
794 
795  contains
796 
797  ! ---------------------------------------------------------------------------------------------
799  function desfr_bfunca(A)
800 
801  implicit none
802 
803  real(dl) :: a
804  real(dl) :: desfr_bfunca
805  real(dl) :: b0
806 
807  call self%solve_designer_equations( params_cache, a, b0, only_b0=.true., success=success )
808 
809  desfr_bfunca = b0
810 
811  end function
812 
813  end subroutine eftcambdesignerfrfindinitialconditions
814 
815  ! ---------------------------------------------------------------------------------------------
817  subroutine eftcambdesignerfrcomputeparametersnumber( self )
818 
819  implicit none
820 
821  class(eftcamb_fr_designer) :: self
822 
823  self%parameter_number = 1
824  self%parameter_number = self%parameter_number +self%DesfRwDE%parameter_number
825 
826  end subroutine eftcambdesignerfrcomputeparametersnumber
827 
828  ! ---------------------------------------------------------------------------------------------
830  subroutine eftcambdesignerfrfeedback( self, print_params )
831 
832  implicit none
833 
834  class(eftcamb_fr_designer) :: self
835  logical, optional :: print_params
837 
838  write(*,*)
839  write(*,'(a,a)') ' Model = ', self%name
840  write(*,'(a,I3)') ' Number of params =' , self%parameter_number
841  ! print model functions informations:
842  if ( self%EFTwDE /= 0 ) then
843  write(*,*)
844  write(*,'(a,I3)') ' EFTwDE =', self%EFTwDE
845  end if
846  write(*,*)
847  write(*,'(a24,F12.6)') ' B0 =', self%B0
848 
849  call self%DesfRwDE%feedback( print_params )
850 
851  end subroutine eftcambdesignerfrfeedback
852 
853  ! ---------------------------------------------------------------------------------------------
855  subroutine eftcambdesignerfrparameternames( self, i, name )
856 
857  implicit none
858 
859  class(eftcamb_fr_designer) :: self
860  integer , intent(in) :: i
861  character(*), intent(out) :: name
862 
863  ! check validity of input:
864  if ( i<=0 .or. i>self%parameter_number ) then
865  write(*,'(a,I3)') 'EFTCAMB error: no parameter corresponding to number ', i
866  write(*,'(a,I3)') 'Total number of parameters is ', self%parameter_number
867  call mpistop('EFTCAMB error')
868  ! the first parameter is B0:
869  else if ( i==1 ) then
870  name = trim('B0')
871  return
872  ! the other parameters are the w_DE parameters:
873  else
874  call self%DesfRwDE%parameter_names( i-1, name )
875  return
876  end if
877 
878  end subroutine eftcambdesignerfrparameternames
879 
880  ! ---------------------------------------------------------------------------------------------
882  subroutine eftcambdesignerfrparameternameslatex( self, i, latexname )
883 
884  implicit none
885 
886  class(eftcamb_fr_designer) :: self
887  integer , intent(in) :: i
888  character(*), intent(out) :: latexname
889 
890  ! check validity of input:
891  if ( i<=0 .or. i>self%parameter_number ) then
892  write(*,'(a,I3)') 'EFTCAMB error: no parameter corresponding to number ', i
893  write(*,'(a,I3)') 'Total number of parameters is ', self%parameter_number
894  call mpistop('EFTCAMB error')
895  ! the first parameter is B0:
896  else if ( i==1 ) then
897  latexname = trim('B_0')
898  return
899  ! the other parameters are the w_DE parameters:
900  else
901  call self%DesfRwDE%parameter_names_latex( i-1, latexname )
902  return
903  end if
904 
905  end subroutine eftcambdesignerfrparameternameslatex
906 
907  ! ---------------------------------------------------------------------------------------------
909  subroutine eftcambdesignerfrparametervalues( self, i, value )
910 
911  implicit none
912 
913  class(eftcamb_fr_designer) :: self
914  integer , intent(in) :: i
915  real(dl), intent(out) :: value
916 
917  ! check validity of input:
918  if ( i<=0 .or. i>self%parameter_number ) then
919  write(*,'(a,I3)') 'EFTCAMB error: no parameter corresponding to number ', i
920  write(*,'(a,I3)') 'Total number of parameters is ', self%parameter_number
921  call mpistop('EFTCAMB error')
922  ! the first parameter is B0:
923  else if ( i==1 ) then
924  value = self%B0
925  return
926  ! the other parameters are the w_DE parameters:
927  else
928  call self%DesfRwDE%parameter_value( i-1, value )
929  return
930  end if
931 
932  end subroutine eftcambdesignerfrparametervalues
933 
934  ! ---------------------------------------------------------------------------------------------
936  subroutine eftcambdesignerfrbackgroundeftfunctions( self, a, eft_par_cache, eft_cache )
937 
938  implicit none
939 
940  class(eftcamb_fr_designer) :: self
941  real(dl), intent(in) :: a
942  type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
943  type(eftcamb_timestep_cache ), intent(inout) :: eft_cache
944 
945  real(dl) :: x, mu
946  integer :: ind
947 
948  x = log(a)
949  call self%EFTOmega%precompute(x, ind, mu )
950 
951  eft_cache%EFTOmegaV = self%EFTOmega%value( x, index=ind, coeff=mu )
952  eft_cache%EFTOmegaP = self%EFTOmega%first_derivative( x, index=ind, coeff=mu )
953  eft_cache%EFTOmegaPP = self%EFTOmega%second_derivative( x, index=ind, coeff=mu )
954  eft_cache%EFTOmegaPPP = self%EFTOmega%third_derivative( x, index=ind, coeff=mu )
955  eft_cache%EFTc = 0._dl
956  eft_cache%EFTLambda = self%EFTLambda%value( x, index=ind, coeff=mu )
957  eft_cache%EFTcdot = 0._dl
958  eft_cache%EFTLambdadot = self%EFTLambda%first_derivative( x, index=ind, coeff=mu )
959 
960  end subroutine eftcambdesignerfrbackgroundeftfunctions
961 
962  ! ---------------------------------------------------------------------------------------------
964  subroutine eftcambdesignerfrsecondordereftfunctions( self, a, eft_par_cache, eft_cache )
965 
966  implicit none
967 
968  class(eftcamb_fr_designer) :: self
969  real(dl), intent(in) :: a
970  type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
971  type(eftcamb_timestep_cache ), intent(inout) :: eft_cache
972 
973  eft_cache%EFTGamma1V = 0._dl
974  eft_cache%EFTGamma1P = 0._dl
975  eft_cache%EFTGamma2V = 0._dl
976  eft_cache%EFTGamma2P = 0._dl
977  eft_cache%EFTGamma3V = 0._dl
978  eft_cache%EFTGamma3P = 0._dl
979  eft_cache%EFTGamma4V = 0._dl
980  eft_cache%EFTGamma4P = 0._dl
981  eft_cache%EFTGamma4PP = 0._dl
982  eft_cache%EFTGamma5V = 0._dl
983  eft_cache%EFTGamma5P = 0._dl
984  eft_cache%EFTGamma6V = 0._dl
985  eft_cache%EFTGamma6P = 0._dl
986 
987  end subroutine eftcambdesignerfrsecondordereftfunctions
988 
989  ! ---------------------------------------------------------------------------------------------
992  function eftcambdesignerfrcomputedtauda( self, a, eft_par_cache, eft_cache )
993 
994  implicit none
995 
996  class(eftcamb_fr_designer) :: self
997  real(dl), intent(in) :: a
998  type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
999  type(eftcamb_timestep_cache ), intent(inout) :: eft_cache
1000 
1001  real(dl) :: eftcambdesignerfrcomputedtauda
1002 
1003  real(dl) :: temp
1004 
1005  temp = eft_cache%grhoa2 +eft_par_cache%grhov*a*a*self%DesfRwDE%integral(a)
1006  eftcambdesignerfrcomputedtauda = sqrt(3/temp)
1007 
1008  end function eftcambdesignerfrcomputedtauda
1009 
1010  ! ---------------------------------------------------------------------------------------------
1012  subroutine eftcambdesignerfrcomputeadotoa( self, a, eft_par_cache, eft_cache )
1013 
1014  implicit none
1015 
1016  class(eftcamb_fr_designer) :: self
1017  real(dl), intent(in) :: a
1018  type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
1019  type(eftcamb_timestep_cache ), intent(inout) :: eft_cache
1020 
1021  eft_cache%grhov_t = eft_par_cache%grhov*self%DesfRwDE%integral(a)
1022  eft_cache%adotoa = sqrt( ( eft_cache%grhom_t +eft_cache%grhov_t )/3._dl )
1023 
1024  end subroutine eftcambdesignerfrcomputeadotoa
1025 
1026  ! ---------------------------------------------------------------------------------------------
1028  subroutine eftcambdesignerfrcomputehubbleder( self, a, eft_par_cache, eft_cache )
1029 
1030  implicit none
1031 
1032  class(eftcamb_fr_designer) :: self
1033  real(dl), intent(in) :: a
1034  type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
1035  type(eftcamb_timestep_cache ), intent(inout) :: eft_cache
1036 
1037  eft_cache%gpiv_t = self%DesfRwDE%value(a)*eft_cache%grhov_t
1038  eft_cache%Hdot = -0.5_dl*( eft_cache%adotoa**2 +eft_cache%gpresm_t +eft_cache%gpiv_t )
1039  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 ) &
1040  & +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) ) &
1041  & +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
1042 
1043  end subroutine eftcambdesignerfrcomputehubbleder
1044 
1045  ! ---------------------------------------------------------------------------------------------
1047  function eftcambdesignerfradditionalmodelstability( self, a, eft_par_cache, eft_cache )
1048 
1049  implicit none
1050 
1051  class(eftcamb_fr_designer) :: self
1052  real(dl), intent(in) :: a
1053  type(eftcamb_parameter_cache), intent(inout) :: eft_par_cache
1054  type(eftcamb_timestep_cache ), intent(inout) :: eft_cache
1055 
1056  logical :: eftcambdesignerfradditionalmodelstability
1057 
1058  eftcambdesignerfradditionalmodelstability = .true.
1059  if ( self%DesfRwDE%value(a) > -1._dl/3._dl ) eftcambdesignerfradditionalmodelstability = .false.
1060 
1061  end function eftcambdesignerfradditionalmodelstability
1062 
1063  ! ---------------------------------------------------------------------------------------------
1064 
1065 end module eftcamb_designer_fr
1066 
1067 !----------------------------------------------------------------------------------------
This module contains the definition of the constant parametrization, inheriting from parametrized_fun...
This module contains the definition of the Taylor expansion parametrization, around a=0...
This module contains the class that can be used for 1D linearly interpolated functions on an equispac...
This module contains the definition of the EFTCAMB caches. These are used to store parameters that ca...
This module contains the definition of the turning point parametrization, inheriting from parametrize...
logical, parameter debugeftcamb
EFTCAMB debug flag.This will turn on printing of many things to aid debugging the code...
Definition: 01_EFT_def.f90:48
This module contains the abstract definition of all the places where EFTCAMB interacts with CAMB in c...
This module contains the definition of the CPL parametrization, inheriting from parametrized_function...
real(dl) function, public zbrent(func, x1, x2, tol, funcZero, succes)
Brent root finding algorithm: This is used to solve numerically the equation func=funcZero by means o...
This module contains the definition of neutral parametrizations that can be used when parametrized fu...
This module contains the definitions of all the EFTCAMB root finding algorithms.
This module contains the definitions of all the EFTCAMB compile time flags.
Definition: 01_EFT_def.f90:25
This module contains the definition of the generalized Jassal-Bagla-Padmanabhan (JBP) parametrization...
This module contains the abstract class for generic parametrizations for 1D functions that are used b...
This module contains the relevant code for designer f(R) models.
subroutine, public zbrac(func, x1, x2, succes, funcZero)
Bracketing subroutine: This subroutine does a outward search for the smallest intervall containing a ...