博文
黑洞质量计算（IDL程序）
|||
Function mcgill08bhmass, z, flux, velocity, $
H0=H0, k=k, Lambda0=Lambda0, Omega_m=Omega_m, q0=q0, $
slope_opt=slope_opt, k_cor=k_cor, $
flux_select=flux_select, velocity_select=velocity_select
;+
;NAME:
; mcgill08bhmass
;PURPOSE:
; Calculate the Black Hole mass
;CALLING SEQUENCE:
; lgMBH=mcgill08bhmass(z, flux, velocity)
;INPUT:
; z
; flux
; velocity
;OPTIONAL KEYWORD INPUT:
; H0 --- Hubble parameter in km/s/Mpc, default is 70
;
; No more thbn 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.5.
; slope_opt --- power law index
; flux ~ (frequency)^(- slope_opt)
; default it,it is 0.
; flux_select
; velocity_select
;OPTIONAL KEYWORD OUTPUT:
; /k_cor ---- if set it,input flux will be not through
; K_correction. default it, input flux will be through
; K_correction.
;METHOD:
; McGill astro-ph 0710.1839v2
;REVISION HISTOTY:
; Original by DL.Wang,12-Dec-2008
;-
;
if not keyword_set(flux_select) then begin
print,' flux selection for luminosity estimator (Unit:erg/s/cm^2)'
print,' 3000 : 3000A continue flux'
print,' 5100n : 5100A nuclei continue flux'
print,' 5100t : 5100A total continue flux'
print,' HB : H Beta emission line flux'
print,' HA : H Alpha emission line flux'
flux_select=''
read,flux_select
endif
if not keyword_set(velocity_select) then begin
print,' velocity selection for velocity estimator (Unit:km/s)'
print,' Mg2s : MgII emission line dispersion'
print,' HBs : H Beta emission line dispersion'
print,' HAs : H Alpha emission line dispersion'
print,' mg2f : MgII emission line FWHM'
print,' HBf : H Beta emission line FWHM'
print,' HAf : H Alpha emission line FWHM'
velocity_select=''
read,velocity_select
endif
case flux_select of
'3000':begin
b=0.470 & lumcal=1.0D44;erg/s
end
'5100n':begin
b=0.518 & lumcal=1.0D44;erg/s
end
'5100t':begin
b=0.690 & lumcal=1.0D44;erg/s
end
'HB':begin
b=0.560 & lumcal=1.0D42;erg/s
end
'HA':begin
b=0.550 & lumcal=1.0D42;erg/s
end
endcase
c=2.0
if flux_select eq '3000' then begin
case velocity_select of
'Mg2s':begin
a=7.207 & aerr=0.052
end
'HBs':begin
a=7.458 & aerr=0.027
end
'HAs':begin
a=7.588 & aerr=0.061
end
'Mg2f':begin
a=6.767 & aerr=0.055
end
'HBf':begin
a=6.803 & aerr=0.069
end
'HAf':begin
a=6.986 & aerr=0.064
end
endcase
endif
if flux_select eq '5100n' then begin
case velocity_select of
'Mg2s':begin
a=7.429 & aerr=0.039
end
'HBs':begin
a=7.680 & aerr=0.021
end
'HAs':begin
a=7.810 & aerr=0.048
end
'Mg2f':begin
a=6.990 & aerr=0.045
end
'HBf':begin
a=7.026 & aerr=0.056
end
'HAf':begin
a=7.209 & aerr=0.055
end
endcase
endif
if flux_select eq '5100t' then begin
case velocity_select of
'Mg2s':begin
a=7.133 & aerr=0.045
end
'HBs':begin
a=7.383 & aerr=0.028
end
'HAs':begin
a=7.514 & aerr=0.060
end
'Mg2f':begin
a=6.693 & aerr=0.053
end
'HBf':begin
a=6.729 & aerr=0.071
end
'HAf':begin
a=6.912 & aerr=0.069
end
endcase
endif
if flux_select eq 'HB' then begin
case velocity_select of
'Mg2s':begin
a=7.150 & aerr=0.053
end
'HBs':begin
a=7.401 & aerr=0.038
end
'HAs':begin
a=7.532 & aerr=0.074
end
'Mg2f':begin
a=6.711 & aerr=0.059
end
'HBf':begin
a=6.747 & aerr=0.077
end
'HAf':begin
a=6.930 & aerr=0.071
end
endcase
endif
if flux_select eq 'HA' then begin
case velocity_select of
'Mg2s':begin
a=6.824 & aerr=0.051
end
'HBs':begin
a=7.074 & aerr=0.042
end
'HAs':begin
a=7.205 & aerr=0.075
end
'Mg2f':begin
a=6.384 & aerr=0.056
end
'HBf':begin
a=6.420 & aerr=0.079
end
'HAf':begin
a=6.603 & aerr=0.073
end
endcase
endif
velcal=1.0D3;km/s
cosmo_param,Omega_m,Lambda0,k,q0
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)'
if not keyword_set(slope_opt) then slope_opt=1.0
k_correct=(1.0+z)^(-1.0+slope_opt)
if keyword_set(k_cor) then begin
Flux_kcor=Flux*k_correct
endif else begin
Flux_kcor=Flux
endelse
dL=lumdist(z,H0=H0,k=k,Omega_m=Omega_m,Lambda0=Lambda0,q0=q0)
factor=3.085678D+24
dL0=dL*factor ;transform Mpc to cm
Lum0=4.0*!DPI*dL0*dL0*Flux_kcor ;erg/s
lgMbh=a+c*alog10(velocity/velcal)+b*alog10(Lum0/lumcal)
return,lgMbh
End
https://blog.sciencenet.cn/blog-456360-658896.html
上一篇：喷流强度计算（IDL程序）
下一篇：高效办公书籍
