/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * This program is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details.

 * You should have received a copy of the GNU General Public License along with
 * this program. If not, see <https://www.gnu.org/licenses/>.
 */

#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[2], ts[2], ts_sigma[2], mu_noise, mu_indirect, tmean;

    n = 1000;
    mu_noise = 0.1;
    mu_indirect = 0.5;
    mu[0] = 1.0;
    mu[1] = 1.0;
    ts[0] = 100.0;
    ts[1] = 120.0;
    tmean = 100.0;
    ts_sigma[0] = PMT_TTS;
    ts_sigma[1] = 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[0] = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"mu-shower")) {
                mu[1] = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"ts")) {
                ts[0] = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"ts-shower")) {
                ts[1] = 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")) {
                ts_sigma[0] = strtod(argv[++i],NULL);
                continue;
            } else if (!strcmp(argv[i]+2,"ts-sigma")) {
                ts_sigma[1] = 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,2,ts,tmean,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,2,ts,tmean,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;
}