// Time-stamp: </wrf/c/2007/union-sphere/us.cc, Tue,  1 May 2007, 19:18:38 EDT, http://wrfranklin.org/>
 
//Compute the area of the union of spheres.
// Read (x,y,z,r)  spheres from cin.
 
// W. Randolph Franklin,  http://wrfranklin.org/
 
//================================================================
// includes
//================================================================
 
#include <boost/multi_array.hpp> 
#include <boost/progress.hpp>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <math.h>
#include <vector>
#include <gsl/gsl_integration.h>
 
//================================================================
// usings
//================================================================
 
using namespace std;
using boost::array;
using boost::timer;
using boost::progress_timer;
using boost::progress_display;
 
 
//================================================================
// macros
//================================================================
 
#define YES
#define WRITEZ   // Write the z values used for the slices?
 
// Print an expression's name then its value, possible followed by a comma or endl.  
// Ex: cout << PRINTC(x) << PRINTN(y);
 
#define PRINT(arg)  #arg "= " << (arg)
#define PRINTC(arg)  (printcolors?redtty:"") << #arg << (printcolors?deftty:"") << "= " << (arg) << ", "
#define PRINTN(arg)  (printcolors?redtty:"") << #arg << (printcolors?deftty:"") << "= " << (arg)  << endl
 
 
// Execute a block of code then print the code and its time to cout.
// Ex:   RT(do_process()) ;
 
#define RT(arg) { boost::progress_timer bpt; {arg;}; cout << "Time to exec " << (printcolors?redtty:"") << #arg << (printcolors?deftty:"") << "is ";  }
// #define RT(arg) { arg;   }
 
 
// Compact for loop
// Ex:   FOR(i,0,3,a[i]=0);
 
#define FOR(foriterator,formin,formaxplus1,forbody) {for(int foriterator=formin; foriterator<formaxplus1; foriterator++) {forbody;};}
 
 
 
//================================================================
// classes
//================================================================
 
class TCircle {
public:
  double x, y, r;
  int sph;          // corresponding sphere
};
 
 
class TSphere {
public:
  double x, y, z, r;
  int cir;           // corresponding circle
};
 
 
 
//================================================================
// typedefs
//================================================================
 
typedef array<double,2> Tpt;
typedef vector<TCircle> TCircles;
typedef vector<TSphere> TSpheres;
 
 
 
//================================================================
// variables
//================================================================
 
const string redtty("\033[1;31m");   // tell tty to switch to bold red
const string bluetty("\033[1;34m");   // tell tty to switch to bold blue
const string greentty("\033[1;32m");   // tell tty to switch to bold green
const string yellowtty("\033[1;33m");   // tell tty to switch to bold yellow
const string deftty("\033[0m");      // tell tty to switch back to default color
const string blackbgtty("\033[40m");      // tell tty to switch background to black
const string defbgtty("\033[47m");      // tell tty to switch background to white

int printcolors(1);   // Use the above colors in PRINTC etc. 

const int debug(0);
const int debugtime(0);
 
const double Pi(3.14159265358979323846264338327950288419716939937510);
 
int NSpheres;
TSpheres Spheres;
 
int NComputeOneSliceArea(0);
int NComputePerimeter(0);
int NComputePerimeterCircles(0);
int NComputePerimeterClose(0);
int NSliceSphere(0);
int NIsPointHidden(0);
int NIsPointHiddenCircle(0);
 
double MaxInputRadius(0);
 
array<array<double, 3>, 2> MolBounds;
array<double,3> MolRange;
 
 
const int MSpheres(5000);
const int MClosePerSphere(100);
array<array<int,MClosePerSphere>,MSpheres> CloseSpheres;
array<int,MSpheres> NClosePerSphere;
int NCloseSpherePairs(0);

const double GivenRelError(0.3);

// No time to debug this:
//vector<vector<int> * > CloseSpheres;   // The spheres close enough to each sphere potentially to intersect it.

#ifdef WRITEZ
ofstream  zstream("z.dat");
#endif

//================================================================
// class output routines
//================================================================
 
