/*+
    VIDE -- Void IDentification and Examination -- ./c_tools/stacking/pruneVoids.cpp
    Copyright (C) 2010-2014 Guilhem Lavaux
    Copyright (C) 2011-2014 P. M. Sutter

    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; version 2 of the License.

    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, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+*/



// Reads in the void catalog and removes any void that could potentially
//   be affected by a mock particle. It does this by computing the longest
//   particle distance within each void and comparing it to the distance
//   of the nearest mock particle. If the void could potentially by rotated
//   to include this particle, we throw out the void.

// This is placed here instead of using the edgeAvoidance option in 
//   stackVoidsZero so that we can optionally filter the entire
//   catalog at once before the stacking phase. This is useful
//   for producing a "clean" void catalog for other purposes.

#include "gsl/gsl_math.h"
#include "gsl/gsl_eigen.h"
#include "string.h"
#include "ctype.h"
#include "stdlib.h"
#include <math.h>
#include <stdio.h>
#include <netcdf>
#include "pruneVoids_conf.h"
#include <vector>
#include "assert.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>

#define LIGHT_SPEED 299792.458
#define MPC2Z 100./LIGHT_SPEED
#define Z2MPC LIGHT_SPEED/100.

#define CENTRAL_VOID 1
#define EDGE_VOID 2

using namespace std;

typedef struct partStruct {
  float x, y, z, vol;
  int nadj, ncnt;
  int *adj;
  int edgeFlag;
} PART;

typedef struct zoneStruct {
  int numPart;
  int *partIDs;
} ZONE2PART;

typedef struct voidZoneStruct {
  int numZones;
  int *zoneIDs;
} VOID2ZONE;

typedef struct voidStruct {
  float vol, coreDens, zoneVol, densCon, voidProb, radius;
  float rescaledCoreDens;
  int voidID, numPart, numZones, coreParticle, zoneNumPart;
  float maxRadius, nearestFlag, centralDen, redshift, redshiftInMpc;
  float nearestEdge;
  float center[3], macrocenter[3];
  int accepted;
  int voidType;
  int parentID, numChildren, level;
  bool isLeaf, hasHighCentralDen;
  gsl_vector *eval;
  gsl_matrix *evec;
  float ellip;
} VOID;

// this defines the expansion function that we will integrate
// Laveaux & Wandelt (2012) Eq. 24
struct my_expan_params { double Om; double w0; double wa; };
double expanFun (double z, void * p) {
  struct my_expan_params * params = (struct my_expan_params *)p;
  double Om = (params->Om);
  double w0 = (params->w0);
  double wa = (params->wa);

  //const double h0 = 1.0;
  const double h0 = 0.71;
  double ez;

  double wz = w0 + wa*z/(1+z);

  ez = Om*pow(1+z,3) + (1.-Om);
  //ez = Om*pow(1+z,3) + pow(h0,2) * (1.-Om)*pow(1+z,3+3*wz);

  ez = sqrt(ez);
  //ez = sqrt(ez)/h0;

  ez = 1./ez;

  return  ez;
}

void openFiles(string outputDir, string sampleName, 
               int numPartTot, int numKept,
               FILE** fpOutput);

void closeFiles(FILE* fpOutput);

void outputVoids(string outputDir, string sampleName, int numPartTot, 
                vector<VOID> voids, 
                bool isObservation, double *boxLen);

