/* pwr.c		Joe Mietus		Oct 7 2008 */

/*
pwr :
Usage: pwr [-r] [-L] lo1 hi1 ... [lo10 hi10]
  Calculate total (and relative) power in fft between lo and hi.
  options :
  [-r] : print ratio of powers to total
  [-L] : print powers on one line
*/


#include <stdio.h>
#include <stdlib.h>

#define MAXBANDS 10

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

    int i, n, nbands, rflag, lflag;
    double lo[MAXBANDS], hi[MAXBANDS], pr[MAXBANDS], tot;
    double freq[2], mag[2], hbin[2], pwr;


    if (argc < 3) {
        usage(argv[0]);
        exit(1);
    }

    rflag = lflag = 0;
    for (i=1, argc--; argv[i][0] == '-'; i++, argc--) {
        switch(argv[i][1]) {
        case 'r': rflag = 1;
                  break;
        case 'L': lflag = 1;
                  break;
        default: usage(argv[0]);
                 exit(1);
        }
    }
        
    if (argc % 2 != 0) {
        usage(argv[0]);
        exit(1);
    }

    if ((nbands = argc/2) > 10) {
        fprintf(stderr, "%s : max power bands = MAXBANDS\n");
        exit(1);
    }

    for (n=0; n<nbands; n++) {
        lo[n] = atof(argv[i++]);
        hi[n] = atof(argv[i++]);
        pr[n] = 0.0;
    }

    tot = 0.0;

    if (scanf("%lf%lf", &freq[0], &mag[0]) != 2)
        exit(2);
    if (scanf("%lf%lf", &freq[1], &mag[1]) != 2)
        exit(2);

    pwr = mag[0]*mag[0];
    tot += pwr;
    for (n=0; n<nbands; n++) {
        if (freq[0] >= lo[n] && freq[0] <= hi[n])
            pr[n] += pwr;
    }
    freq[0] = freq[1];
    mag[0] = mag[1];

    while (scanf("%lf%lf", &freq[1], &mag[1]) == 2) {
        pwr = mag[0]*mag[0];
        tot += pwr;
        for (n=0; n<nbands; n++) {
            if (freq[0] >= lo[n] && freq[0] <= hi[n])
                pr[n] += pwr;
        }
        freq[0] = freq[1];
	mag[0] = mag[1];
    }

    pwr = mag[0]*mag[0];
    tot += pwr;
    for (n=0; n<nbands; n++) {
        if (freq[0] >= lo[n] && freq[0] <= hi[n])
        pr[n] += pwr;
    }

    if (lflag) {
        printf("%g ", tot);
        for (n=0; n<nbands; n++) {
            printf("%g ", pr[n]);
            if (rflag)
                printf("(%.3f) ", pr[n]/tot);
        }
        printf("\n");
    }
    else {
        printf("Total = %g\n", tot);
        for (n=0; n<nbands; n++) {
            printf("%g - %g = %g", lo[n], hi[n], pr[n]);
            if (rflag)
                printf(" (%.3f)\n", pr[n]/tot);
            else
                printf("\n");
        }
    }
}


usage(prog)
char *prog;
{
    fprintf(stderr, "Usage: %s [-r] [-L] lo1 hi1 ... [lo10 hi10]\n", prog);
    fprintf(stderr, " Calculate total (and relative) power");
    fprintf(stderr, " in fft between lo and hi.\n");
    fprintf(stderr, " options :\n");
    fprintf(stderr, " [-r] : print ratio of powers to total\n");
    fprintf(stderr, " [-L] : print powers on one line\n");

}