ostream &operator<<(ostream &o, const Tpt a) {
  o << "(" << a[0] << ", " << a[1] << ')';
  return o;
}
 
 
ostream &operator<<(ostream &o, const array<Tpt,2> a) {
  o << "(" << a[0] << ", " << a[1] << ')';
  return o;
}
 
 
ostream &operator<<(ostream &o, const TCircle a) {
  o << "(" << a.x << ", " << a.y << ", " <<  a.r << ", " << a.sph << ')';
  return o;
}
 
 
ostream &operator<<(ostream &o, const TSphere a) {
  o << "(" << a.x << ", " << a.y << ", " << a.z << ", " << a.r << ", " << a.cir << ')';
  return o;
}
 
 
template<class T> 
ostream &operator<<(ostream &o, const vector<T> v) {
  o << "(";
  for (int i=0; i<v.size(); i++) {
    o << v[i];
    if (i<v.size()-1) o << ',';
  }
  o << ')';
  return o;
}
 
 
//================================================================
//  class arithmetic  routines
//================================================================
 
 
// Multiplying a scalar times a circle multiplies the circle's center point, but doesn't use its radius.
 
inline const Tpt operator*(const double a, const TCircle &b) {
  Tpt c;
  c[0] = a*b.x;
  c[1] = a*b.y;
  return c;
}
 
 
// Subtracting 2 circles subtracts their centers.
 
inline const Tpt operator-(const TCircle &a, const TCircle &b) {
  Tpt c;
  c[0] = a.x - b.x;
  c[1] = a.y - b.y;
  return c;
}
 
 
inline const Tpt operator+(const Tpt &a, const Tpt b) {
  Tpt c(a);
  c[0] += b[0];
  c[1] += b[1];
  return c;
}
 
 
inline const Tpt operator-(const Tpt &a, const Tpt b) {
  Tpt c(a);
  c[0] -= b[0];
  c[1] -= b[1];
  return c;
}
 
 
inline const Tpt operator-(const Tpt &a, const TCircle b) {
  Tpt c(a);
  c[0] -= b.x;
  c[1] -= b.y;
  return c;
}
 
 
template<typename T>
inline const Tpt operator*(const Tpt &a, const T r) {
  Tpt c(a);
  c[0] *= r;
  c[1] *= r;
  return c;
}
 
 
template<typename T>
inline const Tpt operator*(const T r, const Tpt &a) {
  Tpt c(a);
  c[0] *= r;
  c[1] *= r;
  return c;
}
 
 
template<typename T>
inline const array<T,2> operator-(const array<T,2> a) {
  array<T,2> b;
  b[0] = -a[0];
  b[1] = -a[1];
  return b;
}
 
 
template<typename T>   // Add a vector to the negative of each column of an array.
inline const array<T,2> operator-(const T &a, const array<T,2> b) {
  array<T,2> c(-b);
  c[0] += a;
  c[1] += a;
  return c;
}
 
 
 
//================================================================
// general procedures
//================================================================
 
 
inline const Tpt rot90(const Tpt &p) {
  Tpt q;
  q[0] = p[1];
  q[1] = -p[0];
  return q;
}
 
 
#ifdef SWAPNEEDED
template <class T> 
void swap(T a, T b) {
  T c(a);
  a = b;
  b = c;
}
#endif
 
 
template<class T>
inline const T sqr(const T x) {
  return x*x;
}
 
 
inline const double distsq(const Tpt& p, const Tpt& q) {
  return sqr(p[0]-q[0]) + sqr(p[1]-q[1]);
}
 
 
inline const double distsq(const TCircle& p, const TCircle& q) {
  return sqr(p.x-q.x) + sqr(p.y-q.y);
}
 
 
inline const double distsq(const TSphere& p, const TSphere& q) {
  return sqr(p.x-q.x) + sqr(p.y-q.y) + sqr(p.z-q.z);
}
 
 
inline const double distsq(const TCircle& p, const Tpt& q) {
  return sqr(p.x-q[0]) + sqr(p.y-q[1]);
}
 
 
inline const double atanp(const Tpt& v) {
  return atan2(v[1], v[0]);
}
 
 
template <class T> 
inline const int signum(const T x) {
  return x>0?1:x==0?0:-1;
}
 
 
template <class T>    // Reduce an angle to the interval (-Pi,Pi].
inline const T norm1(const T x) {
  T y(x);
  if (y > Pi) y -= 2*Pi;
  else if (y < -Pi) y += 2*Pi;
  return y;
}
 
 
template <class T>    // Reduce an angle to the interval [0,2{i)
inline const T norm2(const T x) {
  T y(x);
  if (y >= 2*Pi) y -= 2*Pi;
  else if (y < 0) y += 2*Pi;
  return y;
}
 
 
 
 
//================================================================
// specific procedures
//================================================================
 
 
void ReadMolecule() {
  Spheres.clear();
  TSphere sphere;
  for (int i=0;;i++) {
    cin >> sphere.x >> sphere.y >> sphere.z >> sphere.r;
    if (!cin.good()) break;
    Spheres.push_back(sphere);
    //  cout << PRINTC(i) << PRINTN(Spheres[i]);
    if (sphere.r > MaxInputRadius) MaxInputRadius = sphere.r;
  }
}
 
 

