CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 #include "CLHEP/GenericFunctions/RungeKuttaClassicalSolver.hh" 00002 #include "CLHEP/GenericFunctions/RKIntegrator.hh" 00003 namespace Classical { 00004 // 00005 // This is the private innards of RungeKuttaSolver 00006 // 00007 class RungeKuttaSolver::Clockwork { 00008 public: 00009 Clockwork(Genfun::GENFUNCTION gH, const PhaseSpace & mphaseSpace):H(gH),phaseSpace(mphaseSpace){} 00010 Genfun::GENFUNCTION H; 00011 const Classical::PhaseSpace & phaseSpace; 00012 Genfun::RKIntegrator integrator; 00013 std::vector<Genfun::Parameter*> startingQ; 00014 std::vector<Genfun::Parameter*> startingP; 00015 Genfun::EnergyFunction *energy; 00016 }; 00017 00018 RungeKuttaSolver::RungeKuttaSolver(Genfun::GENFUNCTION gH, const PhaseSpace & mphaseSpace):c(new Clockwork(gH,mphaseSpace)){ 00019 // 00020 // Dimension (of coords, or phase space) 00021 // 00022 const unsigned int DIM=c->phaseSpace.dim(); 00023 // 00024 // Equations of motion via hamilton's equations: 00025 // 00026 const Classical::PhaseSpace::Component & X=c->phaseSpace.coordinates(); 00027 const Classical::PhaseSpace::Component & P=c->phaseSpace.momenta(); 00028 00029 for (unsigned int i=0;i<DIM;i++) { 00030 Genfun::GENFUNCTION DXDT = c->H.partial(P[i].index()); 00031 c->startingQ.push_back(c->integrator.addDiffEquation(&DXDT,"X",c->phaseSpace.startValue(X[i]))); 00032 } 00033 for (unsigned int i=0;i<DIM;i++) { 00034 Genfun::GENFUNCTION DPDT = -c->H.partial(X[i].index()); 00035 c->startingP.push_back(c->integrator.addDiffEquation(&DPDT,"P",c->phaseSpace.startValue(P[i]))); 00036 } 00037 c->energy=NULL; 00038 00039 } 00040 RungeKuttaSolver::~RungeKuttaSolver(){ 00041 delete c->energy; 00042 delete c; 00043 } 00044 00045 Genfun::GENFUNCTION RungeKuttaSolver::equationOf(const Genfun::Variable & v) const { 00046 return *c->integrator.getFunction(v.index()); 00047 } 00048 Genfun::GENFUNCTION RungeKuttaSolver::hamiltonian() const { 00049 return c->H; 00050 } 00051 const Classical::PhaseSpace & RungeKuttaSolver::phaseSpace() const { 00052 return c->phaseSpace; 00053 } 00054 Genfun::GENFUNCTION RungeKuttaSolver::energy() const { 00055 if (!c->energy) c->energy=new Genfun::EnergyFunction(*this); 00056 return *c->energy; 00057 } 00058 00059 Genfun::Parameter *RungeKuttaSolver::createControlParameter(const std::string & variableName, 00060 double defStartingValue, 00061 double startingValueMin, 00062 double startingValueMax) const { 00063 return c->integrator.createControlParameter(variableName, defStartingValue, startingValueMin, startingValueMax) ; 00064 } 00065 00066 Genfun::Parameter *RungeKuttaSolver::takeQ0(unsigned int index) { 00067 return c->startingQ[index]; 00068 } 00069 Genfun::Parameter *RungeKuttaSolver::takeP0(unsigned int index) { 00070 return c->startingP[index]; 00071 } 00072 00073 }