EFTCAMB  Reference documentation for version 3.0
11_EFTCAMB_IC.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 
25 submodule(gaugeinterface) eftcamb_ic
26 
27 use precision
28 
29 implicit none
30 
31 contains
32 
33  ! ---------------------------------------------------------------------------------------------
35  module subroutine eftcambinitialconditions( y, ev, tau )
36 
37  use thermodata
38 
39  implicit none
40 
41  type(evolutionvars) ev
42  real(dl) :: y(ev%nvar)
43  real(dl) :: tau
44 
45  real(dl) :: a, a2, k, k2, etak, grhob_t, grhoc_t, grhor_t, grhog_t, grhov_t, gpres, grho_matter
46  real(dl) :: dgrho_matter, dgq, dgpnu, grho, adotoa, eft_grhonu, eft_gpinu, eft_grhonudot
47  real(dl) :: eft_gpinudot, grhormass_t, adotdota, cothxor, cs2, dopacity, dgrho, z, dz
48  real(dl) :: clxr, qr, pir, clxg, qg, pig, opacity, vb, clxc, clxb, eftpiedot
49  real(dl) :: wnu_arr(max_nu)
50 
51  integer :: nu_i
52 
53  k = ev%k_buf
54  k2 = ev%k2_buf
55  ! copy variables:
56  a = y(1) ! scale factor
57  etak = y(2) ! conformal time
58  clxc = y(3) ! cdm density perturbations
59  clxb = y(4) ! baryons density perturbations
60  vb = y(5) ! baryons velocity perturbations
61  ! process:
62  a2 = a*a
63  ! compute background densities of different species
64  grhob_t = grhob/a ! 8\pi G_N \rho_b a^2: baryons background density
65  grhoc_t = grhoc/a ! 8\pi G_N \rho_{cdm} a^2: cold dark matter background density
66  grhor_t = grhornomass/a2 ! 8\pi G_N \rho_{\nu} a^2: massless neutrinos background density
67  grhog_t = grhog/a2 ! 8\pi G_N \rho_{\gamma} a^2: radiation background density
68  grhov_t = 0._dl
69 
70  ! start computing background total pressure and total density:
71  gpres = 0._dl
72  grho_matter = grhob_t +grhoc_t
73  ! start computing total density and velocity perturbations:
74  dgrho_matter = grhob_t*clxb +grhoc_t*clxc ! 8\pi G_N \sum(rho_i*clx_i) a^2 : total density perturbation
75  dgq = grhob_t*vb ! 8\pi G_N \sum(rho_i+p_i)*v_i a^2 : total velocity perturbation
76  ! add massive neutrinos to both background and perturbations:
77 
78  dgpnu = 0._dl
79  if (cp%Num_Nu_Massive > 0) then
80  call massivenuvars(ev,y,a,grho_matter,gpres,dgrho_matter,dgq, wnu_arr, dgp=dgpnu)
81  end if
82 
83  ! add radiation, massless neutrinos and Lambda to total background density:
84  grho = grho_matter +grhor_t +grhog_t +grhov_t
85 
86  ! compute gpres: add radiation and massless neutrinos to massive neutrinos
87  gpres = gpres + (grhog_t+grhor_t)/3._dl
88  ! initialize the cache:
89  call ev%eft_cache%initialize()
90  ! start to fill the cache:
91  ev%eft_cache%a = a
92  ev%eft_cache%tau = tau
93  ev%eft_cache%k = k
94  ev%eft_cache%grhom_t = grho
95  ev%eft_cache%gpresm_t = gpres
96  ev%eft_cache%grhob_t = grhob_t
97  ev%eft_cache%grhoc_t = grhoc_t
98  ev%eft_cache%grhor_t = grhor_t
99  ev%eft_cache%grhog_t = grhog_t
100  ! compute the other things:
101  select type ( model => cp%EFTCAMB%model )
102  ! compute the background and the background EFT functions.
103  class is ( eftcamb_full_model )
104  ! background for full models. Here the expansion history is computed from the
105  ! EFT functions. Hence compute them first and then compute the expansion history.
106  call cp%EFTCAMB%model%compute_background_EFT_functions( a, cp%eft_par_cache , ev%eft_cache )
107  call cp%EFTCAMB%model%compute_adotoa( a, cp%eft_par_cache , ev%eft_cache )
108  class is ( eftcamb_designer_model )
109  ! background for designer models. Here the expansion history is parametrized
110  ! and does not depend on the EFT functions. Hence compute first the expansion history
111  ! and then the EFT functions.
112  call cp%EFTCAMB%model%compute_adotoa( a, cp%eft_par_cache , ev%eft_cache )
113  end select
114  ! store adotoa:
115  adotoa = ev%eft_cache%adotoa
116  ! compute massive neutrinos stuff:
117  ! Massive neutrinos mod:
118  if ( cp%Num_Nu_Massive /= 0 ) then
119  do nu_i = 1, cp%Nu_mass_eigenstates
120  eft_grhonu = 0._dl
121  eft_gpinu = 0._dl
122  eft_grhonudot = 0._dl
123  eft_gpinudot = 0._dl
124  grhormass_t=grhormass(nu_i)/a**2
125  call nu_background(a*nu_masses(nu_i),eft_grhonu,eft_gpinu)
126  ev%eft_cache%grhonu_tot = ev%eft_cache%grhonu_tot + grhormass_t*eft_grhonu
127  ev%eft_cache%gpinu_tot = ev%eft_cache%gpinu_tot + grhormass_t*eft_gpinu
128  ev%eft_cache%grhonudot_tot = ev%eft_cache%grhonudot_tot + grhormass_t*(nu_drho(a*nu_masses(nu_i) ,adotoa, eft_grhonu)&
129  & -4._dl*adotoa*eft_grhonu)
130  ev%eft_cache%gpinudot_tot = ev%eft_cache%gpinudot_tot + grhormass_t*(nu_pidot(a*nu_masses(nu_i),adotoa, eft_gpinu )&
131  & -4._dl*adotoa*eft_gpinu)
132  end do
133  end if
134  ! compute pressure dot:
135  ev%eft_cache%gpresdotm_t = -4._dl*adotoa*( grhog_t+grhor_t )/3._dl +ev%eft_cache%gpinudot_tot
136  ! compute remaining quantities related to H:
137  call cp%EFTCAMB%model%compute_H_derivs( a, cp%eft_par_cache , ev%eft_cache )
138  ! store:
139  adotdota = ev%eft_cache%Hdot +ev%eft_cache%adotoa**2
140  ! compute backgrond EFT functions if model is designer:
141  select type ( model => cp%EFTCAMB%model )
142  class is ( eftcamb_designer_model )
143  call cp%EFTCAMB%model%compute_background_EFT_functions( a, cp%eft_par_cache , ev%eft_cache )
144  end select
145  ! compute all other background stuff:
146  call cp%EFTCAMB%model%compute_rhoQPQ( a, cp%eft_par_cache , ev%eft_cache )
147  ! compute second order EFT functions:
148  call cp%EFTCAMB%model%compute_secondorder_EFT_functions( a, cp%eft_par_cache , ev%eft_cache )
149  !
150  cothxor=1._dl/tau
151 
152  ! Get sound speed and ionisation fraction.
153  if ( ev%TightCoupling ) then
154  call thermo( tau, cs2, opacity, dopacity )
155  else
156  call thermo( tau, cs2, opacity )
157  end if
158 
159  ! define total density perturbations:
160  dgrho = dgrho_matter
161 
162  ! compute z and zdot in GR:
163  z = (0.5_dl*dgrho/k + etak)/adotoa
164  dz = -2._dl*adotoa*z + etak
165 
166  if ( ev%no_nu_multpoles ) then
167  !RSA approximation of arXiv:1104.2933, dropping opactity terms in the velocity
168  !Approximate total density variables with just matter terms
169  clxr = -4*dz/k
170  qr = -4._dl/3*z
171  else
172  ! Massless neutrinos
173  clxr = y(ev%r_ix)
174  qr = y(ev%r_ix+1)
175  endif
176 
177  if ( ev%no_phot_multpoles ) then
178  if ( .not. ev%no_nu_multpoles ) then
179  clxg = -4*dz/k -4/k*opacity*( vb +z )
180  qg = -4._dl/3*z
181  else
182  clxg = clxr -4/k*opacity*( vb +z )
183  qg = qr
184  end if
185  else
186  ! Photons
187  clxg = y(ev%g_ix)
188  qg = y(ev%g_ix+1)
189  end if
190 
191  ! add radiation to total density and velocity perturbations:
192  dgrho = dgrho +grhog_t*clxg +grhor_t*clxr ! 8\pi G_N \sum(rho_i*clx_i) a^2 : total density perturbation
193  dgq = dgq +grhog_t*qg +grhor_t*qr ! 8\pi G_N \sum(rho_i+p_i)*v_i a^2 : total velocity perturbation
194 
195  ! recompute z and zdot in GR:
196  z = (0.5_dl*dgrho/k + etak)/adotoa
197  dz = -2._dl*adotoa*z + etak
198 
199  ! complete the cache:
200  ev%eft_cache%z = z
201  ev%eft_cache%dz = dz
202  ev%eft_cache%clxg = clxg
203  ev%eft_cache%clxr = clxr
204  ev%eft_cache%dgpnu = dgpnu
205  ev%eft_cache%dgrho = dgrho
206  ev%eft_cache%dgq = dgq
207 
208  ! compute pi field equation factors:
209  call cp%EFTCAMB%model%compute_pi_factors( a, cp%eft_par_cache , ev%eft_cache )
210 
211  if ( ev%eft_cache%EFTGamma1V /= 0._dl .or. ev%eft_cache%EFTGamma2V /= 0._dl .or. &
212  &ev%eft_cache%EFTGamma3V /= 0._dl .or. ev%eft_cache%EFTGamma4V /= 0._dl .or. &
213  &ev%eft_cache%EFTGamma5V /= 0._dl .or. ev%eft_cache%EFTGamma6V /= 0._dl ) then
214 
215  eftpiedot = 0._dl
216 
217  else
218 
219  eftpiedot = ( ev%eft_cache%EFTcdot +0.75_dl*(2._dl*ev%eft_cache%adotoa*ev%eft_cache%Hdot*a*a*ev%eft_cache%EFTOmegaP**2/(1._dl+ev%eft_cache%EFTOmegaV)&
220  & +2._dl*ev%eft_cache%adotoa**3*a**3*ev%eft_cache%EFTOmegaP*ev%eft_cache%EFTOmegaPP/(1._dl + ev%eft_cache%EFTOmegaV)-a**3*ev%eft_cache%adotoa**3*ev%eft_cache%EFTOmegaP**3/(1._dl+ev%eft_cache%EFTOmegaV)**2))*k*ev%eft_cache%z &
221  & +(ev%eft_cache%EFTc+0.75_dl*ev%eft_cache%adotoa**2*(a*a*ev%eft_cache%EFTOmegaP**2)/(1._dl+ev%eft_cache%EFTOmegaV))*k*ev%eft_cache%dz +2._dl*ev%eft_cache%adotoa*ev%eft_cache%EFTpiE&
222  & +(-ev%eft_cache%dgrho +(ev%eft_cache%grhog_t*ev%eft_cache%clxg+ev%eft_cache%grhor_t*ev%eft_cache%clxr))/4._dl*(+a*ev%eft_cache%adotoa**2*ev%eft_cache%EFTOmegaP +a*ev%eft_cache%Hdot*ev%eft_cache%EFTOmegaP+a*a*ev%eft_cache%adotoa**2*ev%eft_cache%EFTOmegaPP &
223  & -a*a*ev%eft_cache%adotoa**2*ev%eft_cache%EFTOmegaP**2/(1._dl+ev%eft_cache%EFTOmegaV))/(1._dl+ev%eft_cache%EFTOmegaV)&
224  & -ev%eft_cache%adotoa/4._dl*a*ev%eft_cache%EFTOmegaP/(1._dl+ev%eft_cache%EFTOmegaV)*(ev%eft_cache%grhob_t*(-k*(ev%eft_cache%z +vb)-3._dl*ev%eft_cache%adotoa*clxb) + ev%eft_cache%grhoc_t*(-k*ev%eft_cache%z -3._dl*ev%eft_cache%adotoa*clxc))
225 
226  end if
227 
228  if ( ev%eft_cache%EFTpiC +k2*ev%eft_cache%EFTpiD1 +k2**2*ev%eft_cache%EFTpiD2 /= 0._dl ) then
229  y(ev%w_ix) = -cp%eft_par_cache%h0_Mpc*ev%eft_cache%EFTpiE/( ev%eft_cache%EFTpiC +k2*ev%eft_cache%EFTpiD1 +k2**2*ev%eft_cache%EFTpiD2 )
230  y(ev%w_ix+1) = cp%eft_par_cache%h0_Mpc*(ev%eft_cache%EFTpiE/( ev%eft_cache%EFTpiC +k2*ev%eft_cache%EFTpiD1 +k2**2*ev%eft_cache%EFTpiD2 )**2*a*ev%eft_cache%adotoa*(eftpicdotfunction(a,k)+k*k*eftpiddotfunction(a,k))&
231  & -eftpiedot/(ev%eft_cache%EFTpiC +k2*ev%eft_cache%EFTpiD1 +k2**2*ev%eft_cache%EFTpiD2 ) )
232  else
233  y(ev%w_ix) = 0._dl
234  y(ev%w_ix+1) = 0._dl
235  end if
236 
237  ! reset the cache:
238  call ev%eft_cache%initialize()
239 
240  end subroutine eftcambinitialconditions
241 
242  !----------------------------------------------------------------------------------------
244  function eftpicfunction( a, k )
245 
246  implicit none
247 
248  real(dl), intent(in) :: a
249  real(dl), intent(in) :: k
250  real(dl) :: eftpicfunction
251 
252  type(eftcamb_timestep_cache) :: eft_cache
253  integer :: nu_i
254  real(dl) :: grhormass_t, eft_gpinudot, eft_grhonudot, eft_gpinu, eft_grhonu
255  real(dl) :: adotoa, adotdota, grho, grhob_t, grhoc_t, grhor_t, grhog_t, gpres, grho_matter, a2
256 
257  ! initialize the cache:
258  call eft_cache%initialize()
259 
260  ! prepare:
261  a2 = a*a
262 
263  ! compute background densities of different species
264  grhob_t = grhob/a ! 8\pi G_N \rho_b a^2: baryons background density
265  grhoc_t = grhoc/a ! 8\pi G_N \rho_{cdm} a^2: cold dark matter background density
266  grhor_t = grhornomass/a2 ! 8\pi G_N \rho_{\nu} a^2: massless neutrinos background density
267  grhog_t = grhog/a2 ! 8\pi G_N \rho_{\gamma} a^2: radiation background density
268 
269  ! start computing background total pressure and total density:
270  grho_matter = grhob_t +grhoc_t
271  gpres = (grhog_t+grhor_t)/3._dl
272 
273  ! add radiation, massless neutrinos and Lambda to total background density:
274  grho = grho_matter +grhor_t +grhog_t
275 
276  if ( cp%Num_Nu_Massive /= 0 ) then
277  do nu_i = 1, cp%Nu_mass_eigenstates
278  eft_grhonu = 0._dl
279  eft_gpinu = 0._dl
280  grhormass_t=grhormass(nu_i)/a**2
281  call nu_background(a*nu_masses(nu_i),eft_grhonu,eft_gpinu)
282  eft_cache%grhonu_tot = eft_cache%grhonu_tot + grhormass_t*eft_grhonu
283  eft_cache%gpinu_tot = eft_cache%gpinu_tot + grhormass_t*eft_gpinu
284  end do
285  end if
286 
287  grho = grho +eft_cache%grhonu_tot
288  gpres = gpres +eft_cache%gpinu_tot
289 
290  ! start to fill the cache:
291  eft_cache%a = a
292  eft_cache%k = k
293  eft_cache%grhom_t = grho
294  eft_cache%gpresm_t = gpres
295  eft_cache%grhob_t = grhob_t
296  eft_cache%grhoc_t = grhoc_t
297  eft_cache%grhor_t = grhor_t
298  eft_cache%grhog_t = grhog_t
299  ! compute the other things:
300  select type ( model => cp%EFTCAMB%model )
301  ! compute the background and the background EFT functions.
302  class is ( eftcamb_full_model )
303  ! background for full models. Here the expansion history is computed from the
304  ! EFT functions. Hence compute them first and then compute the expansion history.
305  call cp%EFTCAMB%model%compute_background_EFT_functions( a, cp%eft_par_cache , eft_cache )
306  call cp%EFTCAMB%model%compute_adotoa( a, cp%eft_par_cache , eft_cache )
307  class is ( eftcamb_designer_model )
308  ! background for designer models. Here the expansion history is parametrized
309  ! and does not depend on the EFT functions. Hence compute first the expansion history
310  ! and then the EFT functions.
311  call cp%EFTCAMB%model%compute_adotoa( a, cp%eft_par_cache , eft_cache )
312  end select
313  ! store adotoa:
314  adotoa = eft_cache%adotoa
315  ! compute massive neutrinos stuff:
316  ! Massive neutrinos mod:
317  eft_cache%grhonu_tot = 0._dl
318  eft_cache%gpinu_tot = 0._dl
319  if ( cp%Num_Nu_Massive /= 0 ) then
320  do nu_i = 1, cp%Nu_mass_eigenstates
321  eft_grhonu = 0._dl
322  eft_gpinu = 0._dl
323  eft_grhonudot = 0._dl
324  eft_gpinudot = 0._dl
325  grhormass_t=grhormass(nu_i)/a**2
326  call nu_background(a*nu_masses(nu_i),eft_grhonu,eft_gpinu)
327  eft_cache%grhonu_tot = eft_cache%grhonu_tot + grhormass_t*eft_grhonu
328  eft_cache%gpinu_tot = eft_cache%gpinu_tot + grhormass_t*eft_gpinu
329  eft_cache%grhonudot_tot = eft_cache%grhonudot_tot + grhormass_t*(nu_drho(a*nu_masses(nu_i) ,adotoa, eft_grhonu)&
330  & -4._dl*adotoa*eft_grhonu)
331  eft_cache%gpinudot_tot = eft_cache%gpinudot_tot + grhormass_t*(nu_pidot(a*nu_masses(nu_i),adotoa, eft_gpinu )&
332  & -4._dl*adotoa*eft_gpinu)
333  end do
334  end if
335  ! compute pressure dot:
336  eft_cache%gpresdotm_t = -4._dl*adotoa*( grhog_t+grhor_t )/3._dl +eft_cache%gpinudot_tot
337  ! compute remaining quantities related to H:
338  call cp%EFTCAMB%model%compute_H_derivs( a, cp%eft_par_cache , eft_cache )
339  ! store:
340  adotdota = eft_cache%Hdot +eft_cache%adotoa**2
341  ! compute backgrond EFT functions if model is designer:
342  select type ( model => cp%EFTCAMB%model )
343  class is ( eftcamb_designer_model )
344  call cp%EFTCAMB%model%compute_background_EFT_functions( a, cp%eft_par_cache , eft_cache )
345  end select
346  ! compute all other background stuff:
347  call cp%EFTCAMB%model%compute_rhoQPQ( a, cp%eft_par_cache , eft_cache )
348  ! compute second order EFT functions:
349  call cp%EFTCAMB%model%compute_secondorder_EFT_functions( a, cp%eft_par_cache , eft_cache )
350 
351  ! compute Einstein equations factors:
352  call cp%EFTCAMB%model%compute_pi_factors( a, cp%eft_par_cache , eft_cache )
353 
354  eftpicfunction = eft_cache%EFTpiC
355 
356  return
357 
358  end function eftpicfunction
359 
360  !----------------------------------------------------------------------------------------
362  function eftpidfunction(a,k)
363 
364  implicit none
365 
366  real(dl), intent(in) :: a
367  real(dl), intent(in) :: k
368  real(dl) :: eftpidfunction
369 
370  type(eftcamb_timestep_cache) :: eft_cache
371 
372  integer :: nu_i
373  real(dl) :: grhormass_t, eft_gpinudot, eft_grhonudot, eft_gpinu, eft_grhonu
374  real(dl) :: adotoa, adotdota, grho, grhob_t, grhoc_t, grhor_t, grhog_t, gpres, grho_matter, a2
375 
376  ! initialize the cache:
377  call eft_cache%initialize()
378 
379  ! prepare:
380  a2 = a*a
381 
382  ! compute background densities of different species
383  grhob_t = grhob/a ! 8\pi G_N \rho_b a^2: baryons background density
384  grhoc_t = grhoc/a ! 8\pi G_N \rho_{cdm} a^2: cold dark matter background density
385  grhor_t = grhornomass/a2 ! 8\pi G_N \rho_{\nu} a^2: massless neutrinos background density
386  grhog_t = grhog/a2 ! 8\pi G_N \rho_{\gamma} a^2: radiation background density
387 
388  ! start computing background total pressure and total density:
389  grho_matter = grhob_t +grhoc_t
390  gpres = (grhog_t+grhor_t)/3._dl
391 
392  ! add radiation, massless neutrinos and Lambda to total background density:
393  grho = grho_matter +grhor_t +grhog_t
394 
395  if ( cp%Num_Nu_Massive /= 0 ) then
396  do nu_i = 1, cp%Nu_mass_eigenstates
397  eft_grhonu = 0._dl
398  eft_gpinu = 0._dl
399  grhormass_t=grhormass(nu_i)/a**2
400  call nu_background(a*nu_masses(nu_i),eft_grhonu,eft_gpinu)
401  eft_cache%grhonu_tot = eft_cache%grhonu_tot + grhormass_t*eft_grhonu
402  eft_cache%gpinu_tot = eft_cache%gpinu_tot + grhormass_t*eft_gpinu
403  end do
404  end if
405 
406  grho = grho +eft_cache%grhonu_tot
407  gpres = gpres +eft_cache%gpinu_tot
408 
409  ! start to fill the cache:
410  eft_cache%a = a
411  eft_cache%k = k
412  eft_cache%grhom_t = grho
413  eft_cache%gpresm_t = gpres
414  eft_cache%grhob_t = grhob_t
415  eft_cache%grhoc_t = grhoc_t
416  eft_cache%grhor_t = grhor_t
417  eft_cache%grhog_t = grhog_t
418  ! compute the other things:
419  select type ( model => cp%EFTCAMB%model )
420  ! compute the background and the background EFT functions.
421  class is ( eftcamb_full_model )
422  ! background for full models. Here the expansion history is computed from the
423  ! EFT functions. Hence compute them first and then compute the expansion history.
424  call cp%EFTCAMB%model%compute_background_EFT_functions( a, cp%eft_par_cache , eft_cache )
425  call cp%EFTCAMB%model%compute_adotoa( a, cp%eft_par_cache , eft_cache )
426  class is ( eftcamb_designer_model )
427  ! background for designer models. Here the expansion history is parametrized
428  ! and does not depend on the EFT functions. Hence compute first the expansion history
429  ! and then the EFT functions.
430  call cp%EFTCAMB%model%compute_adotoa( a, cp%eft_par_cache , eft_cache )
431  end select
432  ! store adotoa:
433  adotoa = eft_cache%adotoa
434  ! compute massive neutrinos stuff:
435  ! Massive neutrinos mod:
436  eft_cache%grhonu_tot = 0._dl
437  eft_cache%gpinu_tot = 0._dl
438  if ( cp%Num_Nu_Massive /= 0 ) then
439  do nu_i = 1, cp%Nu_mass_eigenstates
440  eft_grhonu = 0._dl
441  eft_gpinu = 0._dl
442  eft_grhonudot = 0._dl
443  eft_gpinudot = 0._dl
444  grhormass_t=grhormass(nu_i)/a**2
445  call nu_background(a*nu_masses(nu_i),eft_grhonu,eft_gpinu)
446  eft_cache%grhonu_tot = eft_cache%grhonu_tot + grhormass_t*eft_grhonu
447  eft_cache%gpinu_tot = eft_cache%gpinu_tot + grhormass_t*eft_gpinu
448  eft_cache%grhonudot_tot = eft_cache%grhonudot_tot + grhormass_t*(nu_drho(a*nu_masses(nu_i) ,adotoa, eft_grhonu)&
449  & -4._dl*adotoa*eft_grhonu)
450  eft_cache%gpinudot_tot = eft_cache%gpinudot_tot + grhormass_t*(nu_pidot(a*nu_masses(nu_i),adotoa, eft_gpinu )&
451  & -4._dl*adotoa*eft_gpinu)
452  end do
453  end if
454  ! compute pressure dot:
455  eft_cache%gpresdotm_t = -4._dl*adotoa*( grhog_t+grhor_t )/3._dl +eft_cache%gpinudot_tot
456  ! compute remaining quantities related to H:
457  call cp%EFTCAMB%model%compute_H_derivs( a, cp%eft_par_cache , eft_cache )
458  ! store:
459  adotdota = eft_cache%Hdot +eft_cache%adotoa**2
460  ! compute backgrond EFT functions if model is designer:
461  select type ( model => cp%EFTCAMB%model )
462  class is ( eftcamb_designer_model )
463  call cp%EFTCAMB%model%compute_background_EFT_functions( a, cp%eft_par_cache , eft_cache )
464  end select
465  ! compute all other background stuff:
466  call cp%EFTCAMB%model%compute_rhoQPQ( a, cp%eft_par_cache , eft_cache )
467  ! compute second order EFT functions:
468  call cp%EFTCAMB%model%compute_secondorder_EFT_functions( a, cp%eft_par_cache , eft_cache )
469 
470  ! compute Einstein equations factors:
471  call cp%EFTCAMB%model%compute_pi_factors( a, cp%eft_par_cache , eft_cache )
472 
473  eftpidfunction = eft_cache%EFTpiD1 + k*k*eft_cache%EFTpiD1
474 
475  return
476 
477  end function eftpidfunction
478 
479  !----------------------------------------------------------------------------------------
482  function eftpicdotfunction(a,k)
483 
484  implicit none
485 
486  real(dl), intent(in) :: a
487  real(dl), intent(in) :: k
488  real(dl) :: eftpicdotfunction
489 
490  real(dl) err
491 
492  eftpicdotfunction = dfridr( eftpicfunction, a, k, 0.03_dl*a, err )
493 
494  return
495 
496  end function eftpicdotfunction
497 
498  !----------------------------------------------------------------------------------------
501  function eftpiddotfunction(a,k)
502 
503  implicit none
504 
505  real(dl), intent(in) :: a
506  real(dl), intent(in) :: k
507  real(dl) :: eftpiddotfunction
508 
509 
510  real(dl) err
511 
512  eftpiddotfunction = dfridr( eftpidfunction, a, k, 0.03_dl*a, err )
513 
514  return
515 
516  end function eftpiddotfunction
517 
518  !----------------------------------------------------------------------------------------
521  function dfridr( func, x, k, h, err )
522 
523  implicit none
524 
525  integer,parameter :: ntab = 100
526  real(dl) dfridr,err,h,x,k,func
527  external func
528 
529  real(dl), parameter :: con=1.4_dl ! decrease of the stepsize.
530  real(dl), parameter :: con2=con*con
531  real(dl), parameter :: big=1.d+30
532  real(dl), parameter :: safe=2._dl
533 
534  integer i,j
535  real(dl) errt, fac, hh
536  real(dl), dimension(ntab,ntab) :: a
537 
538  if (h.eq.0._dl) h = 1.d-8
539 
540  hh=h
541  a(1,1)=(func(x+hh,k)-func(x-hh,k))/(2.0_dl*hh)
542  err=big
543 
544  do 12 i=2,ntab
545  hh=hh/con
546  a(1,i)=(func(x+hh,k)-func(x-hh,k))/(2.0_dl*hh)
547  fac=con2
548  do 11 j=2,i
549  a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1._dl)
550  fac=con2*fac
551  errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
552  if (errt.le.err) then
553  err=errt
554  dfridr=a(j,i)
555  endif
556 11 continue
557  if(abs(a(i,i)-a(i-1,i-1)).ge.safe*err)return
558 12 continue
559  return
560  end function dfridr
561 
562  !----------------------------------------------------------------------------------------
563 
564 end submodule eftcamb_ic
565 
566 !----------------------------------------------------------------------------------------
Generated on Thu Jan 19 2017 12:48:10 for EFTCAMB by   doxygen 1.8.11