/usr/share/gnudatalanguage/astrolib/lumdist.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 | ;+
; NAME:
; LUMDIST
;
; PURPOSE:
; Calculate luminosity distance (in Mpc) of an object given its redshift
; EXPLANATION:
; The luminosity distance in the Friedmann-Robertson-Walker model is
; taken from Caroll, Press, and Turner (1992, ARAA, 30, 499), p. 511
; Uses a closed form (Mattig equation) to compute the distance when the
; cosmological constant is zero. Otherwise integrates the function using
; QSIMP.
; CALLING SEQUENCE:
; result = lumdist(z, [H0 = , k = , Omega_M =, Lambda0 = , q0 = ,/SILENT])
;
; INPUTS:
; z = redshift, positive scalar or vector
;
; OPTIONAL KEYWORD INPUTS:
; /SILENT - If set, the program will not display adopted cosmological
; parameters at the terminal.
; H0: Hubble parameter in km/s/Mpc, default is 70
;
; No more than two of the following four parameters should be
; specified. None of them need be specified -- the adopted defaults
; are given.
; k - curvature constant, normalized to the closure density. Default is
; 0, indicating a flat universe
; Omega_m - Matter density, normalized to the closure density, default
; is 0.3. Must be non-negative
; Lambda0 - Cosmological constant, normalized to the closure density,
; default is 0.7
; q0 - Deceleration parameter, numeric scalar = -R*(R'')/(R')^2, default
; is -0.55
;
; OUTPUTS:
; The result of the function is the luminosity distance (in Mpc) for each
; input value of z.
;
; EXAMPLE:
; (1) Plot the distance of a galaxy in Mpc as a function of redshift out
; to z = 5.0, assuming the default cosmology (Omega_m=0.3, Lambda = 0.7,
; H0 = 70 km/s/Mpc)
;
; IDL> z = findgen(50)/10.
; IDL> plot,z,lumdist(z),xtit='z',ytit='Distance (Mpc)'
;
; Now overplot the relation for zero cosmological constant and
; Omega_m=0.3
; IDL> oplot,z,lumdist(z,lambda=0,omega=0.3),linestyle=1
; COMMENTS:
; (1) Integrates using the IDL Astronomy Version procedure QSIMP. (The
; intrinsic IDL QSIMP function is not called because of its ridiculous
; restriction that only scalar arguments can be passed to the integrating
; function.)
; (2) Can fail to converge at high redshift for closed universes with
; non-zero lambda. This can presumably be fixed by replacing QSIMP with
; an integrator that can handle a singularity
; PROCEDURES CALLED:
; COSMO_PARAM, QSIMP
; REVISION HISTORY:
; Written W. Landsman Raytheon ITSS April 2000
; Avoid integer overflow for more than 32767 redshifts July 2001
; Use double precision J. Moustakas/W. Landsman April 2008
;-
function ldist, z, q0 = q0, lambda0 = lambda0
term1 = (1.+z)^2
term2 = 1.+2.*(q0+lambda0)*z
term3 = z*(2.+z)*lambda0
denom = (term1*term2 - term3)
out = z*0.
good = where(denom GT 0.0, Ngood)
if Ngood GT 0 then out[good] = 1./sqrt(denom[good])
return, out
end
FUNCTION lumdist, z, h0=h0, k = k, Lambda0 = lambda0, Omega_m = Omega_m, $
q0 = q0, Silent = silent
compile_opt idl2
if N_params() eq 0 then begin
print,'Syntax: result = lumdist(z, H0 = ,k=, Lambda0 = ])'
print,'Returns luminosity distance in Mpc'
return, 0.
endif
n = N_elements(z)
cosmo_param,Omega_m,Lambda0, k, q0
; Check keywords
c = 2.99792458D5 ; speed of light in km/s
if N_elements(H0) EQ 0 then H0 = 70
if not keyword_set(silent) then $
print,'LUMDIST: H0:', h0, ' Omega_m:', omega_m, ' Lambda0',lambda0, $
' q0: ',q0, ' k: ', k, f='(A,I3,A,f5.2,A,f5.2,A,f5.2,A,F5.2)'
; For the case of Lambda = 0, we use the closed form from equation 5.238 of
; Astrophysical Formulae (Lang 1998). This avoids terms that almost cancel
; at small q0*z better than the more familiar Mattig formula.
;
if lambda0 EQ 0 then begin
denom = sqrt(1+2*q0*z) + 1 + q0*z
dlum = (c*z/h0)*(1 + z*(1-q0)/denom)
return,dlum
; For non-zero lambda
endif else begin
dlum = z*0.0
for i=0L,N-1 do begin
if z[i] LE 0.0 then dlum[i] = 0.0 else begin
qsimp,'LDIST',0,z[i], lz,q0 = q0, Lambda0 = lambda0
dlum[i] = lz
endelse
endfor
if k GT 0 then $
dlum = sinh(sqrt(k)*dlum)/sqrt(k) $
else if k LT 0 then $
dlum = sin(sqrt(-k)*dlum)/sqrt(-k) > 0
return, c*(1+z)*dlum/h0
endelse
end
|