fix get_matter_power_spectrum for cross-spectra

Author: cmbant

camb/__init__.py

@@ -7,7 +7,7 @@
 __author__ = "Antony Lewis"
 __contact__ = "antony at cosmologist dot info"
 __url__ = "https://camb.readthedocs.io"
-__version__ = "1.3.6"
+__version__ = "1.3.7"
 
 from . import baseconfig
 

camb/cambdll.dll

@@ -446,8 +446,8 @@ def lensed_cl_derivatives(cls, clpp, lmax=None, theta_max=np.pi / 32,
     :param apodize_point_width: if theta_max is set, apodize around the cut using half Gaussian of approx
         width apodize_point_width/lmax*pi
     :param sampling_factor: npoints = int(sampling_factor*lmax)+1
-    :return: array dCL[ix, ell, L], where ix=0,1,2,3 are T, EE, BB, TE and r
-             esult is :math:`d[D^{\rm ix}_\ell]/ d (\log C^{\phi}_L)`
+    :return: array dCL[ix, ell, L], where ix=0,1,2,3 are T, EE, BB, TE and result
+               is :math:`d[D^{\rm ix}_\ell]/ d (\log C^{\phi}_L)`
 
     """
 

camb/correlations.py

@@ -365,8 +365,8 @@ def set_H0_for_theta(self, theta, cosmomc_approx=False, theta_H0_range=(10, 100)
         :param cosmomc_approx: if true, use approximate fitting formula for :math:`z_\star`,
                                if false do full numerical calculation
         :param theta_H0_range: min, max iterval to search for H0 (in km/s/Mpc)
-        :param est_H0: an initial guess for H0 in km/s/Mpc, used in the case comsomc_approx=False.
-        :param iteration_threshold: differnce in H0 from est_H0 for which to iterate, used for cosmomc_approx=False
+        :param est_H0: an initial guess for H0 in km/s/Mpc, used in the case cosmomc_approx=False.
+        :param iteration_threshold: difference in H0 from est_H0 for which to iterate, used for cosmomc_approx=False
         """
 
         if not (0.001 < theta < 0.1):
@@ -425,7 +425,7 @@ def set_cosmology(self, H0: Optional[float] = None, ombh2=0.022, omch2=0.12, omk
         (cosmomc_theta, which is based on a fitting forumula for simple models, or thetastar, which is numerically
         calculated more generally). Note that you must have already set the dark energy model, you can't use
         set_cosmology with theta and then change the background evolution (which would change theta at the calculated
-        H0 value).Likewise the dark energy model cannot depend explicitly on H0.
+        H0 value). Likewise the dark energy model cannot depend explicitly on H0.
 
         :param H0: Hubble parameter today in km/s/Mpc. Can leave unset and instead set thetastar or cosmomc_theta
                   (which solves for the required H0).

camb/model.py

@@ -3,7 +3,7 @@ module config
     use constants, only: const_twopi
     implicit none
 
-    character(LEN=*), parameter :: version = '1.3.6'
+    character(LEN=*), parameter :: version = '1.3.7'
 
     integer :: FeedbackLevel = 0 !if >0 print out useful information about the model
 

fortran/config.f90

@@ -737,9 +737,9 @@ end function CAMBdata_AngularDiameterDistance2
 
     subroutine CAMBdata_AngularDiameterDistance2Arr(this, arr, z1, z2, n)
     class(CAMBdata) :: this
+    integer, intent(in) :: n
     real(dl), intent(out) :: arr(n)
     real(dl), intent(in) :: z1(n), z2(n)
-    integer, intent(in) :: n
     integer i
 
     !$OMP PARALLEL DO DEFAULT(SHARED),SCHEDULE(STATIC)
@@ -3461,7 +3461,8 @@ subroutine Transfer_GetMatterPowerD(State, MTrans, outpower, itf_PK, minkh, dlnk
     real(dl) matpower(MTrans%num_q_trans), kh, kvals(MTrans%num_q_trans), ddmat(MTrans%num_q_trans)
     real(dl) atransfer,xi, a0, b0, ho, logmink,k, h
     integer itf
-    integer :: s1,s2
+    integer :: s1,s2, sign
+    logical log_interp
     real(dl), allocatable :: ratio(:)
 
     s1 = PresentDefault (transfer_power_var, var1)
@@ -3489,10 +3490,19 @@ subroutine Transfer_GetMatterPowerD(State, MTrans, outpower, itf_PK, minkh, dlnk
         atransfer=MTrans%TransferData(s1,ik,itf)*MTrans%TransferData(s2,ik,itf)
         if (State%CP%NonLinear/=NonLinear_none .and. State%CP%NonLinear/=NonLinear_Lens) &
             atransfer = atransfer* ratio(ik)**2 !only one element, this itf
-        matpower(ik) = log(atransfer*k*const_pi*const_twopi*h**3)
+        matpower(ik) = atransfer*k*const_pi*const_twopi*h**3
         !Put in power spectrum later: transfer functions should be smooth, initial power may not be
     end do
-
+    sign = 1
+    log_interp = .true.
+    if (any(matpower <= 0)) then
+        if (all(matpower < 0)) then
+            sign = -1
+        else
+            log_interp = .false.
+        endif
+    endif
+    if (log_interp) matpower = log(sign*matpower)
     call spline(kvals,matpower,MTrans%num_q_trans,cllo,clhi,ddmat)
 
     llo=1
@@ -3527,7 +3537,9 @@ subroutine Transfer_GetMatterPowerD(State, MTrans, outpower, itf_PK, minkh, dlnk
         lastix = lastix+1
     end do
 
-    outpower = exp(max(-30.d0,outpower))
+    if (log_interp) then
+        outpower = sign*exp(max(-30.d0,outpower))
+    end if
     associate(InitPower => State%CP%InitPower)
         do il = 1, npoints
             k = exp(logmink + dlnkh*(il-1))*h