// Swap coords to make the MolZ dimension have the greatest extent
void OrientMolecule() {

  FOR(i, 0, 3, MolBounds[0][i] = std::numeric_limits<double>::max());
  FOR(i, 0, 3, MolBounds[1][i] = std::numeric_limits<double>::min());

  for (int i=0; i<NSpheres; i++) {
    if (Spheres[i].x < MolBounds[0][0]) MolBounds[0][0] = Spheres[i].x;
    if (Spheres[i].x > MolBounds[1][0]) MolBounds[1][0] = Spheres[i].x;
    if (Spheres[i].y < MolBounds[0][1]) MolBounds[0][1] = Spheres[i].y;
    if (Spheres[i].y > MolBounds[1][1]) MolBounds[1][1] = Spheres[i].y;
    if (Spheres[i].z < MolBounds[0][2]) MolBounds[0][2] = Spheres[i].z;
    if (Spheres[i].z > MolBounds[1][2]) MolBounds[1][2] = Spheres[i].z;
  }

  FOR(i, 0, 3, MolRange[i] = MolBounds[1][i] - MolBounds[1][i]);

  if ( MolRange[2] >= MolRange[0] && MolRange[2] >= MolRange[1]) {
    return;
  }
  else if ( MolRange[1] > MolRange[2] && MolRange[1] > MolRange[0]) {  // swap Y and Z
    for (int i=0; i<NSpheres; i++) {
      swap(Spheres[i].y, Spheres[i].z);
    }
    return;
  } else {  // Swap X and Z
    for (int i=0; i<NSpheres; i++) {
      swap(Spheres[i].x, Spheres[i].z);
    }
    return;
  }
}
 


void FindCloseSpheres() {

  int debugfcs(0);

  if (debugfcs) cout << bluetty << "FindCloseSpheres" << deftty << endl;

  //  CloseSpheres.reserve(NSpheres);
  //  CloseSpheres.clear();
  //  if (debugfcs) cout << PRINTN(CloseSpheres.size());

  for (int i=0; i < NSpheres; i++) {
    //    CloseSpheres.push_back(new vector<int>);
    //  if (debugfcs) cout << PRINTN(CloseSpheres[i]->size());
    //    CloseSpheres[i]-> clear();            // Probably not needed.
    NClosePerSphere[i] = 0;
  }
  //  if (debugfcs) cout << PRINTN(CloseSpheres.size());

  for (int i=0; i < NSpheres; i++) {
    for (int j = i+1; j < NSpheres; j++) {
      const double dd(distsq(Spheres[i],Spheres[j]));
      const double rr(sqr(Spheres[i].r + Spheres[j].r));
      // Should also exclude spheres too close for their boundaries to intersect, but this probably doesn't happen with molecules.
      if (dd < rr*1.0000001) {
	//	CloseSpheres[i] -> push_back(j);
	//	CloseSpheres[j] -> push_back(i);
	if (NClosePerSphere[i]>=MClosePerSphere || NClosePerSphere[j]>=MClosePerSphere) {
	  cout << "ERROR overflow in CloseSpheres" << PRINTC(MClosePerSphere) << PRINTC(i) 
	       << PRINTC(j) << PRINTC(NClosePerSphere[i]) << PRINTC(NClosePerSphere[j]) << PRINTC(Spheres[i]) << PRINTN(Spheres[j]);
	  NClosePerSphere[i] = min(NClosePerSphere[i], MClosePerSphere-1);
	  NClosePerSphere[j] = min(NClosePerSphere[j], MClosePerSphere-1);
	}
	CloseSpheres[i][NClosePerSphere[i]++] = j;
	CloseSpheres[j][NClosePerSphere[j]++] = i;
	NCloseSpherePairs++;
	//	if (debugfcs) cout << "Close pair: " << i << ' ' << j << ' ' << PRINTC(Spheres[i]) << PRINTC(Spheres[j]) << PRINTC(CloseSpheres[i]->size()) << PRINTN(CloseSpheres[j]->size());
      }
    }
  }
  cout << PRINTN(NCloseSpherePairs);
}
 
 
 
