//------Variable declaration and initialization------//


	//------Constants------//

	const Double_t pi = TMath::Pi();		// Pi
	const Double_t h = TMath::H();			// Planck's constant
	const Double_t kb = TMath::K();			// Boltzmann's constant
	const Double_t m_nuc = 1.67377437e-27;		// Average mass of a nucleon (proton and neutron)


	//------Atomic parameters------//

	Int_t natoms = 2000;				// Number of atoms
	Int_t nplot = 50;				// Number of atom trajectories to plot
	Int_t nucleon_num = 85;				// Number of nucleons in the atom
	Double_t gam = 6 * 1e6; 			// Linewidth in Hz ( 6 MHz )
	Int_t natoms_mod = 0;				// Number of atoms that go through the cavity mode
	Int_t atom_out; 				// Indicator atom out of molasses range
	Double_t matom;					// Mass of an atom


	//------Atom position & velocity------//

	Double_t vzmodsum = 0; 				// Sum z-velocity of atoms going through the cavity mode
	TVector3 r( 1 ); 				// Position atom
	TVector3 rini( 1 ); 				// Initial position atom in tube reference frame
	TVector3 rmol; 					// Position atom from molasses reference frame
	TVector3 rbeam; 				// Position atom in beam reference frame
	TVector3 v( 1 );				// Velocity atom
	TVector3 vrec_abs( 1 );				// Absorption recoil velocity
	TVector3 vrec_spon( 1 );			// Spontaneous emission recoil velocity
	Double_t vrec; 					// Recoil velocity magnitude
	Double_t t; 					// Time


	//------Vacuum chamber geometry------//

	Double_t molchamrad = 4 * 0.0254;		// Inner radius of the molasses chamber
	Double_t winrad = 1.5 * 0.0254 / 2;		// Diameter of chamber windows
	Double_t motchamsize = 3.680 * 0.0254;		// MOT chamber size
	Double_t hole = 1e-3;				// Hole in mirror diameter
	Double_t tube;					// In vacuo tube length
	Double_t thetam;				// Polar angle corresponding to window aperture


	//------LVIS parameters------//

	Double_t theta0 = 69.075 * pi / 180;		// LVIS initial polar angle
	Double_t phi0 = 45 * pi / 180;			// LVIS initial azimuthal angle
	Double_t v0 = 14; 				// LVIS initial speed
	Double_t dv0 = 2.7;				// LVIS initial longitudinal speed fwhm
	TVector3 r0( 1 );				// Position of LVIS hole center
	Double_t d0; 					// Distance to LVIS hole center
	Double_t dvt0;					// LVIS initial transverse speed width
	Double_t dtheta0;				// LVIS initial polar angle width
	Double_t v0dev;					// Standard deviation LVIS initial speed
	Double_t theta0dev;				// Standard deviation LVIS initial polar angle


	//------Optical cavity parameters------//

	Double_t cavspa = 1e-3;				// Cavity mirror spacing
	Double_t mirdiam = 7e-3;			// Cavity mirror diameter
	Double_t rcavmod = 60 * 1e-6;			// Cavity mode radius
	Double_t zcav = molchamrad + 2.125 * 0.0254;	// Cavity mode position ( 0 , 0 , zcav )
	Double_t radxz;					// Atom radial position from cavity center reference frame


	//------Molasses beams------//

	Int_t nbeams = 4;				// Number of beams
	Double_t lam = 780 * 1e-9;			// Wavelength ( 780 nm )
//	Double_t w0 = 1e-3; 				// Molasses waist
//	Double_t waistdist = 0.04; 			// Distance to molasses waist
	Double_t width = 0;	 			// Average beam width at molasses center
	TVector3 cmol( 0.000 , 0.000 , 0.000 ); 	// Position molasses center
	Double_t dist_atom_mol;				// Distance of atom to molasses beams
	Double_t tot_s;					// Total saturation parameter
	Double_t tot_sr;				// Total scattering rate
	Double_t dt; 					// Molasses region time step
//	Double_t zr;					// Rayleigh length
	TVector3 k[nbeams];				// Wavevectors
	TVector3 pol[nbeams];				// Polarization vectors in beam reference frame
	TVector3 efield;				// Total electric field
	TVector3 w0pos[nbeams];				// Molasses beam waist position from center molasses reference frame
	TRotation rotz[nbeams];				// Rotations to transform from center molasses reference frame to beam coordinate system 
	TRotation rotx[nbeams];


	//------Random numbers------//

	TRandom3 randnum1( 0 );
	TRandom3 randnum2( 0 );
	Double_t randcont1;				// Random number container
	Double_t randcont2;				// Random number container
	Double_t prob[4];				// Storage for random numbers


	//------Misc.------//

	Double_t dtp[2]; 				// Varying precision time step
	Int_t i , j;					// Counters
	Int_t prec;					// Precision chooser
	time_t start , end;				// Time markers
	Double_t timing;				// Timing


	//------Graph definition------//

	TMultiGraph *x_vs_t = new TMultiGraph();

	TMultiGraph *y_vs_t = new TMultiGraph();

	TMultiGraph *z_vs_t = new TMultiGraph();

	TMultiGraph *vx_vs_t = new TMultiGraph();

	TMultiGraph *vy_vs_t = new TMultiGraph();

	TMultiGraph *vz_vs_t = new TMultiGraph();

	TH1D *transvel = new TH1D( "transvel" , "Transverse velocity distribution" , 100 , 0 , 1 );