int main(int argc, char **argv) {

  // initialize arguments
  pruneVoids_info args;
  pruneVoids_conf_params args_params;

  pruneVoids_conf_init(&args);
  pruneVoids_conf_params_init(&args_params);

  args_params.check_required = 0;
  if (pruneVoids_conf_ext (argc, argv, &args, &args_params))
    return 1;

  if (!args.configFile_given) {
    if (pruneVoids_conf_required (&args, PRUNEVOIDS_CONF_PACKAGE))
      return 1;
  } else {
    args_params.check_required = 1;
    args_params.initialize = 0;
    if (pruneVoids_conf_config_file (args.configFile_arg,
                                              &args,
                                              &args_params))
    return 1;
  }

  // initialize cosmology integrator and interpolator
  gsl_function expanF;
  expanF.function = &expanFun;
  struct my_expan_params expanParams;
  expanParams.Om = args.omegaM_arg;
  expanParams.w0 = -1.0;
  expanParams.wa = 0.0;
  expanF.params = &expanParams;
  double result, error;
  size_t nEval;

  int iZ, numZ = 8000;
  double maxZ = 10.0, z, *dL, *redshifts;
  dL = (double *) malloc(numZ * sizeof(double));  
  redshifts = (double *) malloc(numZ * sizeof(double));
  for (iZ = 0; iZ < numZ; iZ++) {
    z = iZ * maxZ/numZ;
    gsl_integration_qng(&expanF, 0.0, z, 1.e-6, 1.e-6, &result, &error, &nEval);
    dL[iZ] = result*LIGHT_SPEED/100.;
    //printf("HERE %e %e\n", z, dL[iZ]);
    redshifts[iZ] = z;
  }
  gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, numZ);
  gsl_interp_init(interp, dL, redshifts, numZ);
  gsl_interp_accel *acc = gsl_interp_accel_alloc();

 
  int i, p, p2, numPartTot, numZonesTot, dummy, iVoid;
  int numVoids;
  FILE *fp;
  PART *part, *voidPart, *flaggedPart;
  ZONE2PART *zones2Parts;
  VOID2ZONE *void2Zones;
  vector<VOID> voids;
  float *temp, junk, voidVol;
  int junkInt, voidID, numPart, numZones, zoneID, partID, maxNumPart;
  int coreParticle, zoneNumPart;
  float coreDens, zoneVol, densCon, voidProb, dist[3], dist2, minDist, maxDist;
  float centralRad, centralDen;
  double nearestEdge, redshift;
  char line[500], junkStr[10];
  string outputDir, sampleName, dataPortion, prefix;
  double ranges[3][2], boxLen[3], mul; 
  double volNormZobov, volNormObs, radius;
  int clock1, clock2, clock3, clock4;
  double interval;
  int periodicX=0, periodicY=0, periodicZ=0;
  string dataPortions[2];

  gsl_eigen_symmv_workspace *eigw = gsl_eigen_symmv_alloc(3);

  numVoids = args.numVoids_arg;

  clock1 = clock();

  // check for periodic box
  periodicX = 0;
  periodicY = 0;
  periodicZ = 0;
  if (!args.isObservation_flag) {
    if ( strchr(args.periodic_arg, 'x') != NULL) {
      periodicX = 1;
      printf("Will assume x-direction is periodic.\n");
    }
    if ( strchr(args.periodic_arg, 'y') != NULL) {
      periodicY = 1;
      printf("Will assume y-direction is periodic.\n");
    }
    if ( strchr(args.periodic_arg, 'z') != NULL) {
      periodicZ = 1;
      printf("Will assume z-direction is periodic.\n");
    }
  }

  // load box size
  printf("\n Getting info...\n");
  netCDF::NcFile f_info(args.extraInfo_arg, netCDF::NcFile::read);
  f_info.getAtt("range_x_min").getValues(&ranges[0][0]);
  f_info.getAtt("range_x_max").getValues(&ranges[0][1]);
  f_info.getAtt("range_y_min").getValues(&ranges[1][0]);
  f_info.getAtt("range_y_max").getValues(&ranges[1][1]);
  f_info.getAtt("range_z_min").getValues(&ranges[2][0]);
  f_info.getAtt("range_z_max").getValues(&ranges[2][1]);

  printf(" Range xmin %e\n", ranges[0][0]);
  printf(" Range xmax %e\n", ranges[0][1]);
  printf(" Range ymin %e\n", ranges[1][0]);
  printf(" Range ymax %e\n", ranges[1][1]);
  printf(" Range zmin %e\n", ranges[2][0]);
  printf(" Range zmax %e\n", ranges[2][1]);

  boxLen[0] = ranges[0][1] - ranges[0][0];
  boxLen[1] = ranges[1][1] - ranges[1][0];
  boxLen[2] = ranges[2][1] - ranges[2][0];
 
  // read in all particle positions
  clock3 = clock();
  printf("\n Loading particles...\n");
  fp = fopen(args.partFile_arg, "r");
  fread(&dummy, 1, 4, fp); 
  fread(&numPartTot, 1, 4, fp);
  fread(&dummy, 1, 4, fp);

  part = (PART *) malloc(numPartTot * sizeof(PART));
  temp = (float *) malloc(numPartTot * sizeof(float));

  volNormZobov = args.volNormZobov_arg;
  volNormObs = args.volNormObs_arg;

  fread(&dummy, 1, 4, fp);
  fread(temp, numPartTot, 4, fp);
  mul = ranges[0][1] - ranges[0][0];
  for (p = 0; p < numPartTot; p++)
    part[p].x = mul*temp[p];
  fread(&dummy, 1, 4, fp);
  fread(&dummy, 1, 4, fp);
  fread(temp, numPartTot, 4, fp);
  mul = ranges[1][1] - ranges[1][0];
  for (p = 0; p < numPartTot; p++) 
    part[p].y = mul*temp[p];
  fread(&dummy, 1, 4, fp);
  fread(&dummy, 1, 4, fp);
  fread(temp, numPartTot, 4, fp);
  mul = ranges[2][1] - ranges[2][0];
  for (p = 0; p < numPartTot; p++) 
    part[p].z = mul*temp[p];

  if (!args.isObservation_flag) {
    for (p = 0; p < numPartTot; p++)  {
      part[p].x += ranges[0][0];
      part[p].y += ranges[1][0];
      part[p].z += ranges[2][0];
    }
  }
  fclose(fp); 

  clock4 = clock();
  interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
  printf(" Read %d particles (%.2f sec)...\n", numPartTot, interval);
  
  // read in desired voids
  clock3 = clock();
  printf(" Loading voids...\n");
  fp = fopen(args.voidDesc_arg ,"r");
  fgets(line, sizeof(line), fp);
  sscanf(line, "%d %s %d %s", &junkInt, junkStr, &junkInt, junkStr);
  fgets(line, sizeof(line), fp);

  voids.resize(numVoids);
  i = 0;
  while (fgets(line, sizeof(line), fp) != NULL) {
    sscanf(line, "%d %d %d %f %f %d %d %f %d %f %f\n", &iVoid, &voidID, 
           &coreParticle, &coreDens, &zoneVol, &zoneNumPart, &numZones, 
           &voidVol, &numPart, &densCon, &voidProb);

    i++;
    voids[i-1].coreParticle = coreParticle;
    voids[i-1].zoneNumPart = zoneNumPart;
    voids[i-1].coreDens = coreDens;
    voids[i-1].zoneVol = zoneVol;
    voids[i-1].voidID = voidID;
    voids[i-1].vol = voidVol;
    voids[i-1].numPart = numPart;
    voids[i-1].numZones = numZones;
    voids[i-1].densCon = densCon;
    voids[i-1].voidProb = voidProb;

    voids[i-1].radius = pow(voidVol/volNormZobov*3./4./M_PI, 1./3.);
    voids[i-1].accepted = 1;

    voids[i-1].isLeaf = true;
    voids[i-1].hasHighCentralDen = false;
    voids[i-1].parentID = -1;  
    voids[i-1].numChildren = 0;
    voids[i-1].level = 0;  

    voids[i-1].eval = gsl_vector_alloc(3);
    voids[i-1].evec = gsl_matrix_alloc(3,3);
    voids[i-1].ellip = 0;
  }
  fclose(fp);

  // load up the zone membership for each void
  printf(" Loading void-zone membership info...\n");
  fp = fopen(args.void2Zone_arg, "r");
  fread(&numZonesTot, 1, 4, fp);
  void2Zones = (VOID2ZONE *) malloc(numZonesTot * sizeof(VOID2ZONE));

  for (iZ = 0; iZ < numZonesTot; iZ++) {
    fread(&numZones, 1, 4, fp);

    void2Zones[iZ].numZones = numZones;
    void2Zones[iZ].zoneIDs = (int *) malloc(numZones * sizeof(int));

    for (p = 0; p < numZones; p++) {
      fread(&void2Zones[iZ].zoneIDs[p], 1, 4, fp);
    }
  }
  fclose(fp);

  // now the particles-zone
  printf(" Loading zone-particle membership info...\n");
  fp = fopen(args.zone2Part_arg, "r");
  fread(&dummy, 1, 4, fp);
  fread(&numZonesTot, 1, 4, fp);
 
  zones2Parts = (ZONE2PART *) malloc(numZonesTot * sizeof(ZONE2PART));
  for (iZ = 0; iZ < numZonesTot; iZ++) {
    fread(&numPart, 1, 4, fp);
    
    zones2Parts[iZ].numPart = numPart;
    zones2Parts[iZ].partIDs = (int *) malloc(numPart * sizeof(int));

    for (p = 0; p < numPart; p++) {
      fread(&zones2Parts[iZ].partIDs[p], 1, 4, fp);
    }
  }

  // and volumes
  printf(" Loading particle volumes...\n");
  fp = fopen(args.partVol_arg, "r");
  fread(&numPartTot, 1, 4, fp);
  for (p = 0; p < numPartTot; p++) {
    fread(&temp[0], 1, 4, fp);
    part[p].vol = temp[0];
  }
  fclose(fp);

  // and finally edge flag info
  printf(" Loading particle edge flags...\n");
  int numFlagged = 0;
  fp = fopen(args.partEdge_arg, "r");
  for (p = 0; p < numPartTot; p++) {
    fscanf(fp, "%d", &part[p].edgeFlag);
    if (part[p].edgeFlag > 0) numFlagged++;
  }
  fclose(fp);

  // copy flagged galaxies to separate array for faster processing later
  flaggedPart = (PART *) malloc(numFlagged * sizeof(PART));
  int iFlag = 0;
  for (p = 0; p < numPartTot; p++) {
    if (part[p].edgeFlag > 0) {
      flaggedPart[iFlag] = part[p];
      iFlag++;
    }
  }

