/** summary.c
 * cycle through Orchestra binary file and print out some useful stats
 * need to set amax = maximum semimajor axis for active objects
 * **************************************************************************/
/* organization of data array
0 = total mass in grams (m)
1 = mass in gas 
2 = mass * rho where rho is mass density in grams**2 / cm**3
3 = radius in cm
4 = luminosity
5 = semimajor axis = a
6, 7, 8 = x, y, z
9, 10, 11 = vx, vy, vz
12 = da/dt or a
13 = ecc
14 = incl
15 = dm/dt or m
*/

/* usage

  gcc -lm summary.c -o summary
  summary filename

*/

#include <stdio.h>      
#include <stdarg.h>
#include <math.h>
#include <time.h>
#include <signal.h>
#include <unistd.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>

int F_GetLine(char *line, int max, FILE *fp);

int main(int argc, char *argv[])
{
  FILE *fp;
  char f1[20];
  char line[120];
  int  i, imax, index, na, nb, nc, nd1, nd2, nt, nw;
  int  il, ns;
  int  j;
  int  len;
  int  nzone[1000000], nbin[1000000]; /* annulus and mass bin for n-bodies or tracers */
  double time=0.;
  double data[10000000];
  double amin, amax;
  double re=6.3781e8;
  double rp=1.1883e8;
  double mp=1.303e25;
  double MAct=0., MInAct=0.;
	            double MSum = 0.;

  imax = 1000000;
  if (argc == 1) {
     printf("No input file\n");
     printf("usage: %s BINFILE\n");
     printf("This code prints summary informtion for records in binary file BINFILE\n");
     exit(0);
  } else {
/* open file */
     printf("#INPUT FILE: %s\n", *++argv);
     sprintf(f1, "%s\0\n", *argv);
     if ((fp = fopen(*argv, "rb")) == NULL) {
        printf("Can't open file %s\n", *argv);
        exit(1);
     } else {
        for (i = 0; i < imax; i++) {
           nb = 0;
           nc = 0;
           nd1 = 0;
           nd2 = 0;
           nt = 0;
           nw = 0;
           len = F_GetLine(line, sizeof(line), fp);
           if (len <= 0) {
              if (il == 0) {
                 printf("# Record    Time      TimeIndex   NAct      AMin      AMax       MAct       MInAct\n");
                 printf("# Time in units of seconds\n");
                 printf("# Semimajor axis of active tracers (AMin or AMax) in units of R-Pluto\n");
                 printf("# Mass of active (MAct) or inactive (MInAct) tracers in units of grams\n");
	         il = 1;
	      }
           } else {
              sscanf(line, "%le %d %d %d %d %d %d %n", &time, &index, &ns, &nb, &nd1, &nt, &nd2, &nc);
/* read n-body info */
              if (nd1 > 0) {
                 nw = fread((char*)nzone,sizeof(int),nb,fp);
                 nw = fread((char*)nbin,sizeof(int),nb,fp);
                 nw = fread((char*)data,sizeof(double),nd1,fp);
/* read tracer info */
                 if (nt > 0) {
                    nw = fread((char*)nzone,sizeof(int),nt,fp);
                    nw = fread((char*)nbin,sizeof(int),nt,fp);
                    nw = fread((char*)data,sizeof(double),nd2,fp);
	            na = 0;
                    amin = 1.e15;
                    amax = 0.;
	            MAct = 0.;
	            MInAct = 0.;
                    for (j = 0; j < nd2; j+=16) {
/* number of active tracers */
	               if (data[j] >= 0.)
	                  MAct += data[j];
	               else
                          MInAct -= data[j] ;
                       if (data[j+5] > 3.07*re && data[j+5] < 1.e3*re) {
                          na += 1;
/* find max and min semimajor axis for tracers */
                          if (data[j+5] > amax)
                             amax = data[j+5] ;
                          if (data[j+5] > 0. && data[j+5] < amin)
                             amin = data[j+5] ;
                       }
                    }
                    amin /= rp;
                    amax /= rp;
                 }
/* scale min and max values by semimajor axis of PC binary */
                 if (i == 0) {
                    printf("# Timestep summary\n");
                    printf("# Time in units of seconds\n");
                    printf("# Semimajor axis of active tracers (AMin or AMax) in units of R-Pluto\n");
                    printf("# Mass of active (MAct) or inactive (MInAct) tracers in units of grams\n");
                    printf("# Record    Time      TimeIndex   NAct      AMin      AMax       MAct       MInAct\n");
	         }
                 printf("%6d %12.4e %9d %9d %8.2f %10.2f %11.2e %11.2e\n", i, time, index, na, amin, amax, MAct, MInAct);
              } else {
                    printf("# Record    Time      TimeIndex   NAct      AMin      AMax       MAct       MInAct\n");
                    printf("# Time in units of seconds\n");
                    printf("# Semimajor axis of active tracers (AMin or AMax) in units of R-Pluto\n");
                    printf("# Mass of active (MAct) or inactive (MInAct) tracers in units of grams\n");
                 exit(1);
              }
           }
        }
     }
  }
  fclose(fp);
  return 0;
}

/* F_GetLine
 *    Input:  string ptr, line length, file ptr
 *       Output: string
 *          Used to read lines from an input parameter file.
 *          */
int
F_GetLine(char *line, int max, FILE *fp)
{
  if (fgets(line, max, fp) == NULL)
    return 0;
  else
    return strlen(line);
} /* end F_GetLine */