#include <Randomize.hh>
#include <G4OpProcessSubType.hh>

#include <RAT/SNOMANOpRayleigh.hh>
using namespace RAT;

#include <CLHEP/Units/PhysicalConstants.h>
using namespace CLHEP;
#include <cmath>

SNOMANOpRayleigh::SNOMANOpRayleigh( const G4String& processName,
                                    G4ProcessType type)
  : G4VDiscreteProcess( processName, type )
{
  SetProcessSubType( fOpRayleigh );

  fPhysicsTable = NULL;
  BuildThePhysicsTable();
}

SNOMANOpRayleigh::~SNOMANOpRayleigh()
{
  if( fPhysicsTable != NULL )
    {
      fPhysicsTable->clearAndDestroy();
      delete fPhysicsTable;
    }
}

G4VParticleChange*
SNOMANOpRayleigh::PostStepDoIt( const G4Track& track,
                                const G4Step& step )
{
  aParticleChange.Initialize( track );
  const G4DynamicParticle* aParticle = track.GetDynamicParticle();

  const double psi = acos( 2.0 * cos( ( acos( 1.0 - 2.0 * G4UniformRand() ) - twopi ) / 3.0 ) );
  const double phi = G4UniformRand() * twopi;

  const G4ThreeVector oldMomentum = aParticle->GetMomentumDirection().unit();
  const G4ThreeVector oldPolarisation = aParticle->GetPolarization().unit();
  const G4ThreeVector e2 = oldMomentum.cross( oldPolarisation );

  // oldPolarisation is Z, oldMomentum is x and e2 y
  G4ThreeVector newMomentum = cos( psi ) * oldPolarisation + sin( psi ) * cos( phi ) * oldMomentum + sin( psi ) * sin( phi ) * e2;
  newMomentum.unit();

  G4ThreeVector newPolarisation = oldPolarisation - cos( psi ) * newMomentum;
  newPolarisation.unit();

  aParticleChange.ProposePolarization( newPolarisation );
  aParticleChange.ProposeMomentumDirection( newMomentum );

  return G4VDiscreteProcess::PostStepDoIt( track, step );
}

G4double
SNOMANOpRayleigh::GetMeanFreePath( const G4Track& track,
                                   G4double,
                                   G4ForceCondition* )
{
  const G4DynamicParticle* particle = track.GetDynamicParticle();
  const G4double photonMomentum = particle->GetTotalMomentum();
  const G4Material* material = track.GetMaterial();
  G4PhysicsOrderedFreeVector* rayleigh = static_cast<G4PhysicsOrderedFreeVector*>( (*fPhysicsTable)( material->GetIndex() ) );

  G4double rsLength = DBL_MAX;
  if( rayleigh != NULL )
    rsLength = rayleigh->Value( photonMomentum );
  return rsLength;
}

G4PhysicsOrderedFreeVector*
SNOMANOpRayleigh::CalculateRayleighMeanFreePaths( G4Material* const material ) const
{
  G4MaterialPropertiesTable* materialProperties = material->GetMaterialPropertiesTable();
  G4double betat;
  if( materialProperties->ConstPropertyExists( "ISOTHERMAL_COMPRESSIBILITY" ) )
    betat = materialProperties->GetConstProperty( "ISOTHERMAL_COMPRESSIBILITY" );
  else
    // If the material doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return
    return NULL;

  G4MaterialPropertyVector* rIndex = materialProperties->GetProperty("RINDEX");
  if( rIndex == NULL )
    // If the material doesn't have a RINDEX property vector then return
    return NULL;

  G4double meanRI = 0.0;
  if( materialProperties->ConstPropertyExists( "RINDEX_MEAN" ) )
    meanRI = materialProperties->GetConstProperty( "RINDEX_MEAN" );
  else // Must calculate the average
    {
      for( size_t iR = 0; iR < rIndex->GetVectorLength(); iR++ )
        meanRI += (*rIndex)[iR];
      meanRI /= static_cast<G4double>( rIndex->GetVectorLength() );
    }

  G4double scaleFactor = 1.0;
  if( materialProperties->ConstPropertyExists( "RS_SCALE_FACTOR" ) )
    scaleFactor= materialProperties->GetConstProperty("RS_SCALE_FACTOR" );

  G4PhysicsOrderedFreeVector* rayleighMeanFreePaths = new G4PhysicsOrderedFreeVector();
  // This calculates the meanFreePath via the Einstein-Smoluchowski formula
  const G4double c1 = scaleFactor * betat * material->GetTemperature() * k_Boltzmann / ( 6.0 * pi );
  const G4double rIndexSquared = meanRI * meanRI;
  for( size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); uRIndex++ )
    {
      const G4double energy = rIndex->Energy( uRIndex );
      const G4double xlambda = h_Planck * c_light / energy;
      const G4double c2 = std::pow( ( twopi / xlambda ), 4 );
      const G4double c3 = std::pow( ( ( rIndexSquared - 1.0 ) * ( rIndexSquared + 2.0 ) / 3.0 ), 2 );

      const G4double meanFreePath = 1.0 / ( c1 * c2 * c3 );

      rayleighMeanFreePaths->InsertValues( energy, meanFreePath );
    }

  return rayleighMeanFreePaths;
}

void
SNOMANOpRayleigh::BuildThePhysicsTable()
{
  if( fPhysicsTable )
    return;
  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
  const G4int numOfMaterials = G4Material::GetNumberOfMaterials();
  fPhysicsTable = new G4PhysicsTable( numOfMaterials );

  for( G4int iMaterial = 0; iMaterial < numOfMaterials; iMaterial++ )
    {
      G4Material* material = (*theMaterialTable)[iMaterial];
      G4MaterialPropertiesTable* materialProperties = material->GetMaterialPropertiesTable();
      G4PhysicsOrderedFreeVector* rayleigh = NULL;
      if( materialProperties != NULL )
        {
          rayleigh = materialProperties->GetProperty( "RSLENGTH" );
          if( rayleigh == NULL )
            rayleigh = CalculateRayleighMeanFreePaths( material );
        }
      fPhysicsTable->insertAt( iMaterial, rayleigh );
    }
}