/* this was used for testing at one point
  // and finally finally adjacencies
  printf(" Loading particle adjacencies...\n");
  fp = fopen(args.partAdj_arg, "r");
  fread(&numPartTot, 1, 4, fp);
  if (numPartTot != numPartTot) {
    printf("NON-MATCHING MOCK INDICES!? %d %d\n", numPartTot, numPartTot);
    exit(-1);
  } 
  int tempInt;
  for (p = 0; p < numPartTot; p++) {
    fread(&tempInt, 1, 4, fp);
    part[p].nadj = tempInt;
    part[p].ncnt = 0;
    if (part[p].nadj > 0)
      part[p].adj = (int *) malloc(part[p].nadj * sizeof(int));
  }
  for (p = 0; p < numPartTot; p++) {
    fread(&tempInt, 1, 4, fp);
    int nin = tempInt;
    if (nin > 0) {
      for (int nAdj = 0; nAdj < nin; nAdj++) {
        int tempAdj;
        fread(&tempAdj, 1, 4, fp);

        // this bit has been readjusted just in case we are 
        //   accidentally still linking to mock particles

        //if (tempAdj < numPartTot) {
          assert(p < tempAdj);
          //if (part[p].ncnt == part[p].nadj) {
          //  printf("OVERFLOW %d\n", p);
          //} else if (part[tempAdj].ncnt == part[tempAdj].nadj) {
          //  printf("OVERFLOW %d\n", tempAdj);
          //} else {
            part[p].adj[part[p].ncnt] = tempAdj;
            part[p].ncnt++;
          if (tempAdj < numPartTot) {
            part[tempAdj].adj[part[tempAdj].ncnt] = p;
            part[tempAdj].ncnt++;
          }
          //}
        //} 
      }
//printf("ADJ %d %d %d %d %d\n", p, nin, part[p].nadj, nAdj, tempInt);
    }
  }
  fclose(fp);
*/

  clock4 = clock();
  interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
  printf(" Read voids (%.2f sec)...\n", interval);

  clock3 = clock();
