/* University College London Dept of P&A Course in C/C++ 3C49 
| all rights reserved 1999 
|
| class: RanCom
|
| 
| This class provides a method to generates an angle in the range 0-180 degrees
| which has the correct probability distribution for Compton scattering.
|
| Author: P.Clarke 
*/

#ifndef  _rancom_cpp
#define  _rancom_cpp 1

#include "RanCom.h"

#define PI 3.14159

//---------------------------------------------
// Needed on some old non-compliant compilers
//typedef int bool;
//#define false 0 ;
//#define true  1 ;



// ----------------------------
// Constructor

RanCom::RanCom( float emassInput, float incidentInput )
{

  // Set the electron mass
  electronMass = emassInput ;

  // Set the incident photon energy
  incidentEnergy = incidentInput ;

}





//------------------------
// User method to generate a random angle in the range [0,180] degrees 

float RanCom::scatter( ) 
{
 
  bool notdone = true ;
  float phi  ;
  float Pmax, Pphi, Ptest ;

  // The maximum of the compton distribution is at phi = 0
  // We need this value
  Pmax = this->probability( 0. ) ;

  while( notdone ) {

    // Generate a random x value in range in [0,180]
    phi = this->uniform() * 180.0   ;

    //Calc the probability function at at phi 
    Pphi = this->probability(phi) ;

    //Generate a random number in the range [0,Pmax]
    Ptest = this->uniform() * Pmax ;

    //Test if Ptest < Pphi
    //If it is then we keep this value of phi and so can stop, 
    //otherwise continue in the loop
    if( Ptest < Pphi ) notdone = false ; 
    
  }

  return phi ;

}



//------------------------
// Private method to calculate the probability distribution as a function of angle
// for Compton scattering into angle phi. Phi is in degrees !

float RanCom::probability( float phiDeg )
{

  float alpha = incidentEnergy/electronMass ;  

  float prob, exp1, exp2, exp3, exp4 ;

  float phi = phiDeg / 360.0 * 2.0 * PI ;


  /* Parts of expression */

  exp1 = alpha * ( 1 - cos(phi) ) ;

  exp2 = pow( 1 / ( 1 + exp1 ),  2 ) ;

  exp3 = pow(exp1, 2) / ( 1 + exp1 ) ;

  exp4 = 1 + pow(cos(phi),2.) + exp3 ;

  prob = exp2 * exp4 ;


  return prob ;

}



//------------------------
// Private method to generate a uniformly distributed random number in [0,1] 

float RanCom::uniform( ) 
{
   int ran ;
   float x ;

   // Use library function to get a number in [0, RAND_MAX] 
   ran = rand() ;

   //Scale to [0,1] and make it a float
   x =  float(ran)/float(RAND_MAX) ;

   return x ;

}



#endif