add a fast likelihood function

This commit adds a fast function to calculate the expected number of PE at a
PMT without numerically integrating over the track. This calculation is *much*
faster than integrating over the track (~30 ms compared to several seconds) and
so we use it during the "quick" minimization phase of the fit to quickly find
the best position.

1 files changed, 120 insertions, 0 deletions

diff --git a/src/misc.c b/src/misc.c
index bff0ba9..82d114b 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -1,6 +1,126 @@
 #include "misc.h"
 #include <math.h>
 #include <stdlib.h> /* for size_t */
+#include <gsl/gsl_sf_gamma.h>
+
+static struct {
+    int n;
+    double f;
+} ln_fact_table[LNFACT_MAX + 1] = {
+    {0,0},
+    {1,0},
+    {2,0.69314718055994529},
+    {3,1.791759469228055},
+    {4,3.1780538303479458},
+    {5,4.7874917427820458},
+    {6,6.5792512120101012},
+    {7,8.5251613610654147},
+    {8,10.604602902745251},
+    {9,12.801827480081469},
+    {10,15.104412573075516},
+    {11,17.502307845873887},
+    {12,19.987214495661885},
+    {13,22.552163853123425},
+    {14,25.19122118273868},
+    {15,27.89927138384089},
+    {16,30.671860106080672},
+    {17,33.505073450136891},
+    {18,36.395445208033053},
+    {19,39.339884187199495},
+    {20,42.335616460753485},
+    {21,45.380138898476908},
+    {22,48.471181351835227},
+    {23,51.606675567764377},
+    {24,54.784729398112319},
+    {25,58.003605222980518},
+    {26,61.261701761002001},
+    {27,64.557538627006338},
+    {28,67.88974313718154},
+    {29,71.257038967168015},
+    {30,74.658236348830158},
+    {31,78.092223553315307},
+    {32,81.557959456115043},
+    {33,85.054467017581516},
+    {34,88.580827542197682},
+    {35,92.136175603687093},
+    {36,95.719694542143202},
+    {37,99.330612454787428},
+    {38,102.96819861451381},
+    {39,106.63176026064346},
+    {40,110.32063971475739},
+    {41,114.03421178146171},
+    {42,117.77188139974507},
+    {43,121.53308151543864},
+    {44,125.3172711493569},
+    {45,129.12393363912722},
+    {46,132.95257503561632},
+    {47,136.80272263732635},
+    {48,140.67392364823425},
+    {49,144.5657439463449},
+    {50,148.47776695177302},
+    {51,152.40959258449735},
+    {52,156.3608363030788},
+    {53,160.3311282166309},
+    {54,164.32011226319517},
+    {55,168.32744544842765},
+    {56,172.35279713916279},
+    {57,176.39584840699735},
+    {58,180.45629141754378},
+    {59,184.53382886144948},
+    {60,188.6281734236716},
+    {61,192.7390472878449},
+    {62,196.86618167289001},
+    {63,201.00931639928152},
+    {64,205.1681994826412},
+    {65,209.34258675253685},
+    {66,213.53224149456327},
+    {67,217.73693411395422},
+    {68,221.95644181913033},
+    {69,226.1905483237276},
+    {70,230.43904356577696},
+    {71,234.70172344281826},
+    {72,238.97838956183432},
+    {73,243.26884900298271},
+    {74,247.57291409618688},
+    {75,251.89040220972319},
+    {76,256.22113555000954},
+    {77,260.56494097186322},
+    {78,264.92164979855278},
+    {79,269.29109765101981},
+    {80,273.67312428569369},
+    {81,278.06757344036612},
+    {82,282.4742926876304},
+    {83,286.89313329542699},
+    {84,291.32395009427029},
+    {85,295.76660135076065},
+    {86,300.22094864701415},
+    {87,304.68685676566872},
+    {88,309.1641935801469},
+    {89,313.65282994987905},
+    {90,318.1526396202093},
+    {91,322.66349912672615},
+    {92,327.1852877037752},
+    {93,331.71788719692847},
+    {94,336.26118197919845},
+    {95,340.81505887079902},
+    {96,345.37940706226686},
+    {97,349.95411804077025},
+    {98,354.53908551944079},
+    {99,359.1342053695754},
+    {100,363.73937555556347},
+};
+
+double lnfact(unsigned int n)
+{
+    /* Returns the logarithm of n!.
+     *
+     * Uses a lookup table to return results for n < 100. */
+
+    if (n <= LNFACT_MAX)
+        return ln_fact_table[n].f;
+
+    return gsl_sf_lnfact(n);
+}
 
 double kahan_sum(double *x, size_t n)
 {