// TODO - TEST NEW TREE BUILDING TECHNIQUE AND REMOVE THIS
/*
  // load voids *again* using Guilhem's code so we can get tree information
  printf(" Re-loading voids and building tree...\n");
  ZobovRep zobovCat;
  if (!loadZobov(args.voidDesc_arg, args.zone2Part_arg, args.void2Zone_arg, 
                 0, zobovCat)) {
       printf("Error loading catalog!\n");
       return -1;
  }
  VoidTree *tree;
  tree = new VoidTree(zobovCat);
  zobovCat.allZones.erase(zobovCat.allZones.begin(), zobovCat.allZones.end());
*/

  // if all of a void's zones also belong to another void,
  // it is a child of that void
  for (iVoid = 0; iVoid < numVoids; iVoid++) {
    voidID = voids[iVoid].voidID;
    int numMyZones = voids[iVoid].numZones;

    for (int iCheck = 0; iCheck < numVoids; iCheck++) {
      int numCheckZones = voids[iCheck].numZones;
      if (numMyZones > numCheckZones) continue;
      if (iVoid == iCheck) continue;

      int checkID = voids[iCheck].voidID;

      bool allZonesMatch = true;
      for (iZ = 0; iZ < void2Zones[voidID].numZones; iZ++) {
        int myZone = void2Zones[voidID].zoneIDs[iZ];
     
        bool foundMatch = false;
        for (int jZ = 0; jZ < void2Zones[checkID].numZones; jZ++) {
          int checkZone = void2Zones[checkID].zoneIDs[jZ];
          foundMatch = myZone == checkZone;
        }

        if (not foundMatch) {
          allZonesMatch = false;
          break;
        }
      }
 
      if (allZonesMatch) {
        voids[iVoid].parentID = checkID;
        voids[iCheck].numChildren++;
      }
    }

  } // end building tree

  // find level in tree
  for (iVoid = 0; iVoid < numVoids; iVoid++) {
    int level = 0;
    int parentIndx = iVoid;
    int parentID = voids[iVoid].parentID;
    while (parentID != -1) {
      level++;
      parentID = voids[parentIndx].parentID;
      // get the index of parent void
      for (int iCheck = 0; iCheck < numVoids; iCheck++){
        if (voids[iCheck].voidID == parentID) {
          parentIndx = iCheck;
          break;
        }
      }
    }
    voids[iVoid].level = level;
  }

