time dependent frequncy in OPAL-cyc implemented

classic/5.0/src/AbsBeamline/RFCavity.cpp

@@ -565,11 +565,18 @@ void RFCavity::initialise(PartBunch *bunch, double &startField, double &endField
 }
 
 // In current version ,this function reads in the cavity voltage profile data from file.
-void RFCavity::initialise(PartBunch *bunch, const double &scaleFactor) {
+void RFCavity::initialise(PartBunch *bunch, const double &scaleFactor, std::shared_ptr<AbstractTimeDependence> freq_atd,
+                          std::shared_ptr<AbstractTimeDependence> ampl_atd, std::shared_ptr<AbstractTimeDependence> phase_atd) {
+
     using Physics::pi;
 
     RefPartBunch_m = bunch;
 
+    /// set the time dependent models
+    setAmplitudeModel(ampl_atd);
+    setPhaseModel(phase_atd);
+    setFrequencyModel(freq_atd);
+
     ifstream in(filename_m.c_str());
     if(!in.good()) {
         throw GeneralClassicException("RFCavity::initialise",
@@ -716,6 +723,18 @@ double RFCavity::getCycFrequency()const {
     return  frequency_m;
 }
 
+
+
+/**
+   \brief used in OPAL-cycl
+
+   Is called from OPAL-cycl and can handle 
+   time dependent frequency, amplitude and phase
+
+   At the moment (test) only the frequence is time
+   dependent
+
+ */
 void RFCavity::getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass, const int chargenumber) {
     using Physics::two_pi;
     using Physics::pi;
@@ -734,6 +753,9 @@ void RFCavity::getMomentaKick(const double normalRadius, double momentum[], cons
     double transit_factor = 0.0;
     double Ufactor = 1.0;
 
+    double frequency_save = frequency_m;
+    frequency_m *= frequency_td_m->getValue(t);
+
     if(gapwidth_m > 0.0) {
         transit_factor = 0.5 * frequency_m * gapwidth_m * 1.0e-3 / (c * beta);
         Ufactor = sin(transit_factor) / transit_factor;
@@ -766,10 +788,10 @@ void RFCavity::getMomentaKick(const double normalRadius, double momentum[], cons
 
     if(PID == 0) {
       *gmsg << "* Cavity " << getName() << " Phase= " << tempdegree << " [deg] transit time factor=  " << Ufactor
-	      << " dE= " << dgam *restMass * 1.0e-6 << " [MeV]"
-	      << " E_kin= " << (gamma - 1.0)*restMass * 1.0e-6 << " [MeV]" << endl;
+            << " dE= " << dgam *restMass * 1.0e-6 << " [MeV]"
+            << " E_kin= " << (gamma - 1.0)*restMass * 1.0e-6 << " [MeV] Time dep freq = " << frequency_td_m->getValue(t) << endl;
     }
-
+    frequency_m = frequency_save;
 }
 
 /* cubic spline subrutine */

classic/5.0/src/AbsBeamline/RFCavity.h

@@ -110,7 +110,13 @@ public:
 
     virtual void initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor);
 
-    virtual void initialise(PartBunch *bunch, const double &scaleFactor);
+//    virtual void initialise(PartBunch *bunch, const double &scaleFactor);
+
+    virtual void initialise(PartBunch *bunch, const double &scaleFactor, std::shared_ptr<AbstractTimeDependence> freq_atd,
+                            std::shared_ptr<AbstractTimeDependence> ampl_atd, std::shared_ptr<AbstractTimeDependence> phase_atd);
+
+
+
 
     virtual void finalise();
 

src/Algorithms/ParallelCyclotronTracker.cpp

@@ -819,28 +819,47 @@ void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {
     double phi0 = elptr->getPhi0();
     *gmsg << "* Initial RF phase (t=0)= " << phi0 << " [deg] " << endl;
 
+
+    /*
+      Setup time dependence and in case of no
+      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
+     */
+
     std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
     std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
     std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;
 
+    std::vector<double>  unityVec;
+    unityVec.push_back(1.);
+    unityVec.push_back(0.);
+    unityVec.push_back(0.);
+    unityVec.push_back(0.);
+    
     if (elptr->getFrequencyModelName() != "") {
       freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
       *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
     }
+    else
+        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
 
     if (elptr->getAmplitudeModelName() != "") {
       ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
       *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
     }
+    else
+        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
 
     if (elptr->getPhaseModelName() != "") {
       phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
       *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
     }
+    else
+        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
 
     // read cavity voltage profile data from file.
-    //    elptr->initialise(itsBunch, 1.0, freq_atd, ampl_atd, phase_atd);
-    elptr->initialise(itsBunch, 1.0); 
+    elptr->initialise(itsBunch, 1.0, freq_atd, ampl_atd, phase_atd);
+
+//    elptr->initialise(itsBunch, 1.0); 
 
     double BcParameter[8];
     for(int i = 0; i < 8; i++)