Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This file contains implementations of procedures [pm_cosmology](@ref pm_cosmology).
19 : !>
20 : !> \author
21 : !> \AmirShahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
22 :
23 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%
26 : #if getSizeUnivNormed_ENABLED
27 : !%%%%%%%%%%%%%%%%%%%%%%%%
28 :
29 : integer(IK) :: sindex(1000), neval_def, nint_def, err_def
30 : real(RKC) :: sinfo(4, 1000), abserr
31 : !reltol = sqrt(epsilon(0._RKC))
32 1508 : if (zplus1 < huge(zplus1)) then
33 : err_def = getQuadErr( getFunc = getIntegrand & ! LCOV_EXCL_LINE
34 : , lb = zplus1 & ! LCOV_EXCL_LINE
35 : , ub = pinf & ! LCOV_EXCL_LINE
36 : , abstol = 0._RKC & ! LCOV_EXCL_LINE
37 : , reltol = reltol & ! LCOV_EXCL_LINE
38 : , qrule = GK21 & ! LCOV_EXCL_LINE
39 : , help = weps & ! LCOV_EXCL_LINE
40 : , integral = sizeUnivNormed & ! LCOV_EXCL_LINE
41 : , abserr = abserr & ! LCOV_EXCL_LINE
42 : , sinfo = sinfo & ! LCOV_EXCL_LINE
43 : , sindex = sindex & ! LCOV_EXCL_LINE
44 : , neval = neval_def & ! LCOV_EXCL_LINE
45 : , nint = nint_def & ! LCOV_EXCL_LINE
46 1508 : )
47 1508 : if (present(err)) then
48 0 : err = err_def
49 1508 : elseif (err_def /= 0_IK) then
50 : error stop SK_"@getQuadErr() failed with err = "//getStr(err_def) ! LCOV_EXCL_LINE
51 : end if
52 1508 : if (present(neval)) neval = neval_def
53 : else
54 0 : sizeUnivNormed = 0._RKC
55 0 : if (present(err)) err = 0_IK
56 0 : if (present(neval)) neval = 0_IK
57 : end if
58 :
59 : contains
60 :
61 377916 : function getIntegrand(zplus1) result(integrand)
62 : real(RKC) , intent(in) :: zplus1
63 : real(RKC) :: integrand
64 : #if Z_ENABLED
65 378 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1)))
66 : #elif ZML_ENABLED
67 377454 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL)))
68 : #elif ZMLR_ENABLED
69 42 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR)))
70 : #elif ZMLRK_ENABLED
71 42 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR, omegaK)))
72 : #else
73 : #error "Unrecognized interface."
74 : #endif
75 377916 : end function
76 :
77 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
78 : #elif getDisLookbackNormed_ENABLED
79 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 :
81 : integer(IK) :: sindex(1000), neval_def, nint_def, err_def
82 : real(RKC) :: sinfo(4, 1000), abserr
83 1684 : if (zplus1 > 1._RKC) then
84 : err_def = getQuadErr( getFunc = getIntegrand & ! LCOV_EXCL_LINE
85 : , lb = 1._RKC & ! LCOV_EXCL_LINE
86 : , ub = zplus1 & ! LCOV_EXCL_LINE
87 : , abstol = 0._RKC & ! LCOV_EXCL_LINE
88 : , reltol = reltol & ! LCOV_EXCL_LINE
89 : , qrule = GK21 & ! LCOV_EXCL_LINE
90 : , integral = disLookbackNormed & ! LCOV_EXCL_LINE
91 : , abserr = abserr & ! LCOV_EXCL_LINE
92 : , sinfo = sinfo & ! LCOV_EXCL_LINE
93 : , sindex = sindex & ! LCOV_EXCL_LINE
94 : , neval = neval_def & ! LCOV_EXCL_LINE
95 : , nint = nint_def & ! LCOV_EXCL_LINE
96 1619 : )
97 1619 : if (present(err)) then
98 56 : err = err_def
99 1563 : elseif (err_def /= 0_IK) then
100 : error stop SK_"@getQuadErr() failed with err = "//getStr(err_def) ! LCOV_EXCL_LINE
101 : end if
102 1619 : if (present(neval)) neval = neval_def
103 : else
104 65 : disLookbackNormed = 0._RKC
105 65 : if (present(err)) err = 0_IK
106 65 : if (present(neval)) neval = 0_IK
107 : end if
108 :
109 : contains
110 :
111 211995 : function getIntegrand(zplus1) result(integrand)
112 : real(RKC) , intent(in) :: zplus1
113 : real(RKC) :: integrand
114 : #if Z_ENABLED
115 9933 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1)))
116 : #elif ZML_ENABLED
117 166866 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL)))
118 : #elif ZMLR_ENABLED
119 9954 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR)))
120 : #elif ZMLRK_ENABLED
121 25242 : integrand = 1._RKC / (zplus1 * sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR, omegaK)))
122 : #else
123 : #error "Unrecognized interface."
124 : #endif
125 211995 : end function
126 :
127 : !%%%%%%%%%%%%%%%%%%%%%%
128 : #elif getDisComNormed_ENABLED
129 : !%%%%%%%%%%%%%%%%%%%%%%
130 :
131 : integer(IK) :: sindex(1000), neval_def, nint_def, err_def
132 : real(RKC) :: sinfo(4, 1000), abserr
133 : !reltol = sqrt(epsilon(0._RKC))
134 76401 : if (zplus1 > 1._RKC) then
135 : err_def = getQuadErr( getFunc = getHubbleParamNormedInv & ! LCOV_EXCL_LINE
136 : , lb = 1._RKC & ! LCOV_EXCL_LINE
137 : , ub = zplus1 & ! LCOV_EXCL_LINE
138 : , abstol = 0._RKC & ! LCOV_EXCL_LINE
139 : , reltol = reltol & ! LCOV_EXCL_LINE
140 : , qrule = GK21 & ! LCOV_EXCL_LINE
141 : , integral = disComNormed & ! LCOV_EXCL_LINE
142 : , abserr = abserr & ! LCOV_EXCL_LINE
143 : , sinfo = sinfo & ! LCOV_EXCL_LINE
144 : , sindex = sindex & ! LCOV_EXCL_LINE
145 : , neval = neval_def & ! LCOV_EXCL_LINE
146 : , nint = nint_def & ! LCOV_EXCL_LINE
147 76138 : )
148 76138 : if (present(err)) then
149 1088 : err = err_def
150 75050 : elseif (err_def /= 0_IK) then
151 : error stop SK_"@getQuadErr() failed with err = "//getStr(err_def) ! LCOV_EXCL_LINE
152 : end if
153 76138 : if (present(neval)) neval = neval_def
154 : else
155 263 : disComNormed = 0._RKC
156 263 : if (present(err)) err = 0_IK
157 263 : if (present(neval)) neval = 0_IK
158 : end if
159 :
160 : contains
161 :
162 3667020 : function getHubbleParamNormedInv(zplus1) result(hubbleParamNormedInv)
163 : real(RKC) , intent(in) :: zplus1
164 : real(RKC) :: hubbleParamNormedInv
165 : #if Z_ENABLED
166 2546628 : hubbleParamNormedInv = 1._RKC / sqrt(getHubbleParamNormedSq(zplus1))
167 : #elif ZML_ENABLED
168 965664 : hubbleParamNormedInv = 1._RKC / sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL))
169 : #elif ZMLR_ENABLED
170 48678 : hubbleParamNormedInv = 1._RKC / sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR))
171 : #elif ZMLRK_ENABLED
172 106050 : hubbleParamNormedInv = 1._RKC / sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR, omegaK))
173 : #else
174 : #error "Unrecognized interface."
175 : #endif
176 3667020 : end function
177 :
178 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
179 : #elif getDisComTransNormed_ENABLED
180 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
181 :
182 : #if Z_ENABLED
183 63729 : disComTransNormed = getDisComNormed(zplus1, reltol, neval, err)
184 : #elif ZML_ENABLED
185 9476 : disComTransNormed = getDisComNormed(zplus1, omegaM, omegaL, reltol, neval, err)
186 : #elif ZMLR_ENABLED
187 476 : disComTransNormed = getDisComNormed(zplus1, omegaM, omegaL, omegaR, reltol, neval, err)
188 : #elif ZMLRK_ENABLED
189 : real(RKC) , parameter :: POS_EPS = epsilon(0._RKC)
190 : real(RKC) , parameter :: NEG_EPS = -POS_EPS
191 876 : disComTransNormed = getDisComNormed(zplus1, omegaM, omegaL, omegaR, omegaK, reltol, neval, err)
192 876 : if (omegaK > POS_EPS) then
193 306 : disComTransNormed = sinh(disComTransNormed * sqrtAbsOmegaK) / sqrtAbsOmegaK
194 570 : elseif (omegaK < NEG_EPS) then
195 266 : disComTransNormed = sin(disComTransNormed * sqrtAbsOmegaK) / sqrtAbsOmegaK
196 : end if
197 : #else
198 : #error "Unrecognized interface."
199 : #endif
200 :
201 : !%%%%%%%%%%%%%%%%%%%%%%
202 : #elif getDisAngNormed_ENABLED
203 : !%%%%%%%%%%%%%%%%%%%%%%
204 :
205 : #if Z_ENABLED
206 18 : disAngNormed = getDisComTransNormed(zplus1, reltol, neval, err) / zplus1
207 : #elif ZML_ENABLED
208 1518 : disAngNormed = getDisComTransNormed(zplus1, omegaM, omegaL, reltol, neval, err) / zplus1
209 : #elif ZMLR_ENABLED
210 18 : disAngNormed = getDisComTransNormed(zplus1, omegaM, omegaL, omegaR, reltol, neval, err) / zplus1
211 : #elif ZMLRK_ENABLED
212 50 : disAngNormed = getDisComTransNormed(zplus1, omegaM, omegaL, omegaR, omegaK, sqrtAbsOmegaK, reltol, neval, err) / zplus1
213 : #else
214 : #error "Unrecognized interface."
215 : #endif
216 :
217 : !%%%%%%%%%%%%%%%%%%%%%%
218 : #elif getDisLumNormed_ENABLED
219 : !%%%%%%%%%%%%%%%%%%%%%%
220 :
221 : #if Z_ENABLED
222 18 : disLumNormed = getDisComTransNormed(zplus1, reltol, neval, err) * zplus1
223 : #elif ZML_ENABLED
224 1518 : disLumNormed = getDisComTransNormed(zplus1, omegaM, omegaL, reltol, neval, err) * zplus1
225 : #elif ZMLR_ENABLED
226 18 : disLumNormed = getDisComTransNormed(zplus1, omegaM, omegaL, omegaR, reltol, neval, err) * zplus1
227 : #elif ZMLRK_ENABLED
228 50 : disLumNormed = getDisComTransNormed(zplus1, omegaM, omegaL, omegaR, omegaK, sqrtAbsOmegaK, reltol, neval, err) * zplus1
229 : #else
230 : #error "Unrecognized interface."
231 : #endif
232 :
233 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
234 : #elif getDisComTransNormedWU10_ENABLED
235 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
236 :
237 : #if Z_ENABLED
238 : real(RKC), parameter :: omegaL = real(OMEGA_L, RKC)
239 : real(RKC), parameter :: omegaM = real(OMEGA_M, RKC)
240 : #define SET_PARAMETER(x) parameter(x)
241 : #elif ZML_ENABLED
242 : #define SET_PARAMETER(x) x
243 : #else
244 : #error "Unrecognized interface."
245 : #endif
246 : real(RKC) :: alpha1, x1, x1Sq, psix1
247 : real(RKC), parameter :: PSI_COEF1 = 2._RKC**(2._RKC/3._RKC)
248 : real(RKC), parameter :: PSI_COEF2 = -PSI_COEF1 / 252._RKC
249 : real(RKC), parameter :: PSI_COEF3 = +PSI_COEF1 / 21060._RKC
250 : real(RKC), parameter :: ONE_THIRD = 1._RKC / 3._RKC
251 : real(RKC), parameter :: ONE_SIXTH = 1._RKC / 6._RKC
252 : real(RKC) :: twiceOmegaL2OmegaM
253 : real(RKC) :: alpha0
254 : real(RKC) :: x0
255 : real(RKC) :: x0Sq
256 : real(RKC) :: psiX0
257 1582 : SET_PARAMETER(twiceOmegaL2OmegaM = 2._RKC * real(omegaL / omegaM, RKC)) ! fpp
258 1582 : SET_PARAMETER(alpha0 = 1._RKC + twiceOmegaL2OmegaM) ! fpp
259 1582 : SET_PARAMETER(x0 = log(alpha0 + sqrt(alpha0**2 - 1._RKC))) ! fpp
260 1582 : SET_PARAMETER(x0Sq = x0**2) ! fpp
261 1582 : SET_PARAMETER(psiX0 = x0**ONE_THIRD * (PSI_COEF1 + x0Sq * (PSI_COEF2 + x0Sq * PSI_COEF3))) ! fpp
262 2664 : CHECK_ASSERTION(__LINE__, 1._RKC <= zplus1, SK_"@getlogDisLum(): The condition `zplus1 >= 0.` must hold. zplus1 = "//getStr(zplus1)) ! fpp
263 5828 : CHECK_ASSERTION(__LINE__, abs(1._RKC - omegaM - omegaL) <= epsilon(0._RKC), SK_"@getlogDisLum(): The condition `1._RKC - omegaM - omegaL == 0._RKC` must hold. omegaM, omegaL = "//getStr([omegaM, omegaL])) ! fpp
264 2664 : alpha1 = 1._RKC + twiceOmegaL2OmegaM / zplus1**3
265 2664 : x1 = log(alpha1 + sqrt(alpha1**2 - 1._RKC))
266 2664 : x1Sq = x1**2
267 2664 : psix1 = x1**ONE_THIRD * (PSI_COEF1 + x1Sq * (PSI_COEF2 + x1Sq * PSI_COEF3))
268 2664 : disComTransNormedWU10 = (psiX0 - psix1) / real(omegaL**ONE_SIXTH * omegaM**ONE_THIRD, RKC)
269 :
270 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
271 : #elif getlogDisLookback_ENABLED || getlogDisLookback_ENABLED
272 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
273 :
274 : #error "Unrecognized interface."
275 :
276 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
277 : #elif getHubbleParamNormedSq_ENABLED
278 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
279 :
280 : #if Z_ENABLED
281 : real(RKC) , parameter :: omegaL = real(OMEGA_L, RKC)
282 : real(RKC) , parameter :: omegaM = real(OMEGA_M, RKC)
283 2626571 : hubbleParamNormedSq = omegaL + zplus1**3 * omegaM
284 : #elif ZML_ENABLED
285 1519360 : CHECK_ASSERTION(__LINE__, 1._RKC <= zplus1, SK_"@getHubbleParamNormedSq(): The condition `1._RKC <= zplus1` must hold. zplus1 = "//getStr(zplus1)) ! fpp
286 4558080 : CHECK_ASSERTION(__LINE__, abs(1._RKC - omegaM - omegaL) < EPS, SK_"@getHubbleParamNormedSq(): The condition `omegaM + omegaL = 1._RKC` must hold. omegaM, omegaL = "//getStr([omegaM, omegaL])) ! fpp
287 1519360 : hubbleParamNormedSq = omegaL + zplus1**3 * omegaM
288 : #elif ZMLR_ENABLED
289 64050 : CHECK_ASSERTION(__LINE__, 1._RKC <= zplus1, SK_"@getHubbleParamNormedSq(): The condition `1._RKC <= zplus1` must hold. zplus1 = "//getStr(zplus1)) ! fpp
290 256200 : CHECK_ASSERTION(__LINE__, abs(1._RKC - omegaM - omegaL - omegaR) < EPS, SK_"@getHubbleParamNormedSq(): The condition `omegaM + omegaL + omegaR = 1._RKC` must hold. omegaM, omegaL, omegaR = "//getStr([omegaM, omegaL, omegaR])) ! fpp
291 64050 : hubbleParamNormedSq = omegaL + zplus1**2 * (zplus1 * (omegaM + zplus1 * omegaR))
292 : #elif ZMLRK_ENABLED
293 136790 : CHECK_ASSERTION(__LINE__, 1._RKC <= zplus1, SK_"@getHubbleParamNormedSq(): The condition `1._RKC <= zplus1` must hold. zplus1 = "//getStr(zplus1)) ! fpp
294 683950 : CHECK_ASSERTION(__LINE__, abs(1._RKC - omegaM - omegaL - omegaR - omegaK) < EPS, SK_"@getHubbleParamNormedSq(): The condition `omegaM + omegaL + omegaR + omegaK = 1._RKC` must hold. omegaM, omegaL, omegaR, omegaK = "//getStr([omegaM, omegaL, omegaR, omegaK])) ! fpp
295 136790 : hubbleParamNormedSq = omegaL + zplus1**2 * (omegaK + zplus1 * (omegaM + zplus1 * omegaR))
296 : #else
297 : #error "Unrecognized interface."
298 : #endif
299 :
300 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
301 : #elif getVolComDiffNormed_ENABLED
302 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
303 :
304 : #if Z_ENABLED
305 63415 : volComDiffNormed = getDisComTransNormed(zplus1, reltol, neval, err)
306 63415 : volComDiffNormed = volComDiffNormed**2 / sqrt(getHubbleParamNormedSq(zplus1))
307 : #elif ZML_ENABLED
308 1662 : volComDiffNormed = getDisComTransNormed(zplus1, omegaM, omegaL, reltol, neval, err)
309 1662 : volComDiffNormed = volComDiffNormed**2 / sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL))
310 : #elif ZMLR_ENABLED
311 162 : volComDiffNormed = getDisComTransNormed(zplus1, omegaM, omegaL, omegaR, reltol, neval, err)
312 162 : volComDiffNormed = volComDiffNormed**2 / sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR))
313 : #elif ZMLRK_ENABLED
314 162 : volComDiffNormed = getDisComTransNormed(zplus1, omegaM, omegaL, omegaR, omegaK, sqrtAbsOmegaK, reltol, neval, err)
315 162 : volComDiffNormed = volComDiffNormed**2 / sqrt(getHubbleParamNormedSq(zplus1, omegaM, omegaL, omegaR, omegaK))
316 : #else
317 : #error "Unrecognized interface."
318 : #endif
319 :
320 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
321 : #elif setVolComDiffNormed_ENABLED
322 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
323 :
324 3388 : volComDiffNormed = disComTransNormedSq / hubbleParamNormed
325 :
326 : !%%%%%%%%%%%%%%%%%%%%%%
327 : #elif getVolComNormed_ENABLED
328 : !%%%%%%%%%%%%%%%%%%%%%%
329 :
330 : #if D_ENABLED
331 : real(RKC), parameter :: COEF = 4._RKC * acos(-1._RKC) / 3._RKC
332 1746 : volComNormed = COEF * disComTransNormed**3
333 : #elif DOS_ENABLED
334 : real(RKC), parameter :: COEF = 2._RKC * acos(-1._RKC)
335 194 : if (omegaK > 0._RKC) then
336 48 : volComNormed = COEF * (disComTransNormed * sqrt(1._RKC + omegaK * disComTransNormed**2) - asinh(sqrtAbsOmegaK * disComTransNormed) / sqrtAbsOmegaK) / omegaK
337 146 : elseif (omegaK < 0._RKC) then
338 50 : volComNormed = COEF * (disComTransNormed * sqrt(1._RKC + omegaK * disComTransNormed**2) - asin (sqrtAbsOmegaK * disComTransNormed) / sqrtAbsOmegaK) / omegaK
339 : else
340 96 : volComNormed = getVolComNormed(disComTransNormed)
341 : end if
342 : #else
343 : #error "Unrecognized interface."
344 : #endif
345 :
346 : #else
347 : !%%%%%%%%%%%%%%%%%%%%%%%%
348 : #error "Unrecognized interface."
349 : !%%%%%%%%%%%%%%%%%%%%%%%%
350 : #endif
351 : #undef SET_PARAMETER
352 : ! use pm_cosmology, only: LOG_LIGHT_SPEED
353 : !#if getLogDVDZ_ENABLED || getLogDVDZ_D1_ENABLED
354 : ! use pm_cosmology, only: omegaM => OMEGA_M
355 : ! use pm_cosmology, only: logH0 => LOG_HUBBLE_CONST
356 : !#define SET_PARAMETER(x) parameter(x)
357 : !#elif getLogDVDZ_HM_ENABLED || getLogDVDZ_HM_D1_ENABLED
358 : !#define SET_PARAMETER(x) x
359 : !#else
360 : !#error "Unrecognized interface."
361 : !#endif
362 : !
363 : !#if getLogDVDZ_D1_ENABLED || getLogDVDZ_HM_D1_ENABLED
364 : ! integer(IK) :: i
365 : !#define GET_ELEMENT(Array) Array(i)
366 : !#elif getLogDVDZ_ENABLED || getLogDVDZ_HM_ENABLED
367 : !#define GET_ELEMENT(Array) Array
368 : !#define ALL
369 : !#else
370 : !#error "Unrecognized interface."
371 : !#endif
372 : !
373 : ! real(RKC) :: logCoef, omegaDE
374 : ! SET_PARAMETER(omegaDE = 1._RKC - real(omegaM,RKC))
375 : ! SET_PARAMETER(logCoef = log(4._RKC * acos(-1._RKC)) + real(LOG_LIGHT_SPEED,RKC) - real(logH0,RKC))
376 : !
377 : !#if CHECK_ENABLED
378 : ! block
379 : ! use pm_err, only: setAsserted
380 : ! use pm_val2str, only: getStr
381 : ! intrinsic :: size
382 : ! call setAsserted(ALL(zplus1 >= 1._RKC), MODULE_NAME//SK_"@getLogDVDZ(): `zplus1 >= 0.` must hold. zplus1 = "//getStr(zplus1))
383 : ! call setAsserted(ALL(logzplus1 >= 0._RKC), MODULE_NAME//SK_"@getLogDVDZ(): `logzplus1 >= 0.` must hold. logzplus1 = "//getStr(logzplus1))
384 : !#if getLogDVDZ_D1_ENABLED || getLogDVDZ_HM_D1_ENABLED
385 : ! call setAsserted(size(zplus1,kind=IK) == size(logzplus1,kind=IK), MODULE_NAME//SK_"@getLogDVDZ(): `size(zplus1) == size(logzplus1)` must hold. size(zplus1), size(logzplus1) = "//getStr([size(zplus1,kind=IK), size(logzplus1,kind=IK)]))
386 : ! call setAsserted(size(zplus1,kind=IK) == size(TwiceLogLumDisMpc,kind=IK), MODULE_NAME//SK_"@getLogDVDZ(): `size(zplus1) == size(TwiceLogLumDisMpc)` must hold. size(zplus1), size(TwiceLogLumDisMpc) = "//getStr([size(zplus1,kind=IK), size(TwiceLogLumDisMpc,kind=IK)]))
387 : !#endif
388 : ! end block
389 : !#endif
390 : !
391 : !
392 : !#if getLogDVDZ_D1_ENABLED || getLogDVDZ_HM_D1_ENABLED
393 : ! do concurrent (i = 1 : size(zplus1, kind = IK))
394 : !#endif
395 : ! GET_ELEMENT(LogDVDZ) = logCoef + GET_ELEMENT(TwiceLogLumDisMpc) - 3 * GET_ELEMENT(logzplus1) - 0.5_RKC * log(real(omegaM,RKC) * GET_ELEMENT(zplus1)**3 + omegaDE)
396 : !#if getLogDVDZ_D1_ENABLED || getLogDVDZ_HM_D1_ENABLED
397 : ! end do
398 : !#endif
399 : !
400 : !#undef SET_PARAMETER
401 : !#undef GET_ELEMENT
402 : !#undef ALL
|