// TODO - TEST NEW TREE BUILDING TECHNIQUE AND REMOVE THIS
/* 
  // copy tree information to our own data structures
  for (iVoid = 0; iVoid < numVoids; iVoid++) {
    voidID = voids[iVoid].voidID;
    voids[iVoid].parentID = tree->getParent(voidID);
    voids[iVoid].numChildren = tree->getChildren(voidID).size();

    // compute level in tree
    int level = 0;
    int parentID = tree->getParent(voidID);
    while (parentID != -1) {
      level++;
      parentID = tree->getParent(parentID);
    }
    voids[iVoid].level = level;
  }
*/

  clock4 = clock();
  interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
  printf(" Building void tree (%.2f sec)...\n", interval);

  printf(" Computing void properties...\n");

  // allocate space for a particle buffer
  maxNumPart = 0;
  for (iVoid = 0; iVoid < numVoids; iVoid++) {
    if (voids[iVoid].numPart > maxNumPart) maxNumPart = voids[iVoid].numPart;
  }
  voidPart = (PART *) malloc(maxNumPart * sizeof(PART));
  
  // main processing of each void
  for (iVoid = 0; iVoid < numVoids; iVoid++) {
    voidID = voids[iVoid].voidID;

    voids[iVoid].center[0] = part[voids[iVoid].coreParticle].x;
    voids[iVoid].center[1] = part[voids[iVoid].coreParticle].y;
    voids[iVoid].center[2] = part[voids[iVoid].coreParticle].z;

    voids[iVoid].voidType = CENTRAL_VOID;
 
    // first load up this void's particles into a buffer 
    clock3 = clock();
    i = 0;
    for (iZ = 0; iZ < void2Zones[voidID].numZones; iZ++) {
      zoneID = void2Zones[voidID].zoneIDs[iZ];

      for (p = 0; p < zones2Parts[zoneID].numPart; p++) {
        partID = zones2Parts[zoneID].partIDs[p];

        // something went haywire
        if (part[partID].vol < 1.e-27 && part[partID].vol > 0.) {
           printf("BAD PART!? %d %d %e", partID, numPartTot, part[partID].vol);
           exit(-1);
        }

        voidPart[i].x = part[partID].x;
        voidPart[i].y = part[partID].y;
        voidPart[i].z = part[partID].z;
        voidPart[i].vol = part[partID].vol;
        voidPart[i].edgeFlag = part[partID].edgeFlag;

        // check for edge contamination
        if (voidPart[i].edgeFlag > 0) {
          voids[iVoid].voidType = EDGE_VOID;
        }

/*
        // testing for edge contamination
        if (part[partID].vol < 1.e-27)  {
           printf("CONTAMINATED!! %d %d\n", iVoid, partID);
        } else {
           //printf("NORMAL!! %d %d %e\n", iVoid, partID, part[partID].vol);
        }
        for (int iAdj = 0; iAdj < part[partID].ncnt; iAdj++) {
          if (part[partID].adj[iAdj] > numPartTot) {
            printf("CONTAMINATED!! %d %d %d\n", iVoid, partID, iAdj);
          } 
        }
*/
        i++;
      }
    } // loading particles

    clock4 = clock();
    interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
    //printf("   %.2f for buffer\n", interval);

    // compute macrocenters
    clock3 = clock();
    double weight = 0.;
    voids[iVoid].macrocenter[0] = 0.;
    voids[iVoid].macrocenter[1] = 0.;
    voids[iVoid].macrocenter[2] = 0.;

    for (p = 0; p < voids[iVoid].numPart; p++) {
      dist[0] = voidPart[p].x - voids[iVoid].center[0];
      dist[1] = voidPart[p].y - voids[iVoid].center[1];
      dist[2] = voidPart[p].z - voids[iVoid].center[2];

      if (periodicX && fabs(dist[0]) > boxLen[0]/2.) 
          dist[0] = dist[0] - copysign(boxLen[0], dist[0]);
      if (periodicY && fabs(dist[1]) > boxLen[1]/2.) 
          dist[1] = dist[1] - copysign(boxLen[1], dist[1]);
      if (periodicZ && fabs(dist[2]) > boxLen[2]/2.) 
          dist[2] = dist[2] - copysign(boxLen[2], dist[2]);

      voids[iVoid].macrocenter[0] += voidPart[p].vol*(dist[0]);
      voids[iVoid].macrocenter[1] += voidPart[p].vol*(dist[1]);
      voids[iVoid].macrocenter[2] += voidPart[p].vol*(dist[2]);
      weight += voidPart[p].vol;
    }
    voids[iVoid].macrocenter[0] /= weight;
    voids[iVoid].macrocenter[1] /= weight;
    voids[iVoid].macrocenter[2] /= weight;
    voids[iVoid].macrocenter[0] += voids[iVoid].center[0];
    voids[iVoid].macrocenter[1] += voids[iVoid].center[1];
    voids[iVoid].macrocenter[2] += voids[iVoid].center[2];

    if (periodicX) {
      if (voids[iVoid].macrocenter[0] > ranges[0][1])
        voids[iVoid].macrocenter[0] = voids[iVoid].macrocenter[0] - boxLen[0];
      if (voids[iVoid].macrocenter[0] < ranges[0][0])
        voids[iVoid].macrocenter[0] = boxLen[0] + voids[iVoid].macrocenter[0];
    }
    if (periodicY) {
      if (voids[iVoid].macrocenter[1] > ranges[1][1])
        voids[iVoid].macrocenter[1] = voids[iVoid].macrocenter[1] - boxLen[1];
      if (voids[iVoid].macrocenter[1] < ranges[1][0])
        voids[iVoid].macrocenter[1] = boxLen[1] + voids[iVoid].macrocenter[1];
    }
    if (periodicZ) {
      if (voids[iVoid].macrocenter[2] > ranges[2][1])
        voids[iVoid].macrocenter[2] = voids[iVoid].macrocenter[2] - boxLen[2];
      if (voids[iVoid].macrocenter[2] < ranges[2][0])
        voids[iVoid].macrocenter[2] = boxLen[2] + voids[iVoid].macrocenter[2];
    }
    clock4 = clock();
    interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
    //printf("   %.2f for macrocenter\n", interval);

    // compute central density
    clock3 = clock();
    centralRad = voids[iVoid].radius/args.centralRadFrac_arg;
    centralDen = 0.;
    int numCentral = 0;
    for (p = 0; p < voids[iVoid].numPart; p++) {
      dist[0] = fabs(voidPart[p].x - voids[iVoid].macrocenter[0]);
      dist[1] = fabs(voidPart[p].y - voids[iVoid].macrocenter[1]);
      dist[2] = fabs(voidPart[p].z - voids[iVoid].macrocenter[2]);

      if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]);
      if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]);
      if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]);

      dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2);
      if (sqrt(dist2) < centralRad) numCentral += 1;
    }
    voids[iVoid].centralDen = numCentral / (volNormObs*4./3. * M_PI * 
                                            pow(centralRad, 3.));

    clock4 = clock();
    interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
    //printf("   %.2f for central density\n", interval);

    // compute maximum extent of void
    clock3 = clock();
    maxDist = 0.;
    for (p = 0; p < voids[iVoid].numPart; p++) {
      dist[0] = fabs(voidPart[p].x - voids[iVoid].macrocenter[0]);
      dist[1] = fabs(voidPart[p].y - voids[iVoid].macrocenter[1]);
      dist[2] = fabs(voidPart[p].z - voids[iVoid].macrocenter[2]);

      if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]);
      if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]);
      if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]);

      dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2);
      if (dist2 > maxDist) maxDist = dist2;
    }
    voids[iVoid].maxRadius = sqrt(maxDist);
    clock4 = clock();
    interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
    //printf("   %.2f for maximum extent\n", interval);
 
    // compute distance from macrocenter to nearest flagged particle
    clock3 = clock();
    if (args.isObservation_flag) {
      minDist = 1.e99;

      for (p = 0; p < numFlagged; p++) {
        dist[0] = voids[iVoid].macrocenter[0] - flaggedPart[p].x;
        dist[1] = voids[iVoid].macrocenter[1] - flaggedPart[p].y;
        dist[2] = voids[iVoid].macrocenter[2] - flaggedPart[p].z;

        dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2);
        if (dist2 < minDist) minDist = dist2;
      }
      voids[iVoid].nearestFlag = sqrt(minDist);
    
    } else {
      voids[iVoid].nearestFlag = 1.e99;
    }
    if (args.isObservation_flag) {
      voids[iVoid].redshiftInMpc = 
                       sqrt(pow(voids[iVoid].macrocenter[0] - boxLen[0]/2.,2) + 
                            pow(voids[iVoid].macrocenter[1] - boxLen[1]/2.,2) + 
                            pow(voids[iVoid].macrocenter[2] - boxLen[2]/2.,2));
      if (args.useComoving_flag) {
        redshift = gsl_interp_eval(interp, dL, redshifts,
                 voids[iVoid].redshiftInMpc, acc);
        voids[iVoid].redshift = redshift;
       } else {
        redshift = voids[iVoid].redshiftInMpc;
        voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.;
      }
    } else {

      voids[iVoid].redshiftInMpc = voids[iVoid].macrocenter[2];
      if (args.useComoving_flag) {
        voids[iVoid].redshift = gsl_interp_eval(interp, dL, redshifts,
                 voids[iVoid].redshiftInMpc, acc);
      } else {
        voids[iVoid].redshift = voids[iVoid].macrocenter[2]/LIGHT_SPEED*100.;
      }

      nearestEdge = 1.e99;
     
      if (!periodicX) {
        nearestEdge = fmin(nearestEdge,
                           fabs(voids[iVoid].macrocenter[0] - ranges[0][0]));
        nearestEdge = fmin(nearestEdge,
                           fabs(voids[iVoid].macrocenter[0] - ranges[0][1]));
      }
      if (!periodicY) {
        nearestEdge = fmin(nearestEdge,
                           fabs(voids[iVoid].macrocenter[1] - ranges[1][0]));
        nearestEdge = fmin(nearestEdge,
                           fabs(voids[iVoid].macrocenter[1] - ranges[1][1]));
      }
      if (!periodicZ) {
        nearestEdge = fmin(nearestEdge,
                           fabs(voids[iVoid].macrocenter[2] - ranges[2][0]));
        nearestEdge = fmin(nearestEdge,
                           fabs(voids[iVoid].macrocenter[2] - ranges[2][1]));
      }
    }

    voids[iVoid].nearestEdge = nearestEdge;

    clock4 = clock();
    interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
    //printf("   %.2f for nearest edge\n", interval);

    // compute eigenvalues and vectors for orientation and shape
    clock3 = clock();
    double inertia[9];
    for (int i = 0; i < 9; i++) inertia[i] = 0.;

    for (int p = 0; p < voids[iVoid].numPart; p++) {
      dist[0] = voidPart[p].x - voids[iVoid].macrocenter[0];
      dist[1] = voidPart[p].y - voids[iVoid].macrocenter[1];
      dist[2] = voidPart[p].z - voids[iVoid].macrocenter[2];

      if (periodicX && fabs(dist[0]) > boxLen[0]/2.)
          dist[0] = dist[0] - copysign(boxLen[0], dist[0]);
      if (periodicY && fabs(dist[1]) > boxLen[1]/2.)
          dist[1] = dist[1] - copysign(boxLen[1], dist[1]);
      if (periodicZ && fabs(dist[2]) > boxLen[2]/2.)
          dist[2] = dist[2] - copysign(boxLen[2], dist[2]);

      for (int i = 0; i < 3; i++) {
      for (int j = i; j < 3; j++) {
        if (i == j) inertia[i*3+j] += dist[0]*dist[0] + dist[1]*dist[1] + 
                                      dist[2]*dist[2];
        inertia[i*3+j] -= dist[i]*dist[j];      
      }
      }
    }
    inertia[1*3+0] = inertia[0*3+1];
    inertia[2*3+0] = inertia[0*3+2];
    inertia[2*3+1] = inertia[1*3+2];

    gsl_matrix_view m = gsl_matrix_view_array(inertia, 3, 3);
    gsl_eigen_symmv(&m.matrix, voids[iVoid].eval, voids[iVoid].evec, eigw);
   
    float a = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,1) +
                        gsl_vector_get(voids[iVoid].eval,2) - 
                        gsl_vector_get(voids[iVoid].eval,0)));
    float b = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,2) +
                        gsl_vector_get(voids[iVoid].eval,0) - 
                        gsl_vector_get(voids[iVoid].eval,1)));
    float c = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,0) +
                        gsl_vector_get(voids[iVoid].eval,1) - 
                        gsl_vector_get(voids[iVoid].eval,2)));
    float ca;
    float cb = c/b;

    float smallest = 1.e99;
    float largest = 0.;
    for (int i = 0; i < 3; i++) {
      if (gsl_vector_get(voids[iVoid].eval,i) < smallest) 
        smallest = gsl_vector_get(voids[iVoid].eval,i);
      if (gsl_vector_get(voids[iVoid].eval,i) > largest) 
        largest = gsl_vector_get(voids[iVoid].eval,i);
    }
    voids[iVoid].ellip = 1.0 - sqrt(sqrt(fabs(smallest/largest)));

    clock4 = clock();
    interval = 1.*(clock4 - clock3)/CLOCKS_PER_SEC;
  } // iVoid

  gsl_eigen_symmv_free(eigw);


  // now filter and categorize the voids based on various criteria
  int numWrong = 0;
  int numHighDen = 0;
  int numCentral = 0;
  int numEdge = 0;
  int numNearZ = 0;
  int numAreParents = 0;
  int numTooSmall = 0;

  printf(" Picking winners and losers...\n");
  printf("  Starting with %d voids\n", (int) voids.size());

  for (iVoid = 0; iVoid < voids.size(); iVoid++) {
    voids[iVoid].accepted = 1;
  }

  // toss out voids that are obviously wrong
  int iGood = 0;
  for (iVoid = 0; iVoid < voids.size(); iVoid++) {

    if (isnan(voids[iVoid].vol) || isinf(voids[iVoid].vol)) {
      numWrong++;
    } else {
      voids[iGood++] = voids[iVoid];
    }
  }
  voids.resize(iGood);
  printf("  1st filter: rejected %d obviously bad\n", numWrong); 

  // toss out voids that are too small
  iGood = 0;
  for (iVoid = 0; iVoid < voids.size(); iVoid++) {
    if (voids[iVoid].radius < args.rMin_arg) {
      numTooSmall++;
    } else {
      voids[iGood++] = voids[iVoid];
    }
  }
  voids.resize(iGood);
  printf("  2nd filter: rejected %d too small\n", numTooSmall); 

  // mark top-level voids
  numAreParents = 0;
  iGood = 0;
  for (iVoid = 0; iVoid < voids.size(); iVoid++) {
    if (voids[iVoid].parentID != -1) {
      numAreParents++;
      voids[iVoid].isLeaf = true;
     } else {
      voids[iVoid].isLeaf = false;
    }
  }

  // mark high-density voids
  for (iVoid = 0; iVoid < voids.size(); iVoid++) {
    if (voids[iVoid].centralDen > args.maxCentralDen_arg) {
      voids[iVoid].accepted = -1;
      voids[iVoid].hasHighCentralDen = true;
      numHighDen++;
    } else {
      voids[iVoid].hasHighCentralDen = false;
    }
  }


  // count voids near survey edges
  for (iVoid = 0; iVoid < voids.size(); iVoid++) {
    if (voids[iVoid].voidType == CENTRAL_VOID) {
      numCentral++;
    } else {
      numEdge++;
    }
  }

  printf("  Number kept: %d (out of %d)\n", (int) voids.size(), numVoids);
  printf("   We have %d edge voids\n", numEdge); 
  printf("   We have %d central voids\n", numCentral); 
  printf("   We have %d high central density\n", numHighDen); 
  printf("   We have %d that are not leaf nodes\n", numAreParents); 
  
  outputDir = string(args.outputDir_arg);
  sampleName = (string(args.sampleName_arg)+".out");

  dataPortions[0] = "central";
  dataPortions[1] = "all";

  printf(" Output fully untrimmed catalog...\n");
  outputVoids(outputDir, sampleName,
                numPartTot,
                voids,
                args.isObservation_flag, boxLen);

  clock2 = clock();
  printf(" Time: %f sec (for %d voids)\n", 
         (1.*clock2-clock1)/CLOCKS_PER_SEC, numVoids);
  clock1 = clock();


  printf("Done!\n");
} // end main


