Add imf function (#19)

* Intial work, received values do not confirm with IDL values currently

* Add imf function

* TODO.md: Remove 'imf' from list.
* docs/src/ref.md: Add "imf" entry to the manual.
* src/utils.jl: Add imf entry to "utils".
* src/imf.jl: Contains code of imf function.
* test/util-tests.jl: Include tests for imf.

Values now match (with the normal difference in least significant digits) with those
of IDL AstroLib.

* Implement recommended changes and modify test

* One more trick up the sleeve

Change sum() to dot() as only scalar product is involved

Author: TestSubjector

TODO.md

@@ -75,7 +75,6 @@ Missing in AstroLib.jl
 * `helio`
 * `hor2eq`
 * `imcontour`
-* `imf`
 * `planet_coords`
 * `qdcb_grid`
 * `tdb2tdt`

docs/src/ref.md

@@ -142,6 +142,7 @@ julia> AstroLib.planets["saturn"].mass
 [`deredd()`](@ref),
 [`flux2mag()`](@ref),
 [`gal_uvw()`](@ref),
+[`imf()`](@ref),
 [`ismeuv()`](@ref),
 [`kepler_solver()`](@ref),
 [`lsf_rotate()`](@ref),

src/imf.jl

@@ -0,0 +1,76 @@
+# This file is a part of AstroLib.jl. License is MIT "Expat".
+
+function _imf{T<:AbstractFloat}(mass::AbstractVector{T}, expon::AbstractVector{T},
+                                mass_range::AbstractVector{T})
+    ne_comp = length(expon)
+    if length(mass_range) != ne_comp + 1
+        error("Length of array mass_range is not one more than that of expon")
+    end
+    integ = Vector{T}(ne_comp)
+    joint = ones(T, ne_comp)
+    for i = 1:ne_comp
+        if expon[i] != -1
+            integ[i] = (mass_range[i+1]^(1 + expon[i]) - mass_range[i]^(1 + expon[i]))/
+                       (1 + expon[i])
+        else
+            integ[i] = log(mass_range[i+1]/mass_range[i])
+        end
+        if i != 1
+            joint[i] = joint[i-1]*(mass_range[i]^(expon[i-1] - expon[i]))
+        end
+    end
+    norm = joint./(dot(integ, joint))
+    psi = zeros(mass)
+    for i = 1:ne_comp
+        test = find(mass_range[i].< mass.<mass_range[i+1])
+        if length(test) !=0
+            psi[test] = norm[i].*(mass[test].^expon[i])
+        end
+    end
+    return psi
+end
+
+"""
+    imf(mass, expon, mass_range) -> psi
+
+### Purpose ###
+
+Compute an N-component power-law logarithmic initial mass function (IMF).
+
+### Explanation ###
+
+The function is normalized so that the total mass distribution equals
+one solar mass.
+
+### Arguments ###
+
+* `mass`: mass in units of solar mass, vector.
+* `expon`: power law exponent, vector. The number of values in expon equals
+  the number of different power-law components in the IMF.
+* `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 and positive.
+
+### Output ###
+
+* `psi`: mass function, number of stars per unit logarithmic mass interval
+  evaluated for supplied masses.
+
+### Example ###
+
+Show 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.
+
+```julia
+julia> imf([3], [-1.35], [0.1, 110]) / 3
+1-element Array{Float64,1}:
+ 0.0129414
+```
+
+### Notes ###
+
+Code of this function is based on IDL Astronomy User's Library.
+"""
+imf(mass::AbstractVector{<:Real}, expon::AbstractVector{<:Real}, mass_range::AbstractVector{<:Real}) =
+    _imf(float(mass), float(expon), float(mass_range))

src/utils.jl

@@ -70,6 +70,9 @@ export helio_jd
 include("helio_rv.jl")
 export helio_rv
 
+include("imf.jl")
+export imf
+
 include("ismeuv.jl")
 export ismeuv
 

test/utils-tests.jl

@@ -257,6 +257,15 @@ end
 @test helio_rv.([0.1, 0.9], 0, 1, 0, 100, 0.6, 45) ≈
     [-45.64994926111004, 89.7820347358485]
 
+# Test imf
+@testset "imf" begin
+    @test_throws ErrorException imf([5], [-6.75], [0.9])
+    @test imf([0.1, 0.01], [-0.6, -1], [ 0.007, 1.8, 110] ) ≈
+        [0.49627714725007616, 1.9757149090208912   ]
+    @test imf.([[3],[5]], [[-1.35], [-0.6, -1.7]], [[0.1, 100], [0.007, 1.8, 110]]) ≈
+        [0.038948937298846235, 0.027349915327755464]
+end
+
 # Test ismeuv
 @testset "ismeuv" begin
     @test ismeuv(304, 1e20) ≈ 58.30508020244554