The ParaMonte Documentation Website
Current view: top level - kernel - Cosmology_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: Serial Kernel - Code Coverage Report Lines: 40 40 100.0 %
Date: 2021-01-08 13:03:42 Functions: 6 6 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !> \brief This module contains procedures and constants for cosmological calculations.
      44             : !> \author Amir Shahmoradi
      45             : 
      46             : module Cosmology_mod
      47             : 
      48             :     use Constants_mod, only: RK, PI
      49             : 
      50             :     implicit none
      51             : 
      52             :     character(*), parameter :: MODULE_NAME = "@Cosmology_mod"
      53             : 
      54             :     ! Cosmological constants
      55             : 
      56             :     real(RK), parameter :: LIGHT_SPEED = 3.e5_RK                                    !< @public LIGHT_SPEED is the speed of light (Km/s).
      57             :     real(RK), parameter :: HUBBLE_CONST = 7.1e1_RK                                  !< @public HUBBLE_CONST is the Hubble constant in units of km/s/MPc.
      58             :     real(RK), parameter :: HUBBLE_TIME_GYRS = 13.8_RK                               !< @public hubble time (liddle 2003, page 57) in units of gyrs
      59             :     real(RK), parameter :: INVERSE_HUBBLE_CONST = 1._RK / HUBBLE_CONST              !< @public inverse of HUBBLE_CONST: 0.014084507042254
      60             :     real(RK), parameter :: LS2HC = LIGHT_SPEED / HUBBLE_CONST                       !< @public the speed of light in units of km/s divided by the Hubble constant.
      61             :     real(RK), parameter :: LOGLS2HC = log(LIGHT_SPEED) - log(HUBBLE_CONST)          !< @public log speed of light in units of km/s divided by the Hubble constant.
      62             :     real(RK), parameter :: MPC2CM = 3.09e24_RK                                      !< @public 1 Mega Parsec = MPC2CM centimeters.
      63             :     real(RK), parameter :: LOGMPC2CMSQ4PI   = log(4._RK*PI) + 2._RK*log(MPC2CM)     !< @public log(MegaParsec2centimeters).
      64             :     real(RK), parameter :: LOG10MPC2CMSQ4PI = log10(4._RK*PI) + 2._RK*log10(MPC2CM) !< @public log10(MegaParsec2centimeters).
      65             :     real(RK), parameter :: OMEGA_DE = 0.7_RK                                        !< @public Dark Energy density.
      66             :     real(RK), parameter :: OMEGA_DM = 0.3_RK                                        !< @public Dark Matter density.
      67             : 
      68             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      69             : 
      70             : contains
      71             : 
      72             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      73             : 
      74             :     !> \brief
      75             :     !> Return the natural logarithm of the differential comoving volume of cosmos.
      76             :     !>
      77             :     !> @param[in]   zplus1              : The redshift plus 1.
      78             :     !> @param[in]   logzplus1           : The logarithm of redshift plus 1.
      79             :     !> @param[in]   twiceLogLumDisMpc   : The term representing the logarithm of the luminosity distance squared.
      80             :     !>
      81             :     !> \return
      82             :     !> `logdvdz` : The natural logarithm of the differential comoving volume of cosmos.
      83           1 :     pure function getlogdvdz(zplus1,logzplus1,twiceLogLumDisMpc) result(logdvdz)
      84             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
      85             :         !DEC$ ATTRIBUTES DLLEXPORT :: getlogdvdz
      86             : #endif
      87             :         use Constants_mod, only: RK, PI
      88             :         implicit none
      89             :         real(RK), intent(in)    :: zplus1, logzplus1, twiceLogLumDisMpc
      90             :         real(RK), parameter     :: LOG_COEF = log(4*PI*LS2HC)
      91             :         real(RK)                :: logdvdz
      92           1 :         logdvdz = LOG_COEF + twiceLogLumDisMpc - ( 2*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) )
      93           2 :     end function getlogdvdz
      94             : 
      95             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      96             : 
      97             :     !> \brief
      98             :     !> Return the approximate logarithm of the cosmological luminosity distance in units of MPc.
      99             :     !>
     100             :     !> @param[in]   zplus1 : The redshift plus 1.
     101             :     !>
     102             :     !> \return
     103             :     !> `logLumDisWicMpc` : The approximate logarithm of the cosmological luminosity distance.
     104             :     !>
     105             :     !> \warning
     106             :     !> The approximation is accurate with a relative error of 0.001 for any redshift above 0.1, or `zplus1 >= 1.1`.
     107             :     !! Note that for redshifts less than 0.1, the error in the calculated luminosity distance grows to more than 0.001.
     108             :     !! This algorithm should therefore not be used for zplus1<0.1.
     109             :     !>
     110             :     !> \remark
     111             :     !> The distance is calculated according to approximate algorithm of Wickramasinghe & Okwatta (2010).
     112          16 :     pure function getLogLumDisWicMpc(zplus1) result(logLumDisWicMpc)
     113             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     114             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLogLumDisWicMpc
     115             : #endif
     116           1 :         use Constants_mod, only: RK
     117             :         implicit none
     118             :         real(RK), intent(in)    :: zplus1   ! redshift + 1
     119             :         real(RK)                :: logLumDisWicMpc
     120          16 :         real(RK)                :: alpha1, x1, x1Sq, psix1
     121             :         real(RK), parameter     :: TWICE_OMEGA_DE_OVER_OMEGA_DM = 2._RK * OMEGA_DE / OMEGA_DM
     122             :         real(RK), parameter     :: ALPHA0 = 1._RK + TWICE_OMEGA_DE_OVER_OMEGA_DM
     123             :         real(RK), parameter     :: X0 = log( ALPHA0 + sqrt(ALPHA0**2-1._RK) )
     124             :         real(RK), parameter     :: X0Sq = X0**2
     125             :         real(RK), parameter     :: PSI_COEF1 = 2._RK**(2._RK/3._RK)
     126             :         real(RK), parameter     :: PSI_COEF2 = -PSI_COEF1 / 252._RK
     127             :         real(RK), parameter     :: PSI_COEF3 = +PSI_COEF1 / 21060._RK
     128             :         real(RK), parameter     :: PSIX0 = X0**0.333333333333333_RK * ( PSI_COEF1 + X0Sq * (PSI_COEF2 + X0Sq*PSI_COEF3) )
     129             :         real(RK), parameter     :: LOG_COEF = log(LS2HC / (OMEGA_DE**0.1666666666666667_RK*OMEGA_DM**0.3333333333333333_RK))
     130          16 :         alpha1          = 1._RK + TWICE_OMEGA_DE_OVER_OMEGA_DM / zplus1**3
     131          16 :         x1              = log( alpha1 + sqrt( alpha1**2 - 1._RK ) )
     132          16 :         x1Sq            = x1**2
     133          16 :         psix1           = x1**0.333333333333333_RK * ( PSI_COEF1 + x1Sq * (PSI_COEF2 + x1Sq*PSI_COEF3) )
     134          16 :         logLumDisWicMpc = LOG_COEF +  log(zplus1*(PSIX0-psix1))
     135          32 :     end function getLogLumDisWicMpc
     136             : 
     137             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     138             : 
     139             :     !> \brief
     140             :     !> Return the approximation to the cosmological luminosity distance.
     141             :     !>
     142             :     !> @param[in]   zplus1 : The redshift plus 1.
     143             :     !>
     144             :     !> \return
     145             :     !> `ldiswickram` : The approximate cosmological luminosity distance.
     146             :     !>
     147             :     !> The approximation is accurate with a relative error of 0.001 for any redshift above 0.1.
     148             :     !>
     149             :     !> \remark
     150             :     !> This function is the same as [getLogLumDisWicMpc()](@ref getloglumdiswicmpc) function,
     151             :     !> except for the fact that it return the natural value, as opposed to the natural logarithm.
     152             :     !> It is kept only for legacy reasons and should not be used in new code.
     153           1 :     pure function ldiswickram(zplus1)
     154             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     155             :         !DEC$ ATTRIBUTES DLLEXPORT :: ldiswickram
     156             : #endif
     157          16 :         use Constants_mod, only: RK
     158             :         implicit none
     159             :         real(RK), intent(in)    :: zplus1
     160           1 :         real(RK)                :: ldiswickram,alpha,alpha0,x,x0,psi,psi0
     161           1 :         alpha=1._RK+2*OMEGA_DE/(OMEGA_DM*zplus1**3)
     162           1 :         alpha0=1._RK+2*OMEGA_DE/OMEGA_DM
     163           1 :         x=log(alpha+sqrt(alpha*alpha-1._RK))
     164           1 :         x0=log(alpha0+sqrt(alpha0*alpha0-1._RK))
     165           1 :         psi=x**0.333333333333333_RK*(1.5874010519682_RK-6.2992105236833e-3_RK*x*x+7.5375168659459e-5_RK*x**4)
     166           1 :         psi0=x0**0.333333333333333_RK*(1.5874010519682_RK-6.2992105236833e-3_RK*x0*x0+7.5375168659459e-5_RK*x0**4)
     167           1 :         ldiswickram=LS2HC*zplus1*(psi0-psi)/(OMEGA_DE**0.1666666666666667_RK*OMEGA_DM**0.333333333333333_RK)
     168           2 :     end function ldiswickram
     169             : 
     170             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     171             : 
     172             :     !> \brief
     173             :     !> Return the cosmological lookback time in GYrs at the given redshift for the assumed cosmological parameters.
     174             :     !>
     175             :     !> @param[in]   zplus1              : The redshift plus 1.
     176             :     !> @param[in]   maxRelativeError    : The maximum tolerance for error in the numerical integration (**optional**, default = 1.e-6).
     177             :     !> @param[in]   nRefinement         : The number of refinements in the Romberg numerical integration (**optional**, default = 5).
     178             :     !> \return
     179             :     !> `lookBackTime` : The cosmological lookback time in GYrs at the given redshift.
     180             :     !>
     181             :     !> \remark
     182             :     !> The rest of the required parameters are taken from the module:
     183             :     !> - `HUBBLE_TIME_GYRS` : The Hubble time in units of Gyrs (Giga Years). It should be a number close to 13.8 Gyrs (Liddle, 2003, Page 57).
     184             :     !>
     185             :     !> \remark
     186             :     !> The integrations are performed using the Romberg integration method.
     187             :     !>
     188             :     !> \author
     189             :     !> Amir Shahmoradi, Sunday 2:31 PM, January 6, 2013, IFS, The University of Texas at Austin.
     190        1306 :     function getLookBackTime(zplus1,maxRelativeError,nRefinement) result(lookBackTime)
     191             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     192             :         !DEC$ ATTRIBUTES DLLEXPORT :: getLookBackTime
     193             : #endif
     194             :         use, intrinsic :: iso_fortran_env, only: output_unit
     195           1 :         use Integration_mod, only: doQuadRombClosed, ErrorMessage
     196             :         use Constants_mod, only: RK, IK
     197             :         implicit none
     198             :         real(RK)    , intent(in)            :: zplus1
     199             :         real(RK)    , intent(in), optional  :: maxRelativeError
     200             :         integer(IK) , intent(in), optional  :: nRefinement
     201             :         real(RK)    , parameter             :: ZPLUS1_MIN = 1._RK
     202             :         real(RK)                            :: lookBackTime
     203        1306 :         real(RK)                            :: maxRelativeErrorDefault, relerr
     204             :         integer(IK)                         :: neval, ierr, nRefinementDefault
     205             : 
     206        1305 :         maxRelativeErrorDefault = 1.e-6_RK; if (present(maxRelativeError)) maxRelativeErrorDefault = maxRelativeError
     207        1306 :         nRefinementDefault = 7_IK; if (present(nRefinement)) nRefinementDefault = nRefinement
     208             : 
     209             :         call doQuadRombClosed   ( getFunc           = getLookBackTimeDensity    &
     210             :                                 , lowerLim          = ZPLUS1_MIN                &
     211             :                                 , upperLim          = zplus1                    &
     212             :                                 , maxRelativeError  = maxRelativeErrorDefault   &
     213             :                                 , nRefinement       = nRefinementDefault        &
     214             :                                 , integral          = lookBackTime              &
     215             :                                 , relativeError     = relerr                    &
     216             :                                 , numFuncEval       = neval                     &
     217             :                                 , ierr              = ierr                      &
     218        1306 :                                 )
     219        1306 :         if (ierr/=0_IK) then
     220             :             ! LCOV_EXCL_START
     221             :             write(output_unit,"(A)") ErrorMessage(ierr)
     222             :             error stop
     223             :             ! LCOV_EXCL_STOP
     224             :         end if
     225        1306 :         lookBackTime = HUBBLE_TIME_GYRS * lookBackTime
     226             : 
     227        1306 :     end function getLookBackTime
     228             : 
     229             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     230             : 
     231             :     !> \brief
     232             :     !> Return the **differential** (w.r.t. `z`) cosmological lookback time in GYrs at the given redshift for the assumed cosmological parameters.
     233             :     !>
     234             :     !> @param[in]   zplus1              : The redshift plus 1.
     235             :     !> \return
     236             :     !> `lookBackTimeDnesity` : The cosmological lookback time in GYrs at the given redshift.
     237             :     !>
     238             :     !> \remark
     239             :     !> The rest of the required parameters are taken from the module:
     240             :     !> - `OMEGA_DE` : the dark energy density,
     241             :     !> - `OMEGA_DM` : the dark matter density.
     242             :     !>
     243             :     !> \remark
     244             :     !> This function could have become an internal function within [getLookBackTime](@ref getlookbacktime).
     245             :     !> However, the library yields segmentation fault error when compiled and run on the Windows Subsystem 
     246             :     !> for Linux Ubuntu with GFortran. As such, it is implemented as an independent function.
     247             :     !>
     248             :     !> \author
     249             :     !> Amir Shahmoradi, Sunday 2:31 PM, January 6, 2013, IFS, The University of Texas at Austin.
     250       61242 :     pure function getLookBackTimeDensity(zplus1) result(lookBackTimeDnesity)
     251        1306 :         use Constants_mod, only: RK
     252             :         implicit none
     253             :         real(RK), intent(in)    :: zplus1
     254             :         real(RK)                :: lookBackTimeDnesity
     255       61242 :         lookBackTimeDnesity = 1._RK / ( zplus1 * sqrt(OMEGA_DM * zplus1**3 + OMEGA_DE) )
     256      122484 :     end function getLookBackTimeDensity
     257             : 
     258             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     259             : 
     260             :     !> \brief
     261             :     !> Return the derivative of the age of the Universe, w.r.t. redshift for a given input redshift + 1.
     262             :     !>
     263             :     !> @param[in] zplus1 : The redshift plus 1.
     264             :     !>
     265             :     !> \return
     266             :     !> `universeAgeDerivative` : The derivative of the age of the Universe, w.r.t. redshift for a given input `zplus1`.
     267             :     !>
     268             :     !> \remark
     269             :     !> The rest of the required parameters are taken from the module:
     270             :     !> - `INVERSE_HUBBLE_CONST`
     271             :     !> - `OMEGA_DE`
     272             :     !> - `OMEGA_DM`
     273             :     !>
     274             :     !> \remark
     275             :     !> The integrations are performed using the Romberg integration method.
     276             :     !>
     277             :     !> \author
     278             :     !> Amir Shahmoradi, Sunday 2:31 PM, January 6, 2013, IFS, The University of Texas at Austin.
     279        1297 :     pure function getUniverseAgeDerivative(zplus1) result(universeAgeDerivative)
     280             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     281             :         !DEC$ ATTRIBUTES DLLEXPORT :: getUniverseAgeDerivative
     282             : #endif
     283       61242 :         use Constants_mod, only: RK
     284             :         implicit none
     285             :         real(RK), intent(in)    :: zplus1
     286             :         real(RK)                :: universeAgeDerivative
     287        1297 :         universeAgeDerivative = INVERSE_HUBBLE_CONST / ( zplus1 * sqrt(OMEGA_DM * zplus1**3 + OMEGA_DE) )
     288        1297 :     end function getUniverseAgeDerivative
     289             : 
     290             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     291             : 
     292             : end module Cosmology_mod

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021