#include <iostream>
#include <iomanip>
#include <string>
#include <math.h>

using namespace std;


double unirandom(float min, float max);




// =====================================================================================
//  Nucleon class
// =====================================================================================

class nucleon_t
{
public:

  nucleon_t(float x=0., float y=0., float z=0.)
    : wounded(false)
  { set_xyz(x,y,z); }
  
  void set_xyz(float x, float y, float z)
  { pos[0]=x; pos[1]=y; pos[2]=z; }

  void set_wounded(bool val=true) {wounded=val;}
  bool is_wounded() { return wounded; }

  float x() { return pos[0]; }
  float y() { return pos[1]; }
  float z() { return pos[2]; }

  static nucleon_t generate_in_sphere(float Rmax);

private:
  float pos[3];
  bool wounded;

};

nucleon_t nucleon_t::generate_in_sphere(float Rmax)
{
  nucleon_t pos;

  do {

    // generate nucleons in cube with edge length 2*Rmax
    pos.set_xyz (unirandom(-Rmax, Rmax),
                 unirandom(-Rmax, Rmax),
                 unirandom(-Rmax, Rmax));
    
    // try again if nucleon is outside sphere with radius Rmax
  } while  ( sqrt(pos.x()*pos.x() + pos.y()*pos.y() + pos.z()*pos.z()) > Rmax );

  return pos;  
};


ostream& operator<<(ostream& os, nucleon_t& nucleon)
{
  os << "(" << setprecision(3) << nucleon.x() 
     << "," << setprecision(3) << nucleon.y() 
     << "," << setprecision(3) << nucleon.z() 
     << ")";

  return os;
};



// =====================================================================================
//  Nucleus class
// =====================================================================================

class nucleus_t
{
public:
  
  nucleus_t(int Z, int A);
  ~nucleus_t();

  nucleon_t& nucleon(int i) { return nucleons[i]; }

  int A() { return a; }
  int N() { return a-z; }
  int Z() { return z; }
  
  float R() {return pow(A(), 1./3.) * 1.4; }

  string symbol(); // symbol of element
      
private:

  int a;
  int z;
  
  nucleon_t *nucleons;
};

nucleus_t::nucleus_t(int ZZ, int AA)
  : z(ZZ), a(AA)
{
  
  nucleons = new nucleon_t[A()];

  for (int i=0; i<A(); i++) {
    nucleons[i] = nucleon_t::generate_in_sphere(R());
  }

}

nucleus_t::~nucleus_t()
{
  delete[] nucleons;
}

string nucleus_t::symbol()
{
  switch (Z()) {
  case  6: return "C";
  case 79: return "Au";
  case 82: return "Pb";
  default: return "XX";
  }
}



ostream& operator<<(ostream& os, nucleus_t& nuc)
{
  os << nuc.symbol() << nuc.A() << "(Z=" << nuc.Z() << ",N=" << nuc.N() << ")";

  return os;
};




// function to calculate a uniformly distributed random number 
// between min and max
double unirandom(float min, float max)
{
  return min + (max-min)*random() / (pow(2., 31.)-1.);
}


int main()
{

  // ----------------------------------------------------------------------
  // testing the nucleus and nucleon classes
  // ----------------------------------------------------------------------
  nucleus_t Pb208 ( 82, 208 );
  nucleus_t Au197 ( 79, 197 );
  nucleus_t C14   (  6,  14 );

  cout << Pb208 << " R=" << Pb208.R() << endl;
  cout << Au197 << " R=" << Au197.R() << endl;
  cout << C14   << " R=" << C14.R() << endl;

  for (int i=0; i<C14.A(); i++) {
    cout << C14.nucleon(i) << endl;
  }


  // ----------------------------------------------------------------------
  // run the actual simulation
  // ----------------------------------------------------------------------

  for (int ev=0; ev<1000; ev++) {

    // random impact parameter in 2D
    float bx = unirandom(-20., 20.);
    float by = unirandom(-20., 20.);
    
    // create two nuclei
    nucleus_t nucl1(82,208); 
    nucleus_t nucl2(82,208); 

    // initialize number of collisions
    int ncoll = 0;

    // loop over all pairs of nucleons
    for (int i=0; i< nucl1.A(); i++) {
      for (int j=0; j< nucl2.A(); j++) {

        float dx =  nucl1.nucleon(i).x() - nucl2.nucleon(j).x() + bx;
        float dy =  nucl1.nucleon(i).y() - nucl2.nucleon(j).y() + by;
    
        float dist = hypot(dx,dy);
    
        // 42 mb = 42 E-31 m^2 = 4.2 E-30 m^2 = 4.2 fm^2
        if (dist <= sqrt(4.2 / M_PI)) {
          ncoll++;
          nucl1.nucleon(i).set_wounded();
          nucl2.nucleon(j).set_wounded();
        }

      }
    }
    
    if (ncoll==0) 
      continue;

    // calculate number of participating nucleons
    int npart = 0;
    for (int i=0; i< nucl1.A(); i++) {
      if (nucl1.nucleon(i).is_wounded()) {
        npart++;
      }
    }

    for (int j=0; j< nucl2.A(); j++) {
      if (nucl2.nucleon(j).is_wounded()) {
        npart++;
      }
    }


    cout << "b = "       << setw(5) << hypot(bx,by)  
         << "  Npart = " << setw(5) << npart
         << "  Ncoll = " << setw(5) << ncoll
         << endl;
  }
  

}