#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "isofile.h"
#include "mymalloc.h"
#include <gp.h>
#define TINY 1.e-14
#define TINY_FLT 1.e-5


void OutputIsolated(float ***level, int lx, int ly, int lz)
{
  int x,y,z;
  int number_non_zero=0, number_non_one=0,number_tiny=0, number_qone=0;
  double sum,IsolatedMass=0.,IsolatedHole=0.;
  double epsArround,epsCenter;
  float min_non_zero = 2.;
  float max_non_one = 0.;
  double TotalMass =0.;

  epsArround= GP_dfetch("epsArround",0.5,"WARNING");
  epsCenter = GP_dfetch("epsCenter",0.5,"WARNING");
 for (x = 0; x < lx - 1; x++)
    for (y = 0; y < ly - 1; y++)
      for (z = 0; z < lz - 1; z++) {
	TotalMass += (double) level[x][y][z];
	if((x > 0)&&(y > 0)&&(z > 0)&&(x < (lx-1))&&(y < (ly-1))&&(z < (lz-1))) {
	  if( ((double) level[x][y][z]) < epsCenter)  {
	    sum =  (double) level[x+1][y  ][z  ];
	    sum +=  (double) level[x-1][y  ][z  ];
	    sum +=  (double) level[x  ][y+1][z  ];
	    sum +=  (double) level[x  ][y-1][z  ];
	    sum +=  (double) level[x  ][y  ][z+1];
	    sum +=  (double) level[x  ][y  ][z-1];
	    if(sum < epsArround) {
	      IsolatedMass +=  (double) level[x][y][z]; 
	      if( ((double) level[x][y][z]) > TINY_FLT )  number_tiny++;
	      if( ( level[x][y][z]) > 0. )  {
		number_non_zero++;
		if ( level[x][y][z] < min_non_zero )   min_non_zero = level[x][y][z];
	      }
	    }
	  }
	  if( ((double) level[x][y][z]) > 1. -  epsCenter)  {
	    sum =  (double) level[x+1][y  ][z  ];
	    sum +=  (double) level[x-1][y  ][z  ];
	    sum +=  (double) level[x  ][y+1][z  ];
	    sum +=  (double) level[x  ][y-1][z  ];
	    sum +=  (double) level[x  ][y  ][z+1];
	    sum +=  (double) level[x  ][y  ][z-1];  
	    if(sum > 6. - epsArround) {
	      IsolatedHole +=  1. -(double) level[x][y][z]; 
	      if( ((double) level[x][y][z]) < 1. - TINY_FLT )  number_qone++;
	      if( ( level[x][y][z]) <  1. )  {
		number_non_one++;
		if ( level[x][y][z] > max_non_one ) max_non_one = level[x][y][z];
	      }
	    }
	  }
	}
      }
 // printf("%g \n",IsolatedMass);
#if 0
 {
   double  fraction_not_tiny = number_tiny/(float) ((lx-2)*(ly-2)*(lz-2));
   double  fraction_non_zero= number_non_zero/(float) ((lx-2)*(ly-2)*(lz-2));
   //   double  cell_size = 1./(float) ((lx-2)*(ly-2)*(lz-2));
   double  fraction_isolated_mass = IsolatedMass/TotalMass;
   double mean_mass_in_isolated_cell = IsolatedMass/number_non_zero;
   printf("%g %g %g %g %g %g %g \n",fraction_not_tiny,fraction_non_zero,
	fraction_isolated_mass,min_non_zero,TotalMass,mean_mass_in_isolated_cell,IsolatedMass);
 exit(1);
}
#else
 {
   double  fraction_not_tiny = number_qone/(float) ((lx-2)*(ly-2)*(lz-2));
   double  fraction_non_zero= number_non_one/(float) ((lx-2)*(ly-2)*(lz-2));
   //   double  cell_size = 1./(float) ((lx-2)*(ly-2)*(lz-2));
   double  fraction_isolated_mass = IsolatedHole/TotalMass;
   double mean_mass_in_isolated_cell = IsolatedHole/number_non_zero;
   printf("%g %g %g %g %g %g %g \n",fraction_not_tiny,fraction_non_zero,
	fraction_isolated_mass,min_non_zero,TotalMass,mean_mass_in_isolated_cell,IsolatedHole);
 exit(1);
}
#endif
}

// 36000 18089 -1789368629 0.00332378 0 995123