; Computes and return the observed IMFs and CFMs by
; Salpeter (1995), Kroupa(2001), Chabrier (2003,2005)
; they come out such that N(M=1M_sol)=1 
; written by Christoph Federrath 2016

; m_in: mass input in M_sol
function imf_chabrier, m_in, YEAR=year, SYSTEM=system, _EXTRA=extra

    if not keyword_set(year) then year = 2005
    if not keyword_set(system) then system = 0

    n = n_elements(m_in)
    imf = dblarr(n)
    m = dblarr(n)
    m[0:n-1] = m_in[0:n-1]

    ; Chabrier (2005)
    if year eq 2005 then begin
        ind = where(m lt 1d0) ; < 1 M_sol
        if not system then begin
            ; Equation (1) in Chabrier (2005)
            imf[ind] = 9.3d-2 * exp(-((alog10(m[ind]/2.0d-1))^2/(2d0*(5.5d-1)^2)))
        endif else begin
            ; Equation (2) in Chabrier (2005)
            imf[ind] = 7.6d-2 * exp(-((alog10(m[ind]/2.5d-1))^2/(2d0*(5.5d-1)^2)))
        endelse
        ind = where(m ge 1d0) ; >= 1 M_sol
        imf[ind] = 4.1d-2 * m[ind]^(-1.35d0)
    endif

    ; Chabrier (2003)
    if year eq 2003 then begin
        ind = where(m lt 1d0) ; < 1 M_sol
        if not system then begin
            ; Equation (17) in Chabrier (2003)
            imf[ind] = 1.58d-1 * exp(-((alog10(m[ind]/7.9d-2))^2/(2d0*(6.9d-1)^2)))
        endif else begin
            ; Equation (18) in Chabrier (2003)
            imf[ind] = 8.6d-2 * exp(-((alog10(m[ind]/2.2d-1))^2/(2d0*(5.7d-1)^2)))
        endelse
        ind = where(m ge 1d0) ; >= 1 M_sol
        imf[ind] = 4.1d-2 * m[ind]^(-1.35d0)
    endif

    ; normalise Chabrier IMF to 1 at M=1
    ind = where(m ge 1d0)
    imf /= imf[ind[0]]

    return, imf

end ; ===============================================

; m_in: mass input in M_sol
function imf_kroupa, m_in, BINARY_CORRECT=binary_correct, YEAR_2008=year_2008, _EXTRA=extra

    if not keyword_set(binary_correct) then binary_correct = 0
    
    n = n_elements(m_in)
    imf = dblarr(n)
    m = dblarr(n)
    m[0:n-1] = m_in[0:n-1]

    if not keyword_set(year_2008) then begin
        ; here normalised to 1 for M=1
        if not binary_correct then begin
            ; Equations (1) and (2) in Kroupa (2001)
            A1 = 0.5^(-1.0)
            A0 = A1*0.08^(-1.0)
            ind = where(m ge 0.01d0 and m lt 0.08d0) ; 0.01 <= M < 0.08
            imf[ind] = m[ind]^(-0.3+1.0) * A0
            ind = where(m ge 0.08d0 and m lt 0.50d0) ; 0.08 <= M < 0.50
            imf[ind] = m[ind]^(-1.3+1.0) * A1
            ind = where(m ge 0.50d0 and m lt 1.00d0) ; 0.50 <= M < 1.00
            imf[ind] = m[ind]^(-2.3+1.0)
            ind = where(m ge 1.00d0) ; M >= 1.00
            imf[ind] = m[ind]^(-2.3+1.0)
        endif else begin
            ; Equations (1) and (6) in Kroupa (2001)
            A1 = 0.5^(-0.9)
            A0 = A1*0.08^(-1.5)
            ind = where(m ge 0.01d0 and m lt 0.08d0) ; 0.01 <= M < 0.08
            imf[ind] = m[ind]^(-0.3+1.0) * A0
            ind = where(m ge 0.08d0 and m lt 0.50d0) ; 0.08 <= M < 0.50
            imf[ind] = m[ind]^(-1.8+1.0) * A1
            ind = where(m ge 0.50d0 and m lt 1.00d0) ; 0.50 <= M < 1.00
            imf[ind] = m[ind]^(-2.7+1.0)
            ind = where(m ge 1.00d0) ; M >= 1.00
            imf[ind] = m[ind]^(-2.3+1.0)
        endelse
    endif else begin
        ; Equation (1) in Kroupa (2008)
        A1 = 0.5^(-1.0)
        ind = where(m ge 0.08d0 and m lt 0.50d0) ; 0.08 <= M < 0.5
        imf[ind] = m[ind]^(-1.3+1.0) * A1
        ind = where(m ge 0.50d0 and m lt 150d0) ; 0.5 <= M < 150
        imf[ind] = m[ind]^(-2.3+1.0)
    endelse

    return, imf

