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