/** 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;
  FILE *ft;
  char f1[20], f2[20];
  char line[120];
  int  i, imax, index, na, nb, nc, nd1, nd2, nt, nw;
  int  nl, ns;
  int  iprt, step;
  int  j;
  int  nin=0, nout=0;
  int  len;
  int  nzone[1000000], nbin[1000000]; /* annulus and mass bin for n-bodies or tracers */
  double time=0., time0=-1.;
  double data[10000000];
  double amax=0.;
  double yr=3.1556926e7;
  double frac=0., integer=0., now=0.;

  step = 10;
  if (argc < 4) {
     printf("No input file or step information\n");
     printf("usage: extract BINFILE n1 n2 (with n2 < n1)\n");
     printf("This code prints physical parameters for n-bodies and tracers\n");
     exit(0);
  } else {
/* open file */
     printf("Using INPUT FILE: %s\n", *++argv);
     strcpy(f2, *argv);
     if ((fp = fopen(*argv, "r")) == NULL) {
        printf("Can't open file %s\n", *argv);
        exit(1);
     } else {
        printf("input timeindex is %s\n", *++argv);
        imax = atoi(*argv) ;
        printf("input step is %s\n", *++argv);
        step = atoi(*argv) ;
        for (i = 0; i < imax; i++) {
           now = ((double) i) / ((double) step) ;
           frac = modf( now, &integer) ;
           if (frac == 0. || i == imax) {
              iprt = 0 ;
              sprintf(f1, "%s.%d\0", f2, i) ;
              ft = fopen(f1,"w") ;
           } else {
             iprt = -1 ;
           }
           nb = 0;
           nc = 0;
           nd1 = 0;
           nd2 = 0;
           nt = 0;
           nw = 0;
           nin = 0;
           nout = 0;
           len = F_GetLine(line, sizeof(line), fp);
           if (len <= 0) {
              exit(1);
           } else {
              sscanf(line, "%le %d %d %d %d %d %d %n", &time, &index, &ns, &nb, &nd1, &nt, &nd2, &nc);
              fprintf(ft, "# Mass in grams, distances in cm, velocities in cm/sec\n");
              fprintf(ft, "# time:  %e sec %e yr\n", time, time/yr);
              if (iprt == 0) 
                 printf("%d %e %s\n", i, time, f1);
/* 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);
                 na = 0;
	         nl = 0;
                 for (j = 0; j < nd1; j+=16) {
/* if data[j] non-zero and data[j+5] < amax1, n-body is active; otherwise not active 
 * (data array is all zeros OR semimajor axis has exceeded maximum) 
*/
                    if (data[j] != 0.)
                       na +=1 ;
                    if (data[j] < 0.)
                       nl +=1 ;
                 }
                 if (iprt == 0) {
                    fprintf(ft, "# Number of active n-bodies: %d number lost: %d\n", na, nl);
                    fprintf(ft, "#    Mass        Number        Radius     SemiMajorAxis        X             Y             Z           XDOT           YDOT          ZDOT            ECC          INCL\n");
                    for (j = 0; j < nd1; j+=16) {
                       if (data[j] != 0.)
                          fprintf(ft, "%13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e\n", data[j], data[j+1], data[j+3], data[j+5], data[j+6], data[j+7], data[j+8], data[j+9], data[j+10], data[j+11], data[j+13], data[j+14]);
                    }
                    fprintf(ft, "# Number of active tracers: %d\n", nt);
                    fprintf(ft, "#    Mass        Number        Radius     SemiMajorAxis        X             Y             Z           XDOT           YDOT          ZDOT            ECC          INCL\n");

                 }
/* 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);
                    for (j = 0; j < nd2; j+=16) {
                       if (iprt == 0) {
                          fprintf(ft, "%13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e %13.7e\n", data[j], data[j+1], data[j+3], data[j+5], data[j+6], data[j+7], data[j+8], data[j+9], data[j+10], data[j+11], data[j+13], data[j+14]);
                       }
                    }
                 }
              } else {
                 exit(1);
              }
           }
           if (iprt == 0)
              fclose(ft);
        }
     }
  }
  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 */