#include "G4MonopoleFieldSetup.hh"
#include "G4MonopoleFieldMessenger.hh"

#include "G4MagneticField.hh"
#include "G4UniformMagField.hh"
#include "G4FieldManager.hh"
#include "G4TransportationManager.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4MonopoleEquation.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4ChordFinder.hh"

// #include "G4ExplicitEuler.hh"
// #include "G4ImplicitEuler.hh"
// #include "G4SimpleRunge.hh"
// #include "G4SimpleHeum.hh"
#include "G4ClassicalRK4.hh"
// #include "G4HelixExplicitEuler.hh"
// #include "G4HelixImplicitEuler.hh"
// #include "G4HelixSimpleRunge.hh"
// #include "G4CashKarpRKF45.hh"
// #include "G4RKG3_Stepper.hh"

#include "G4SystemOfUnits.hh"
#include "G4TransportationManager.hh"
#include "G4PropagatorInField.hh"


//G4MonopoleFieldSetup* G4MonopoleFieldSetup::fMonopoleFieldSetup=0;

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4MonopoleFieldSetup::G4MonopoleFieldSetup()
 : fFieldManager(0),
   fChordFinder(0),
   fUsualChordFinder(0),
   fMonopoleChordFinder(0),
   fEquation(0),
   fMonopoleEquation(0),
   fMagneticField(0),
   fStepper(0),
   fMonopoleStepper(0),
   fMinStep(0.0),
   fZmagFieldValue(4*tesla),
   //fZmagFieldValue(0.000001*tesla),
   fMonopoleFieldMessenger(0)
{

  fMonopoleFieldMessenger = new G4MonopoleFieldMessenger(this);
  fFieldManager = GetGlobalFieldManager();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

// G4MonopoleFieldSetup* G4MonopoleFieldSetup::GetMonopoleFieldSetup()
// {
//    if (0 == fMonopoleFieldSetup)
//    {
//      static G4ThreadLocal G4MonopoleFieldSetup theInstance;
//      fMonopoleFieldSetup = &theInstance;
//    }
   
//    return fMonopoleFieldSetup;
// }

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4MonopoleFieldSetup::~G4MonopoleFieldSetup()
{
  delete fMonopoleFieldMessenger;
  if(fMagneticField) delete fMagneticField;
  //  if(fChordFinder)   delete fChordFinder;
  if(fUsualChordFinder)   delete fUsualChordFinder;
  if(fMonopoleChordFinder)   delete fMonopoleChordFinder;
  if(fStepper)       delete fStepper;
  if(fMonopoleStepper)  delete fMonopoleStepper;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void G4MonopoleFieldSetup::SetZMagFieldValue (G4double val)
{
  //set new magnetic field value and rebuild the field
  fZmagFieldValue = val;

  ConstructMagField();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void G4MonopoleFieldSetup::ConstructMagField()
{
  //apply a global uniform magnetic field along Z axis  
  if (fMagneticField) { delete fMagneticField; }  //delete the existing magn field

  if (fZmagFieldValue != 0.)     // create a new one if non nul
    {
      // G4cout << "Go to create new field ..." << G4endl;
      fMagneticField = new G4UniformMagField(G4ThreeVector(0., 0., 
                                                           fZmagFieldValue));        
      InitialiseAll();
    }
   else
    {
      // G4cout << "Set field = 0 ..." << G4endl;
      fMagneticField = 0;
    }

  fFieldManager->SetDetectorField(fMagneticField);
  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void G4MonopoleFieldSetup::InitialiseAll()
{
  fEquation = new G4Mag_UsualEqRhs(fMagneticField); 
  fMonopoleEquation = new G4MonopoleEquation(fMagneticField);
 
  fMinStep     = 0.00001*mm ; // minimal step of 1 mm is default
                                        
  fMonopoleStepper = new G4ClassicalRK4( fMonopoleEquation, 8 ); 
  fStepper = new G4ClassicalRK4( fEquation );                                         
                                         
  fUsualChordFinder = new G4ChordFinder( fMagneticField, fMinStep, fStepper);      
  fMonopoleChordFinder = new G4ChordFinder( fMagneticField, fMinStep, 
                                            fMonopoleStepper);
  
  G4double deltaChord        = 0.00001*mm;
  fMonopoleChordFinder->SetDeltaChord( deltaChord );

  G4double deltaOneStep      = 0.0001*mm;
  fFieldManager->SetAccuraciesWithDeltaOneStep(deltaOneStep);

  G4double deltaIntersection = 0.0001*mm;
  fFieldManager->SetDeltaIntersection(deltaIntersection);

  G4TransportationManager* transportManager =
                        G4TransportationManager::GetTransportationManager();

  G4PropagatorInField* fieldPropagator =
                                   transportManager->GetPropagatorInField();

  G4double epsMin            = 1e-6*mm;
  G4double epsMax            = 1e-6*mm;

  fieldPropagator->SetMinimumEpsilonStep(epsMin);
  fieldPropagator->SetMaximumEpsilonStep(epsMax);


  SetStepperAndChordFinder(0);
}