#include "solid_angle.h"
#include <math.h>
#include <stdio.h>

struct solid_angle_results {
    double L;
    double r0;
    double omega;
} solid_angle_results[] = {
    {0.5,0.0,3.4732594},
    {0.5,0.2,3.4184435},
    {0.5,0.4,3.2435434},
    {0.5,0.6,2.9185178},
    {0.5,0.8,2.4122535},
    {0.5,1.0,1.7687239},
    {0.5,1.2,1.1661307},
    {0.5,1.4,0.7428889},
    {0.5,1.6,0.4841273},
    {0.5,1.8,0.3287007},
    {0.5,2.0,0.2324189}
};

int isclose(double a, double b, double rel_tol, double abs_tol)
{
    /* Returns 1 if a and b are "close". This algorithm is taken from Python's
     * math.isclose() function.
     *
     * See https://www.python.org/dev/peps/pep-0485/. */
    return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol);
}

int test_solid_angle(char *err)
{
    /* Tests the get_solid_angle() function. */
    int i;
    double pmt[3] = {0,0,0};
    double pos[3] = {0,0,1};
    double n[3] = {0,0,1};
    double r = 1.0;
    double solid_angle;

    solid_angle = get_solid_angle(pos,pmt,n,r);

    if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-9, 0)) {
        sprintf(err, "solid angle = %.2f, but expected %.2f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
        return 1;
    }

    for (i = 0; i < sizeof(solid_angle_results)/sizeof(struct solid_angle_results); i++) {
        pos[0] = solid_angle_results[i].r0*r;
        pos[2] = solid_angle_results[i].L*r;

        solid_angle = get_solid_angle(pos,pmt,n,r);

        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-4, 0)) {
            sprintf(err, "solid angle = %.2f, but expected %.2f", solid_angle, solid_angle_results[i].omega);
            return 1;
        }
    }

    return 0;
}

int main(int argc, char **argv)
{
    char err[256];

    if (!test_solid_angle(err)) {
        printf("[\033[92mok\033[0m] test_solid_angle\n");
        return 0;
    } else {
        printf("[\033[91mfail\033[0m] test_solid_angle: %s\n", err);
        return 1;
    }
}