/usr/share/gnudatalanguage/astrolib/imf.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 124 125 126 127 128 129 | function imf, mass, expon, mass_range
;+
; NAME:
; IMF
; PURPOSE:
; Compute an N-component power-law logarithmic initial mass function
; EXPLANTION:
; The function is normalized so that the total mass distribution
; equals one solar mass.
;
; CALLING SEQUENCE:
; psi = IMF( mass, expon, mass_range )
;
; INPUTS:
; mass - mass in units of solar masses (scalar or vector)
; Converted to floating point if necessary
; expon - power law exponent, usually negative, scalar or vector
; The number of values in expon equals the number of different
; power-law components in the IMF
; A Saltpeter IMF has a scalar value of expon = -1.35
; mass_range - vector containing the mass upper and lower limits of the
; IMF and masses where the IMF exponent changes. The number
; of values in mass_range should be one more than in expon.
; The values in mass_range should be monotonically increasing.
;
; OUTPUTS
; psi - mass function, number of stars per unit logarithmic mass interval
; evaluated for supplied masses
;
; NOTES:
; The mass spectrum f(m) giving the number of stars per unit mass
; interval is related to psi(m) by m*f(m) = psi(m). The normalization
; condition is that the integral of psi(m) between the upper and lower
; mass limit is unity.
;
; EXAMPLE:
; (1) Print the number of stars per unit mass interval at 3 Msun
; for a Salpeter (expon = -1.35) IMF, with a mass range from
; 0.1 MSun to 110 Msun.
;
; IDL> print, imf(3, -1.35, [0.1, 110] ) / 3
;
; (2) Lequex et al. (1981, A & A 103, 305) describes an IMF with an
; exponent of -0.6 between 0.007 Msun and 1.8 Msun, and an
; exponent of -1.7 between 1.8 Msun and 110 Msun. Plot
; the mass spectrum f(m)
;
; IDL> m = [0.01,0.1,indgen(110) + 1 ] ;Make a mass vector
; IDL> expon = [-0.6, -1.7] ;Exponent Vector
; IDL> mass_range = [ 0.007, 1.8, 110] ;Mass range
; IDL> plot,/xlog,/ylog, m, imf(m, expon, mass_range ) / m
;
; METHOD
; IMF first calculates the constants to multiply the power-law
; components such that the IMF is continuous at the intermediate masses,
; and that the total mass integral is one solar mass. The IMF is then
; calculated for the supplied masses. Also see Scalo (1986, Fund. of
; Cosmic Physics, 11, 1)
;
; PROCEDURES CALLED:
; None
; REVISION HISTORY:
; Written W. Landsman August, 1989
; Set masses LE mass_u rather than LT mass_u August, 1992
; Major rewrite to accept arbitrary power-law components April 1993
; Convert EXPON to float if necessary W. Landsman March 1996
; Remove call to DATATYPE, V5.3 version W. Landsman August 2000
;-
On_error,2
if N_params() LT 3 then begin
print,'Syntax - psi = IMF( mass, expon, mass_range)'
return,-1
endif
Ncomp = N_elements(expon)
if N_elements( mass_range) NE Ncomp + 1 then message, $
'ERROR - Mass Range Vector must have ' + strtrim(Ncomp+1,2) + ' components'
if ( min(mass_range) LE 0 ) then message, $
'ERROR - Mass range Vector must be positive definite'
npts = N_elements(mass)
if ( npts LT 1 ) then begin
message, 'Mass vector (first parameter) has not been defined',/CON
return,0
endif
if size(mass,/TNAME) NE 'DOUBLE' then mass = float(mass) ;Make sure not integer
if size(expon,/TNAME) NE 'DOUBLE' then expon = float(expon)
; Get normalization constants for supplied power-law exponents
integ = fltarr(ncomp)
;Compute the unnormalized integral over each power law section
for i = 0, Ncomp-1 do begin
if ( expon[i] NE -1 ) then integ[i] = $
(mass_range[i+1]^(1+expon[i]) - mass_range[i]^(1+expon[i]))/(1+expon[i]) $
else integ[i] = alog(mass_range[i+1]/mass_range[i])
endfor
; Insure continuity where the power law functions meet
joint = fltarr(ncomp)
joint[0] = 1
if ncomp GT 1 then for i = 1,ncomp-1 do begin
joint[i] = joint[i-1]*mass_range[i]^( expon[i-1] - expon[i] )
endfor
norm = fltarr(ncomp)
norm[0] = 1./ total(integ*joint)
if ncomp GT 1 then for i = 1,ncomp-1 do norm[i] = norm[0]*joint[i]
f = mass*0.
for i = 0, Ncomp-1 do begin
test = where( (mass GT mass_range[i]) and (mass LE mass_range[i+1]), Ntest )
if ( Ntest GT 0 ) then f[test] = norm[i]*mass[test]^(expon[i])
endfor
return,f
end
|