end ; ===============================================

; m_in: mass input in M_sol
function imf_salpeter, m_in, _EXTRA=extra

    n = n_elements(m_in)
    imf = dblarr(n)
    m = dblarr(n)
    m[0:n-1] = m_in[0:n-1]

    ; Equation (5) in Salpeter (1995)
    ind = where(m ge 10d0^(-0.4) and m le 100d0)
    imf[ind] = 0.03*m[ind]^(-1.35)

    ; normalise Salpeter IMF to 1 at M=1
    ind = where(m ge 1d0)
    imf /= imf[ind[0]]

    return, imf

end ; ===============================================

; get the mean mass of a given IMF
function get_mean_mass, m_in, imf_in
    
    dims = n_elements(m_in)
    m = m_in[0:dims-2]
    dm = (shift(m_in,-1))[0:dims-2] - m
    imf = imf_in[0:dims-2] / m
    mean_mass = total(m*imf*dm,/double) / total(imf*dm,/double)
    return, mean_mass
    
end ; ===============================================

; ############################# MAIN ##################################

; type = ['chabrier2005single', 'chabier2005system']
pro imf_obs, TYPE=type, MASSES=masses, IMF=imf, CMF=cmf, NBINS=nbins, MEAN_MASS=mean_mass, $
    _EXTRA=extra

    if not keyword_set(type) then begin
        print, "imf_obs: keyword TYPE must be one of ['chabrier2005single', 'chabier2005system']"
        stop
    endif
    if not keyword_set(nbins) then nbins = 1000

    ; note that we fix the minimum and maximum mass
    min_mass = 1d-2 ; in M_sol
    max_mass = 1d3 ; in M_sol
    masses = 10d0^(dindgen(nbins)/double(nbins-1)*alog10(max_mass/min_mass)+alog10(min_mass))

    if type eq 'chabrier2005single' then begin
        imf = imf_chabrier(masses,YEAR=2005,SYSTEM=0)
    endif
    if type eq 'chabrier2005system' then begin
        imf = imf_chabrier(masses,YEAR=2005,SYSTEM=1)
    endif
    if type eq 'chabrier2003single' then begin
        imf = imf_chabrier(masses,YEAR=2003,SYSTEM=0)
    endif
    if type eq 'chabrier2003system' then begin
        imf = imf_chabrier(masses,YEAR=2003,SYSTEM=1)
    endif
    
    if type eq 'kroupa2001' then begin
        imf = imf_kroupa(masses,BINARY_CORRECT=0)
    endif
    if type eq 'kroupa2001binary_corrected' then begin
        imf = imf_kroupa(masses,BINARY_CORRECT=1)
    endif
    if type eq 'kroupa2008' then begin
        imf = imf_kroupa(masses,YEAR_2008=1)
    endif

    if type eq 'salpeter1995' then begin
        imf = imf_salpeter(masses)
    endif

    ; get the mean mass for this IMF
    mean_mass = get_mean_mass(masses,imf)

    cmf = total(imf,/double,/cumul) / total(imf,/double)

end ; ===============================================