#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <errno.h> /* for errno */
#include <string.h> /* for strerror() */
#include "likelihood.h"

void usage(void)
{
    fprintf(stderr,"Usage: ./test-time-pdf\n");
    fprintf(stderr,"  -n            number of points\n");
    fprintf(stderr,"  --mu-noise    Expected number of noise PE\n");
    fprintf(stderr,"  --mu-indirect Expected number of PE from scattered/reflected light\n");
    fprintf(stderr,"  --mu-shower   Expected number of PE from the electromagnetic shower\n");
    fprintf(stderr,"  --ts          Expected time of direct light\n");
    fprintf(stderr,"  --ts-shower   Expected time of shower light\n");
    fprintf(stderr,"  --t-mean      Starting time for scattered/reflected light\n");
    fprintf(stderr,"  --t-sigma     Standard deviation of direct light\n");
    fprintf(stderr,"  --ts-sigma    Standard deviation of shower light\n");
    fprintf(stderr,"  -h            Display this help message\n");
    exit(1);
}

int main(int argc, char **argv)
{
    size_t i, n;
    double *x;
    double mu_noise, mu_indirect, mu_direct, mu_shower, ts, ts_shower, tmean, tsigma, ts_sigma;

    n = 1000;
    mu_noise = 0.1;
    mu_indirect = 0.5;
    mu_direct = 1.0;
    mu_shower = 1.0;
    ts = 100.0;
    ts_shower = 120.0;
    tmean = 100.0;
    tsigma = PMT_TTS;
    ts_sigma = 10.0;

    for (i = 1; i < argc; i++) {
        if (!strncmp(argv[i], "--", 2)) {
            if (!strcmp(argv[i]+2,"mu-noise")) {
                mu_noise = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"mu-indirect")) {
                mu_indirect = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"mu-direct")) {
                mu_direct = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"mu-shower")) {
                mu_shower = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"ts")) {
                ts = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"ts-shower")) {
                ts_shower = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"t-mean")) {
                tmean = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"t-sigma")) {
                tsigma = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"ts-sigma")) {
                ts_sigma = strtod(argv[++i],NULL);
                continue;
            }
        } else if (argv[i][0] == '-') {
            switch (argv[i][1]) {
            case 'n':
                n = atoi(argv[++i]);
                break;
            case 'h':
                usage();
            default:
                fprintf(stderr, "unrecognized option '%s'\n", argv[i]);
                exit(1);
            }
        }
    }

    x = malloc(n*sizeof(double));

    for (i = 0; i < n; i++) {
        x[i] = GTVALID*i/(n-1);
    }

    FILE *pipe = popen("graph -T X --bitmap-size 2000x2000 -X 'Time (ns)' -Y Probability", "w");

    if (!pipe) {
        fprintf(stderr, "error running graph command: %s\n", strerror(errno));
        exit(1);
    }

    for (i = 0; i < n; i++) {
        fprintf(pipe, "%.10f %.10f\n", x[i], time_pdf(x[i],mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma));
    }
    fprintf(pipe, "\n\n");

    for (i = 0; i < n; i++) {
        fprintf(pipe, "%.10f %.10f\n", x[i], time_cdf(x[i],mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma));
    }
    fprintf(pipe, "\n\n");

    if (pclose(pipe)) {
        fprintf(stderr, "error closing graph command: %s\n", strerror(errno));
        exit(1);
    }

    free(x);

    return 0;
}