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

using namespace std;

struct nucleon_pos
{
  float x;
  float y;
  float z;
};

struct nucleus
{
  int A;
  int Z;
  
  float R; // radius in fm

  nucleon_pos *nucleon;
};

double unirandom(float min, float max);


nucleon_pos generate_nucleon(float R)
{
  nucleon_pos pos;

  do {

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

  return pos;  
};


void initialize (struct nucleus *nucl, int z, int a)
{
  nucl->A = a;
  nucl->Z = z;
  nucl->R = pow(nucl->A, 1./3.) * 1.4;
  
  nucl->nucleon = new nucleon_pos[nucl->A];

  for (int i=0; i<nucl->A; i++) {
    nucl->nucleon[i] = generate_nucleon(nucl->R);
  }

}


// 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()
{

  nucleus Pb208;
  initialize(&Pb208, 82, 208);

  nucleus Au197;
  initialize(&Au197, 79, 197);

  nucleus C14;
  initialize(&C14, 6, 14);

  cout << "C14: Z=" << C14.Z 
       << "  N=" << (C14.A-C14.Z)
       << "  A=" << C14.A 
       << "  R=" << C14.R << "fm" 
       << endl;

  for (int i=0; i<C14.A; i++) {
    cout << setw(10) << C14.nucleon[i].x << "   "
         << setw(10) << C14.nucleon[i].z << "   "
         << setw(10) << C14.nucleon[i].y << endl;

  }
  
  


}