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 3 : 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 3 : logdvdz = LOG_COEF + twiceLogLumDisMpc - ( 2*logzplus1 + 0.5_RK*log(OMEGA_DM*zplus1**3+OMEGA_DE) )
93 6 : 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 48 : 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 3 : use Constants_mod, only: RK
117 : implicit none
118 : real(RK), intent(in) :: zplus1 ! redshift + 1
119 : real(RK) :: logLumDisWicMpc
120 48 : 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 48 : alpha1 = 1._RK + TWICE_OMEGA_DE_OVER_OMEGA_DM / zplus1**3
131 48 : x1 = log( alpha1 + sqrt( alpha1**2 - 1._RK ) )
132 48 : x1Sq = x1**2
133 48 : psix1 = x1**0.333333333333333_RK * ( PSI_COEF1 + x1Sq * (PSI_COEF2 + x1Sq*PSI_COEF3) )
134 48 : logLumDisWicMpc = LOG_COEF + log(zplus1*(PSIX0-psix1))
135 96 : 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 3 : 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 48 : use Constants_mod, only: RK
158 : implicit none
159 : real(RK), intent(in) :: zplus1
160 3 : real(RK) :: ldiswickram,alpha,alpha0,x,x0,psi,psi0
161 3 : alpha=1._RK+2*OMEGA_DE/(OMEGA_DM*zplus1**3)
162 3 : alpha0=1._RK+2*OMEGA_DE/OMEGA_DM
163 3 : x=log(alpha+sqrt(alpha*alpha-1._RK))
164 3 : x0=log(alpha0+sqrt(alpha0*alpha0-1._RK))
165 3 : psi=x**0.333333333333333_RK*(1.5874010519682_RK-6.2992105236833e-3_RK*x*x+7.5375168659459e-5_RK*x**4)
166 3 : psi0=x0**0.333333333333333_RK*(1.5874010519682_RK-6.2992105236833e-3_RK*x0*x0+7.5375168659459e-5_RK*x0**4)
167 3 : ldiswickram=LS2HC*zplus1*(psi0-psi)/(OMEGA_DE**0.1666666666666667_RK*OMEGA_DM**0.333333333333333_RK)
168 6 : 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 3918 : 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 3 : 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 3918 : real(RK) :: maxRelativeErrorDefault, relerr
204 : integer(IK) :: neval, ierr, nRefinementDefault
205 :
206 3915 : maxRelativeErrorDefault = 1.e-6_RK; if (present(maxRelativeError)) maxRelativeErrorDefault = maxRelativeError
207 3918 : 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 3918 : )
219 3918 : 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 3918 : lookBackTime = HUBBLE_TIME_GYRS * lookBackTime
226 :
227 3918 : 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 183726 : pure function getLookBackTimeDensity(zplus1) result(lookBackTimeDnesity)
251 3918 : use Constants_mod, only: RK
252 : implicit none
253 : real(RK), intent(in) :: zplus1
254 : real(RK) :: lookBackTimeDnesity
255 183726 : lookBackTimeDnesity = 1._RK / ( zplus1 * sqrt(OMEGA_DM * zplus1**3 + OMEGA_DE) )
256 367452 : 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 3891 : 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 183726 : use Constants_mod, only: RK
284 : implicit none
285 : real(RK), intent(in) :: zplus1
286 : real(RK) :: universeAgeDerivative
287 3891 : universeAgeDerivative = INVERSE_HUBBLE_CONST / ( zplus1 * sqrt(OMEGA_DM * zplus1**3 + OMEGA_DE) )
288 3891 : end function getUniverseAgeDerivative
289 :
290 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
291 :
292 : end module Cosmology_mod
|
