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


#define MAXOBS 5000
#define MAXARG 4
#define MAXNb 20
#define MAXNc 5

int main(int argc, char *argv[]) 
{
    FILE *in;
    int N = 0, i, j;
    double ObsJD[MAXOBS],ObsMag[MAXOBS];
    double SumJD = 0,SumMag = 0;
    char fmt[255];
    int rising = 0;
    int compact = 1;
    double theta = 1.0;

    double X = 0, X2 = 0, sigma2;
    double phase, s2, f;
    int bin;

    int Nb = 5, Nc = 2, NbNc = Nb * Nc;
    double x[MAXNc][MAXNb], n[MAXNc][MAXNb];

    double f0 = atof(argv[1]), f1 = atof(argv[2]);
    double df = atof(argv[3]);

    if (argc == MAXARG)
    { // Use STDIN
      in = stdin;
    }
    else
    {
//fprintf(stderr,"File name = %s\n",argv[MAXARG]);
      in = fopen(argv[MAXARG],"r");
    }

    while ( !feof(in) ) 
    {
      if(fscanf(in,"%lf%lf",&(ObsJD[N]),&(ObsMag[N])))
      {
//printf("%d: %.4f %.3f\n", (N + 1), ObsJD[N],ObsMag[N]);
        if (ObsJD[N] && ObsMag[N])
        {
          SumJD += ObsJD[N];
          SumMag += ObsMag[N];
          ++N;
        }
      }
    }

    sprintf(fmt,"%%.%df %%.3f\n",abs((int)(log10(df))) + 2);

//fprintf(stderr,"%d %.4f %.3f %.3f -- %s\n", N, SumJD/N, SumMag/N, df, fmt);
//printf(fmt, f0, df);
//printf(fmt, f1, df);

    for (i = 0; i < N; i++)
    {
      ObsJD[i] -= SumJD/N;
      ObsMag[i] -= SumMag/N;
      X += ObsMag[i];
      X2 += ObsMag[i] * ObsMag[i];
    }

    sigma2 = (X2 - X * X/N)/(N - 1);

//printf("%.6f %.6f %.6f\n\n", X, X2, sigma2);

    for (f = f0; f <= f1; f += df)
    {
      // For each frequency
      // Initialize arrays
      for (j = 0; j < Nc; j++) 
      {
        for (i = 0; i < Nb; i++) 
        {
          x[j][i] = 0;
          n[j][i] = 0;
        }
      }

      // Loop through all observations
      for (i = 1; i <= N; i++) 
      {
//printf("\ni%d-%.5f ",i,df);
        for (j = 0; j < Nc; j++) 
        {
          phase = ObsJD[i] * f + j/NbNc;
          bin = (int)(Nb * (phase - floor(phase)));
//printf("j%d-%.5f:%d %.3f | ",j,df,bin,phase);
          x[j][bin] += ObsMag[i];
          n[j][bin] += 1;
        }
      }

      s2 = 0;
      for (j = 0; j < Nc; j++) 
      {
        for (i = 0; i < Nb; i++) 
        {
          if (n[j][i] > 0) 
          {
            s2 += x[j][i] * x[j][i] / n[j][i];
          }
        }
      }
      s2 = (Nc * X2 - s2)/Nc/(N - Nb);

      if (compact)
      { // Show only extrema of the statistic
        if (s2 < theta)
        { // statistic descends
          if (rising > 0)
          {// statistic was rising, so has reached minimum
            printf(fmt, f - df, theta/sigma2);
          }
          rising = -1;
        }
        else if (s2 > theta)
        {// statistic rises
          if (rising < 0)
          {// statistic was descending, so has reached maximum
            printf(fmt, f - df, theta/sigma2);
          }
          rising = 1;
        }
        else
          rising = 0;

        theta = s2;
      }
      else
        printf(fmt, f, s2/sigma2);

      // End for each frequency
    }

    if (compact)
    {
      printf(fmt, f - df, theta/sigma2);
    }

    return EXIT_SUCCESS;
}