inline const int IsPointHidden(const double z, const TCircles &C, const int &nc,  const int &ic, const int &jc, const Tpt &p) {
 
  // Is point p hidden by any circle other than ic or jc?   Ok for ic==jc.   p is on circles ic and jc.
  // p has only (x,y); z is a separate var.

  NIsPointHidden++;
 
  const int is(C[ic].sph);
  const int nclose(NClosePerSphere[is]);
  for (int jjc=0; jjc < nclose; jjc++) {
    NIsPointHiddenCircle++;
    const TSphere &s(Spheres[CloseSpheres[is][jjc]]);
    if (s.cir==ic || s.cir == jc) continue;
    if (sqr(p[0]-s.x) + sqr(p[1]-s.y) + sqr(z-s.z) <= sqr(s.r)) return 1;
  }
  return 0;
}
 
 
 
 
 
 
 
double ComputePerimeter(const double z, const TCircles &C, const int nc) {
  NComputePerimeter++;
  int NXsects(0);
  vector<Tpt> Xsects;
  int NVisXsects(0);
  int NVisRtPt(0);
  double Perimeter(0);    // of the union.  Output value.
  int NPairs(0);

  //   Find and process intersections
 
  for (int ic=0; ic < nc ; ic++) {
    NComputePerimeterCircles++;
    const TCircle &c(C[ic]);
    if (c.r < 0) exit(2);
    const int is(c.sph);
    const int nclose(NClosePerSphere[is]);
    for (int jjc=0; jjc < nclose; jjc++) {
      NComputePerimeterClose++;
      const int jc = Spheres[CloseSpheres[is][jjc]].cir;
      if (jc<ic) continue;
      NPairs++;
      array<Tpt,2> Intersections;
 
      const double x2(distsq(c, C[jc]));
      const double R(c.r);
      const double r(C[jc].r);
 
      if (x2 <= sqr(R-r)) {
	continue;
      }
      if (x2 >= sqr(R+r)) {
	continue;
      }
 
      const double a((sqr(R) - sqr(r))/(2*x2) + 0.5);
      const Tpt mid((1-a)*c + a*(C[jc]));
      const double ybx(sqrt(sqr(R)/x2-sqr(a)));
      const Tpt y(ybx*(C[jc]-c));
      const Tpt yrot90(rot90(y));
 
      Intersections[0] = mid + yrot90;
      Intersections[1] = mid - yrot90;
      NXsects += 2;
      for (int ii=0; ii<2; ii++) {
	if (IsPointHidden(z, C, nc, ic, jc, Intersections[ii])) continue;
	NVisXsects++;
 
	const double a1(atanp(Intersections[ii]-c));
	const double a2(atanp(Intersections[ii]-C[jc]));
	const double da1(norm1(a2-a1));
	const double da2(r*norm2(a2) - R*norm2(a1));
	const double p(signum(da1)*da2);
	Perimeter += p;
      }
    }
  }
 
  //  Process far right x pts
 
  int NExt(0);
  for (int ic=0; ic < nc ; ic++) {
    const TCircle &c(C[ic]);
    NExt++;
    const Tpt p = {c.x, c.y};
    const double R(c.r);
    Tpt RtPt(p);
    RtPt[0] += R;     // This is the point on the right extreme of the circle.
    if (!IsPointHidden(z, C, nc, ic, ic, RtPt)) {
      NVisRtPt++;
      Perimeter += 2*Pi* Spheres[c.sph].r;
    }
  }
  return Perimeter;
}
 
 
array<double,2> ComputeZBounds() {
  array<double,2> b = {std::numeric_limits<double>::max(), 
		       std::numeric_limits<double>::min()};
  for (int i=0; i<NSpheres; i++) {
    if (Spheres[i].z < b[0]) b[0] = Spheres[i].z;
    if (Spheres[i].z > b[1]) b[1] = Spheres[i].z;
  }
  return b;
}
 
 
 
// Intersect a constant Z plane with the sphere to return a circle.   If the plane does not intersect the circle, then return the circle's radius as NaN.
 
TCircle SliceSphere(const TSphere &s, const double z) {
  NSliceSphere++;
  TCircle c;
  const double dd(sqr(s.r) - sqr(s.z - z));
  c.x = s.x;
  c.y = s.y;
  // I would have liked to use the fact that sqrt(-1) returns NaN, but using NaN is too slow.
  c.r = (dd<0) ? -999.0 : sqrt(dd);    
#ifdef NONO
  cout << bluetty << "SliceSphere: " << PRINTN(c);
#endif
  return c;
}