// ----------------------------------------------------------------------------
void openFiles(string outputDir, string sampleName, 
               int numPartTot, int numKept,
               FILE** fpOutput) {

  *fpOutput = fopen((outputDir+"/voidDatabase_"+sampleName).c_str(), "w");
  //fprintf(*fpZobov, "%d particles, %d voids.\n", numPartTot, numKept);
  fprintf(*fpOutput, "# Void ID, void type, center (x,y,z) (Mpc/h), volume (normalized), volume(Mpc/h^3), radius (Mpc/h), redshift, RA, Dec, density contrast, max extent, nearest edge, num part, parent ID, tree level, num children, central density, core particle, core density, zone vol, zone num part, num zones, void probability, ellipticity, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n");

} // end openFiles


// ----------------------------------------------------------------------------
void closeFiles(FILE* fpOutput) {

  fclose(fpOutput);
  //fclose(fpCenters);
  //fclose(fpBarycenter);
  //fclose(fpDistances);
  //fclose(fpShapes);
  //fclose(fpExtra);
  //fclose(fpSkyPositions);

} // end closeFile

// ----------------------------------------------------------------------------
void outputVoids(string outputDir, string sampleName, int numPartTot,
                vector<VOID> voids, 
                bool isObservation, double *boxLen) {

  int iVoid;
  VOID outVoid;
  FILE *fp, *fpOutput;

  openFiles(outputDir, sampleName,
              numPartTot, voids.size(),
              &fpOutput);


  for (iVoid = 0; iVoid < voids.size(); iVoid++) {
    outVoid = voids[iVoid];


     double phi = atan2(outVoid.macrocenter[1]-boxLen[1]/2.,
                        outVoid.macrocenter[0]-boxLen[0]/2.);
     if (phi < 0) phi += 2.*M_PI;
     double RA = phi * 180./M_PI;
           
     double theta = acos((outVoid.macrocenter[2]-boxLen[2]/2.) / 
                          outVoid.redshiftInMpc);
     double dec = (M_PI/2. - theta) * 180./M_PI;
 
     fprintf(fpOutput, "%d %d %e %e %e %e %e %e %e %e %e %e %e %e %d %d %d %d %e %d %e %e %d %d %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", 
             outVoid.voidID, 
             outVoid.voidType,
             outVoid.macrocenter[0],
             outVoid.macrocenter[1],
             outVoid.macrocenter[2],
             outVoid.vol, 
             4./3.*M_PI*pow(outVoid.radius, 3),
             outVoid.radius,
             outVoid.redshift, 
             RA,                  
             dec,
             outVoid.densCon, 
             outVoid.maxRadius,
             outVoid.nearestFlag,
             outVoid.numPart, 
             outVoid.parentID,
             outVoid.level,
             outVoid.numChildren,
             outVoid.centralDen,
             outVoid.coreParticle, 
             outVoid.coreDens, 
             outVoid.zoneVol, 
             outVoid.zoneNumPart, 
             outVoid.numZones, 
             outVoid.voidProb,
             outVoid.ellip,
             gsl_vector_get(outVoid.eval, 0),
             gsl_vector_get(outVoid.eval, 1),
             gsl_vector_get(outVoid.eval, 2),
             gsl_matrix_get(outVoid.evec, 0 ,0),
             gsl_matrix_get(outVoid.evec, 1 ,0),
             gsl_matrix_get(outVoid.evec, 2 ,0),
             gsl_matrix_get(outVoid.evec, 0 ,1),
             gsl_matrix_get(outVoid.evec, 1 ,1),
             gsl_matrix_get(outVoid.evec, 2 ,1),
             gsl_matrix_get(outVoid.evec, 0 ,2),
             gsl_matrix_get(outVoid.evec, 1 ,2),
             gsl_matrix_get(outVoid.evec, 2 ,2)
            );

  } // end iVoid

  closeFiles(fpOutput);

} // end outputVoids