tidy errors etc

camb/dark_energy.py

@@ -105,17 +105,18 @@ class AxionEffectiveFluid(DarkEnergyModel):
     Not well tested, but should serve to demonstrate how to make your own custom classes.
     """
     _fields_ = [
-        ("w_n", c_double),
-        ("om", c_double),
-        ("a_c", c_double),
-        ("theta_i", c_double)]
+        ("w_n", c_double, "effective equation of state parameter"),
+        ("fde_zc", c_double, "energy density fraction at z=zc"),
+        ("zc", c_double, "decay transition redshift (not same as peak of energy density fraction)"),
+        ("theta_i", c_double, "initial condition field value")]
+
     _fortran_class_name_ = 'TAxionEffectiveFluid'
     _fortran_class_module_ = 'DarkEnergyFluid'
 
-    def set_params(self, w_n, om, a_c, theta_i=None):
+    def set_params(self, w_n, fde_zc, zc, theta_i=None):
         self.w_n = w_n
-        self.om = om
-        self.a_c = a_c
+        self.fde_zc = fde_zc
+        self.zc = zc
         if theta_i is not None:
             self.theta_i = theta_i
 

camb/tests/camb_test.py

@@ -92,7 +92,7 @@ def testAssigments(self):
                                   'cosmomc_theta', 'YHe', 'wa', 'cs2', 'H0', 'mnu', 'Alens', 'TCMB', 'ns',
                                   'nrun', 'As', 'nt', 'r', 'w', 'omch2'})
         params2 = camb.get_valid_numerical_params(dark_energy_model='AxionEffectiveFluid')
-        self.assertEqual(params2.difference(params), {'om', 'w_n', 'a_c', 'theta_i'})
+        self.assertEqual(params2.difference(params), {'fde_zc', 'w_n', 'zc', 'theta_i'})
 
     def testBackground(self):
         pars = camb.CAMBparams()
@@ -704,5 +704,5 @@ def test_memory(self):
 
     def test_quintessence(self):
         pars = camb.set_params( ombh2=0.022, omch2=0.122,
-                               dark_energy_model='Quintessence', H0=67)
+                               dark_energy_model='EarlyQuintessence', H0=67)
         camb.get_background(pars)

fortran/DarkEnergyFluid.f90

@@ -23,10 +23,12 @@ module DarkEnergyFluid
     !This is an example, it's not supposed to be a rigorous model!  (not very well tested)
     type, extends(TDarkEnergyModel) :: TAxionEffectiveFluid
         real(dl) :: w_n = 1._dl !Effective equation of state when oscillating
-        real(dl) :: Om = 0._dl !Omega of the early DE component today (assumed to be negligible compared to omega_lambda)
-        real(dl) :: a_c  !transition scale factor
+        real(dl) :: fde_zc = 0._dl ! energy density fraction at a_c (not the same as peak dark energy fraction)
+        real(dl) :: zc  !transition redshift (scale factor a_c)
         real(dl) :: theta_i = const_pi/2 !Initial value
-        real(dl), private :: pow, omL, acpow, freq, n !cached internally
+        !om is Omega of the early DE component today (assumed to be negligible compared to omega_lambda)
+        !omL is the lambda component of the total dark energy omega
+        real(dl), private :: a_c, pow, om, omL, acpow, freq, n !cached internally
     contains
     procedure :: ReadParams =>  TAxionEffectiveFluid_ReadParams
     procedure, nopass :: PythonClass => TAxionEffectiveFluid_PythonClass
@@ -151,9 +153,12 @@ subroutine TAxionEffectiveFluid_ReadParams(this, Ini)
     class(TIniFile), intent(in) :: Ini
 
     call this%TDarkEnergyModel%ReadParams(Ini)
+    if (Ini%HasKey('AxionEffectiveFluid_a_c')) then
+        error stop 'AxionEffectiveFluid inputs changed to AxionEffectiveFluid_fde_zc and AxionEffectiveFluid_zc'
+    end if
     this%w_n  = Ini%Read_Double('AxionEffectiveFluid_w_n')
-    this%om  = Ini%Read_Double('AxionEffectiveFluid_om')
-    this%a_c  = Ini%Read_Double('AxionEffectiveFluid_a_c')
+    this%fde_zc  = Ini%Read_Double('AxionEffectiveFluid_fde_zc')
+    this%zc  = Ini%Read_Double('AxionEffectiveFluid_zc')
     call Ini%Read('AxionEffectiveFluid_theta_i', this%theta_i)
 
     end subroutine TAxionEffectiveFluid_ReadParams
@@ -184,10 +189,14 @@ subroutine TAxionEffectiveFluid_Init(this, State)
 
     select type(State)
     class is (CAMBdata)
-        this%is_cosmological_constant = this%om==0
+        this%is_cosmological_constant = this%fde_zc==0
         this%pow = 3*(1+this%w_n)
-        this%omL = State%Omega_de - this%om !Omega_de is total dark energy density today
+        this%a_c = 1/(1+this%zc)
         this%acpow = this%a_c**this%pow
+        !Omega in early de at z=0
+        this%om = 2*this%fde_zc/(1-this%fde_zc)*&
+            (State%grho_no_de(this%a_c)/this%a_c**4/State%grhocrit + State%Omega_de)/(1 + 1/this%acpow)
+        this%omL = State%Omega_de - this%om !Omega_de is total dark energy density today
         this%num_perturb_equations = 2
         if (this%w_n < 0.9999) then
             ! n <> infinity

fortran/DarkEnergyQuintessence.f90

@@ -267,7 +267,7 @@ subroutine TQuintessence_PerturbationEvolve(this, ayprime, w, w_ix, &
     ayprime(w_ix)= vq
     ayprime(w_ix+1) = - 2*adotoa*vq - k*z*phidot - k**2*clxq - a**2*clxq*this%Vofphi(phi,2)
 
-        end subroutine TQuintessence_PerturbationEvolve
+    end subroutine TQuintessence_PerturbationEvolve
 
     ! Early Quintessence example, axion potential from e.g. arXiv: 1908.06995
 
@@ -297,7 +297,6 @@ end function TEarlyQuintessence_VofPhi
 
 
     subroutine TEarlyQuintessence_Init(this, State)
-    use FileUtils
     use Powell
     class(TEarlyQuintessence), intent(inout) :: this
     class(TCAMBdata), intent(in), target :: State
@@ -308,7 +307,6 @@ subroutine TEarlyQuintessence_Init(this, State)
     real(dl), parameter :: splZero = 0._dl
     real(dl) lastsign, da_osc, last_a, a_c
     real(dl) initial_phi, initial_phidot, a2
-    Type(TTextFile) Fout
     real(dl), dimension(:), allocatable :: sampled_a, phi_a, phidot_a, fde
     integer npoints, tot_points, max_ix
     logical has_peak
@@ -348,30 +346,32 @@ subroutine TEarlyQuintessence_Init(this, State)
             end do
             call Timer%WriteTime('Timing for fitting')
         end if
-        if (.True.) then
-            if (this%DebugLevel>0) call Timer%Start()
-            !Minimize in log f, log m
-            param_min(1) = log(0.001_dl)
-            param_min(2) = log(1d-57)
-            param_max(1) = log(1e5_dl)
-            param_max(2) = log(1d-50)
-            if (Minimize%BOBYQA(this, match_fde_zc, 2, 5, log_params,param_min, &
-                param_max, 0.8_dl,1e-4_dl,this%DebugLevel,2000)) then
-                if (Minimize%Last_bestfit > 1e-3) then
-                    print *, 'ERROR finding good solution for fde, zc - best fit = ', Minimize%Last_bestfit
-                end if
-                this%f = exp(log_params(1))
-                this%m = exp(log_params(2))
-                if (this%DebugLevel>0) then
-                    call this%calc_zc_fde(fzero, xzero)
-                    print *, 'matched outputs Bobyqa zc, fde = ', fzero, xzero
-                end if
-            else
-                print *, 'ERROR finding solution for fde, zc'
+        if (this%DebugLevel>0) call Timer%Start()
+        !Minimize in log f, log m
+        param_min(1) = log(0.001_dl)
+        param_min(2) = log(1d-57)
+        param_max(1) = log(1e5_dl)
+        param_max(2) = log(1d-50)
+        if (Minimize%BOBYQA(this, match_fde_zc, 2, 5, log_params,param_min, &
+            param_max, 0.8_dl,1e-4_dl,this%DebugLevel,2000)) then
+            if (Minimize%Last_bestfit > 1e-3) then
+                global_error_flag = error_darkenergy
+                global_error_message= 'TEarlyQuintessence ERROR converging solution for fde, zc'
+                print *, 'last-bestfit= ', Minimize%Last_bestfit 
+                return
+            end if
+            this%f = exp(log_params(1))
+            this%m = exp(log_params(2))
+            if (this%DebugLevel>0) then
+                call this%calc_zc_fde(fzero, xzero)
+                print *, 'matched outputs Bobyqa zc, fde = ', fzero, xzero
             end if
-            if (this%DebugLevel>0) call Timer%WriteTime('Timing for BOBYQA')
+        else
+            global_error_flag = error_darkenergy
+            global_error_message= 'TEarlyQuintessence ERROR finding solution for fde, zc'
+            return
         end if
-
+        if (this%DebugLevel>0) call Timer%WriteTime('Timing for BOBYQA')
     end if
 
     this%dloga = (-this%log_astart)/(this%npoints-1)
@@ -449,7 +449,6 @@ subroutine TEarlyQuintessence_Init(this, State)
 
     ind=1
     afrom=this%log_astart
-    !call Fout%CreateOpenFile('z:\test.txt')
     do i=1, npoints-1
         aend = this%log_astart + this%dloga*i
         ix = i+1
@@ -468,15 +467,12 @@ subroutine TEarlyQuintessence_Init(this, State)
             lastsign= y(2)
         end if
 
-        !Define fde is ratio of energy density to total assuming the neutrinos fully relativistic
-        !adotrad = sqrt((this%grhog+this%grhornomass+sum(this%grhormass(1:this%CP%Nu_mass_eigenstates)))/3)
-
+        !Define fde as ratio of early dark energy density to total
         fde(ix) = 1/((this%state%grho_no_de(sampled_a(ix)) +  this%frac_lambda0*this%State%grhov*a2**2) &
             /(a2*(0.5d0* phidot_a(ix)**2 + a2*this%Vofphi(y(1),0))) + 1)
         if (max_ix==0 .and. ix > 2 .and. fde(ix)< fde(ix-1)) then
             max_ix = ix-1
         end if
-        !call Fout%Write(exp(aend), phi_a(ix), phidot_a(ix), fde(ix))
         if (sampled_a(ix)*(exp(this%dloga)-1)*this%min_steps_per_osc > da_osc) then
             !Step size getting too big to sample oscillations well
             exit
@@ -518,11 +514,8 @@ subroutine TEarlyQuintessence_Init(this, State)
         if (max_ix==0 .and. this%fde(ix)< this%fde(ix-1)) then
             max_ix = ix-1
         end if
-        !call Fout%Write(aend, this%phi_a(ix), this%phidot_a(ix), this%fde(ix))
     end do
 
-    !call Fout%Close()
-
     call spline(this%sampled_a,this%phi_a,tot_points,splZero,splZero,this%ddphi_a)
     call spline(this%sampled_a,this%phidot_a,tot_points,splZero,splZero,this%ddphidot_a)
     call spline(this%sampled_a,this%fde,tot_points,splZero,splZero,this%ddfde)
@@ -540,10 +533,12 @@ subroutine TEarlyQuintessence_Init(this, State)
         this%zc = 1/a_c-1
         this%fde_zc = this%fdeAta(a_c)
     else
-        print *, 'NO PEAK FINAL'
+        if (this%DebugLevel>0) write(*,*) 'TEarlyQuintessence: NO PEAK '
         this%zc = -1
     end if
-    print *, 'zc, fde used', this%zc, this%fde_zc
+    if (this%DebugLevel>0) then
+        write(*,*) 'TEarlyQuintessence zc, fde used', this%zc, this%fde_zc
+    end if
 
     end subroutine TEarlyQuintessence_Init
 
@@ -610,8 +605,10 @@ function match_fde_zc(this, x)
     this%f = exp(x(1))
     this%m = exp(x(2))
     call this%calc_zc_fde(zc, fde_zc)
-    match_fde_zc = (log(this%fde_zc)-log(fde_zc))**2 + (zc/this%zc-1)**2
-
+
+    match_fde_zc = (log(this%fde_zc)-log(fde_zc))**2 + (log(zc)-log(this%zc))**2
+    !print *, 'f, m, zc, fde_zc', this%f, this%m, zc, fde_zc, match_fde_zc
+
     end function match_fde_zc
 
     subroutine calc_zc_fde(this, z_c, fde_zc)

fortran/PowellConstrainedMinimize.f90

@@ -2,7 +2,7 @@
     !F90 translated version of plato.asu.edu/ftp/other_software/bobyqa.zip
     !Powell 2009 Method http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2009_06.pdf
 
-    !AL Sept 2012: translated to F90, changed to input "funkk" class function argument
+    !AL Sept 2012/2020: translated to F90, changed to input "funkk" class function argument
     !to minimize rather than calfun subroutine as original
 
     MODULE Powell
@@ -41,7 +41,7 @@ function BOBYQA (this, obj, funk, N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,MAXFUN)
     use iso_c_binding
     class(TBOBYQA) this
     class(*), target :: obj
-    real(dp), external :: funk !a class function f(obj,x)
+    real(dp), external :: funk !a class function f(obj,x) where obj is any class instance
     integer, intent(in)::n, npt, maxfun, iPrint
     real(dp), intent(in) :: RHOBEG, RHOEND
     logical BOBYQA

fortran/inidriver.f90

@@ -13,6 +13,7 @@ program driver
     if (InputFile == '') error stop 'No parameter input file'
 
     call CAMB_CommandLineRun(InputFile)
+    deallocate(InputFile) ! Just so no memory leaks in valgrind
 
     end program driver