inline double ComputeOneSliceArea(const double z, void *params) {
#ifdef DEBUGCOS
  int debugcos(0);
  if (debugcos) cout << bluetty << "Starting ComputeOneSliceArea:" << deftty << PRINTN(z);
#endif
  TCircles C(NSpheres);
  int nc(0);
  NComputeOneSliceArea++;
  if (NComputeOneSliceArea%500 == 0) cout << PRINTN(NComputeOneSliceArea);
  // The slicing plane may not intersect some spheres.  In that case, the returned circle has a
  // radius of NaN.   Store such circles.

#ifdef WRITEZ
  zstream << z << endl;
#endif

  for (int ic=0; ic<NSpheres; ic++) {
    C[nc] = SliceSphere(Spheres[ic], z);
    //    if (debugcos) cout << PRINTC(nc) << PRINTN(C[nc]);
    if (C[nc].r < 0) Spheres[ic].cir = -999;
    else {
      Spheres[ic].cir = nc;
      C[nc].sph = ic;
      nc++;
    }
  }
#ifdef DEBUGCOS
  if (debugcos) {
    cout << PRINTN(nc);
    FOR(i, 0, NSpheres, cout << PRINTC(i) << PRINTC(Spheres[i]) << PRINTN(C[i]));
  }
#endif
  double p;
  if (nc==0) p=0;
  else p = ComputePerimeter(z, C, nc);
#ifdef DEBUGCOS
  if (debugcos) cout << bluetty << "Ending ComputeOneSliceArea:" << deftty << PRINTN(p);
#endif
  return p;
}


double ComputeArea() {
  //  cout << bluetty << "ComputeArea:" << deftty << endl;
  array<double,2> zbounds(ComputeZBounds());
  int n(NSpheres*10+1);   // Initial number of slices.   Must be odd for Simpson's rule.
  const double zl(zbounds[0] - MaxInputRadius);
  const double zh(zbounds[1] + MaxInputRadius);
  cout << PRINTC(zl) << PRINTN(zh);

  const int Mws(10000);
  gsl_integration_workspace * w 
    = gsl_integration_workspace_alloc (Mws);
       
  double Area, error;
  double alpha = 1.0;
     
  gsl_function F;
  F.function = &ComputeOneSliceArea ;
  F.params = &alpha;

  const double abs_err(12*MaxInputRadius*MaxInputRadius*GivenRelError);
     
  //  gsl_integration_qags (&F, zl, zh, abs_err, GivenRelError, Mws, w, &area, &error); 

  gsl_integration_qag (&F, zl, zh, abs_err, GivenRelError, Mws, GSL_INTEG_GAUSS15, w, &Area, &error); 

  cout << PRINTC(abs_err) << PRINTC(GivenRelError) << PRINTC(Area) << PRINTC(error) << PRINTN(w->size);
  gsl_integration_workspace_free (w);
  return Area;
}


int main(int argc, char** argv) {
  boost::progress_timer t;  // start timing; value will print at end.
  boost::timer t2;
  cout << blackbgtty << yellowtty << endl << "       us - area of union of spheres\n " << defbgtty << deftty << endl;

  ReadMolecule();
  NSpheres = Spheres.size();
  if (NSpheres > MSpheres) {
    cout << PRINTC(NSpheres) << " > " << PRINTC(MSpheres) << ". Please recompile\n";
    exit(1);
  }
  cout << PRINTC(MaxInputRadius) << PRINTN(NSpheres);
  OrientMolecule();
  FindCloseSpheres();
  double Area=ComputeArea();
  cout << PRINTC(NComputeOneSliceArea) << PRINTC(NComputePerimeter) 
       << PRINTC(NComputePerimeterCircles) << PRINTC(NComputePerimeterClose) 
       <<PRINTC(NIsPointHidden) << PRINTC(NIsPointHiddenCircle)
       << PRINTC(NSliceSphere) << PRINTN(Area);

  printcolors = 0;
  cerr << PRINTC(NSpheres) << PRINTC(MaxInputRadius) << PRINTN(NCloseSpherePairs);
  cerr << PRINTC(GivenRelError) << PRINTC(NComputeOneSliceArea) << PRINTN(NComputePerimeter);
  cerr << PRINTC(NComputePerimeterCircles) << PRINTC(NComputePerimeterClose) << PRINTN(NSliceSphere);
  const double Time(t2.elapsed());
  cerr << PRINTC(NIsPointHidden) << PRINTC(NIsPointHiddenCircle) 
       << PRINTC(Time)
       << "    " << PRINTN(Area);

  cout << "\nTotal time is ";  // following value will print  on exit
}