RFCavity.h 14.8 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
#ifndef CLASSIC_RFCavity_HH
#define CLASSIC_RFCavity_HH

// ------------------------------------------------------------------------
// $RCSfile: RFCavity.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: RFCavity
//   Defines the abstract interface for an accelerating structure.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------

24

gsell's avatar
gsell committed
25
#include "AbsBeamline/Component.h"
26
#include "Algorithms/AbstractTimeDependence.h"
gsell's avatar
gsell committed
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
#include "Physics/Physics.h"

class Fieldmap;

// Class RFCavity
// ------------------------------------------------------------------------
/// Interface for RF cavity.
//  Class RFCavity defines the abstract interface for RF cavities.


class RFCavity: public Component {

public:

    enum CavityType { SW, SGSW };
    /// Constructor with given name.
    explicit RFCavity(const std::string &name);

    RFCavity();
    RFCavity(const RFCavity &);
    virtual ~RFCavity();

    /// Apply visitor to RFCavity.
50
    virtual void accept(BeamlineVisitor &) const override;
gsell's avatar
gsell committed
51 52 53 54 55 56 57 58 59 60 61 62 63

    /// Get RF amplitude.
    virtual double getAmplitude() const = 0;

    /// Get RF frequencey.
    virtual double getFrequency() const = 0;

    /// Get RF phase.
    virtual double getPhase() const = 0;

    void dropFieldmaps();

    /// Set the name of the field map
64
    virtual void setFieldMapFN(std::string fmapfn);
gsell's avatar
gsell committed
65

66
    virtual std::string getFieldMapFN() const;
gsell's avatar
gsell committed
67

68 69 70 71
    virtual void setAmplitudem(double vPeak);
    virtual double getAmplitudem() const;
    virtual void setAmplitudeError(double vPeakError);
    virtual double getAmplitudeError() const;
gsell's avatar
gsell committed
72

73
    virtual void setFrequencym(double freq);
gsell's avatar
gsell committed
74 75 76

    void setFrequency(double freq);

77
    virtual double getFrequencym() const ;
gsell's avatar
gsell committed
78

79
    virtual void setPhasem(double phase);
gsell's avatar
gsell committed
80

81
    virtual double getPhasem() const;
82 83
    double getPhasem(double t) const;

84 85
    virtual void setPhaseError(double phaseError);
    virtual double getPhaseError() const;
gsell's avatar
gsell committed
86 87 88 89 90

    void setCavityType(std::string type);

    std::string getCavityType() const;

91
    virtual void setFast(bool fast);
gsell's avatar
gsell committed
92

93
    virtual bool getFast() const;
gsell's avatar
gsell committed
94

95
    virtual void setAutophaseVeto(bool veto = true);
gsell's avatar
gsell committed
96

97
    virtual bool getAutophaseVeto() const;
gsell's avatar
gsell committed
98

99
    virtual double getAutoPhaseEstimate(const double & E0, const double & t0, const double & q, const double & m);
100
    virtual double getAutoPhaseEstimateFallback(double E0, double t0, double q, double m);
gsell's avatar
gsell committed
101

102
    virtual std::pair<double, double> trackOnAxisParticle(const double & p0,
103 104 105
                                                          const double & t0,
                                                          const double & dt,
                                                          const double & q,
106
                                                          const double & mass,
kraus's avatar
kraus committed
107
                                                          std::ofstream *out = NULL);
gsell's avatar
gsell committed
108

109
    virtual void addKR(int i, double t, Vector_t &K) override;
gsell's avatar
gsell committed
110

111
    virtual void addKT(int i, double t, Vector_t &K) override;
gsell's avatar
gsell committed
112

113 114 115
    virtual bool apply(const size_t &i,
                       const double &t,
                       Vector_t &E,
116
                       Vector_t &B) override;
gsell's avatar
gsell committed
117

118 119 120 121
    virtual bool apply(const Vector_t &R,
                       const Vector_t &P,
                       const double &t,
                       Vector_t &E,
122
                       Vector_t &B) override;
123

124 125 126 127
    virtual bool applyToReferenceParticle(const Vector_t &R,
                                          const Vector_t &P,
                                          const double &t,
                                          Vector_t &E,
128
                                          Vector_t &B) override;
129

130
    virtual void initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) override;
131

132
    virtual void initialise(PartBunchBase<double, 3> *bunch,
133 134 135
                            std::shared_ptr<AbstractTimeDependence> freq_atd,
                            std::shared_ptr<AbstractTimeDependence> ampl_atd,
                            std::shared_ptr<AbstractTimeDependence> phase_atd);
gsell's avatar
gsell committed
136

137
    virtual void finalise() override;
gsell's avatar
gsell committed
138

139
    virtual bool bends() const override;
gsell's avatar
gsell committed
140

141
    virtual void goOnline(const double &kineticEnergy) override;
gsell's avatar
gsell committed
142

143
    virtual void goOffline() override;
gsell's avatar
gsell committed
144

145 146
    virtual void setDesignEnergy(const double& ekin, bool changeable = true)  override;
    virtual double getDesignEnergy() const override;
147

gsell's avatar
gsell committed
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
    void setRmin(double rmin);

    void setRmax(double rmax);

    void setAzimuth(double angle);

    void setPerpenDistance(double pdis);

    void setGapWidth(double gapwidth);

    void setPhi0(double phi0);

    virtual double getRmin() const;

    virtual double getRmax() const;

    virtual double getAzimuth() const;

    virtual double getCosAzimuth() const;

    virtual double getSinAzimuth() const;

    virtual double getPerpenDistance() const;

    virtual double getGapWidth() const;

    virtual double getPhi0() const;

176
    virtual void setComponentType(std::string name) override;
gsell's avatar
gsell committed
177

178
    virtual std::string getComponentType()const override;
gsell's avatar
gsell committed
179 180 181

    virtual double getCycFrequency()const;

182
    void getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass,const int chargenumber);
gsell's avatar
gsell committed
183 184 185

    double spline(double z, double *za);

186
    virtual ElementBase::ElementType getType() const override;
gsell's avatar
gsell committed
187

188
    virtual void getDimensions(double &zBegin, double &zEnd) const override;
gsell's avatar
gsell committed
189

190
    virtual bool isInside(const Vector_t &r) const override;
191

192 193 194
    void setAmplitudeModel(std::shared_ptr<AbstractTimeDependence> time_dep);
    void setAmplitudeModelName(std::string name);
    std::string getAmplitudeModelName();
195

196 197 198
    void setPhaseModel(std::shared_ptr<AbstractTimeDependence> time_dep);
    void setPhaseModelName(std::string name);
    std::string getPhaseModelName();
199

200 201 202
    void setFrequencyModel(std::shared_ptr<AbstractTimeDependence> time_dep);
    void setFrequencyModelName(std::string name);
    std::string getFrequencyModelName();
203

204
    virtual double getElementLength() const override;
205
    virtual void getElementDimensions(double &begin,
206
                                      double &end) const override;
207

208 209
    virtual CoordinateSystemTrafo getEdgeToBegin() const override;
    virtual CoordinateSystemTrafo getEdgeToEnd() const override;
210

211
protected:
212 213 214 215 216 217
    std::shared_ptr<AbstractTimeDependence> phase_td_m;
    std::string phase_name_m;
    std::shared_ptr<AbstractTimeDependence> amplitude_td_m;
    std::string amplitude_name_m;
    std::shared_ptr<AbstractTimeDependence> frequency_td_m;
    std::string frequency_name_m;
218

gsell's avatar
gsell committed
219
    std::string filename_m;             /**< The name of the inputfile*/
220

gsell's avatar
gsell committed
221
    double scale_m;              /**< scale multiplier*/
222 223 224 225 226
    double scaleError_m;         /**< additive scale error*/
    double phase_m;              /**< phase shift of time varying field (rad)*/
    double phaseError_m;         /**< phase shift error (rad)*/
    double frequency_m;          /**< Read in frequency of time varying field(Hz)*/

227 228 229 230
    bool fast_m;
    bool autophaseVeto_m;

    double designEnergy_m;
231

232
    Fieldmap* fieldmap_m;
233
    double length_m;
gsell's avatar
gsell committed
234
    double startField_m;         /**< starting point of field(m)*/
235 236

private:
gsell's avatar
gsell committed
237 238 239 240 241 242 243 244 245 246 247
    CavityType type_m;

    double rmin_m;
    double rmax_m;
    double angle_m;
    double sinAngle_m;
    double cosAngle_m;
    double pdis_m;
    double gapwidth_m;
    double phi0_m;

248 249 250
    std::unique_ptr<double[]> RNormal_m;
    std::unique_ptr<double[]> VrNormal_m;
    std::unique_ptr<double[]> DvDr_m;
gsell's avatar
gsell committed
251 252
    int num_points_m;

253 254
    double getdE(const int & i,
                 const std::vector<double> & t,
gsell's avatar
gsell committed
255 256 257 258 259
                 const double & dz,
                 const double & phi,
                 const double & frequency,
                 const std::vector<double> & F) const;

260
    double getdT(const int & i,
gsell's avatar
gsell committed
261 262 263
                 const std::vector<double> & E,
                 const double & dz,
                 const double mass) const;
264

gsell's avatar
gsell committed
265
    double getdA(const int & i,
266
                 const std::vector<double> & t,
gsell's avatar
gsell committed
267 268 269 270 271
                 const double & dz,
                 const double & frequency,
                 const std::vector<double> & F) const;

    double getdB(const int & i,
272
                 const std::vector<double> & t,
gsell's avatar
gsell committed
273 274 275 276 277 278 279 280 281 282
                 const double & dz,
                 const double & frequency,
                 const std::vector<double> & F) const;

    // Not implemented.
    void operator=(const RFCavity &);
};

inline
double RFCavity::getdE(const int & i,
283
                       const std::vector<double> & t,
gsell's avatar
gsell committed
284 285 286 287
                       const double & dz,
                       const double & phi,
                       const double & frequency,
                       const std::vector<double> & F) const {
288
    return dz / (frequency * frequency * (t[i] - t[i-1]) * (t[i] - t[i-1])) *
gsell's avatar
gsell committed
289 290 291 292 293
        (frequency * (t[i] - t[i-1]) * (F[i] * sin(frequency * t[i] + phi) - F[i-1] * sin(frequency * t[i-1] + phi)) +
         (F[i] - F[i-1]) * (cos(frequency * t[i] + phi) - cos(frequency * t[i-1] + phi)));
}

inline
294
double RFCavity::getdT(const int & i,
gsell's avatar
gsell committed
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322
                       const std::vector<double> & E,
                       const double & dz,
                       const double mass) const {
    double gamma1  = 1. + (19. * E[i-1] + 1. * E[i]) / (20. * mass);
    double gamma2  = 1. + (17. * E[i-1] + 3. * E[i]) / (20. * mass);
    double gamma3  = 1. + (15. * E[i-1] + 5. * E[i]) / (20. * mass);
    double gamma4  = 1. + (13. * E[i-1] + 7. * E[i]) / (20. * mass);
    double gamma5  = 1. + (11. * E[i-1] + 9. * E[i]) / (20. * mass);
    double gamma6  = 1. + (9. * E[i-1] + 11. * E[i]) / (20. * mass);
    double gamma7  = 1. + (7. * E[i-1] + 13. * E[i]) / (20. * mass);
    double gamma8  = 1. + (5. * E[i-1] + 15. * E[i]) / (20. * mass);
    double gamma9  = 1. + (3. * E[i-1] + 17. * E[i]) / (20. * mass);
    double gamma10 = 1. + (1. * E[i-1] + 19. * E[i]) / (20. * mass);
    return dz *
        (1. / sqrt(1. - 1. / (gamma1 * gamma1)) +
         1. / sqrt(1. - 1. / (gamma2 * gamma2)) +
         1. / sqrt(1. - 1. / (gamma3 * gamma3)) +
         1. / sqrt(1. - 1. / (gamma4 * gamma4)) +
         1. / sqrt(1. - 1. / (gamma5 * gamma5)) +
         1. / sqrt(1. - 1. / (gamma6 * gamma6)) +
         1. / sqrt(1. - 1. / (gamma7 * gamma7)) +
         1. / sqrt(1. - 1. / (gamma8 * gamma8)) +
         1. / sqrt(1. - 1. / (gamma9 * gamma9)) +
         1. / sqrt(1. - 1. / (gamma10 * gamma10))) / (10. * Physics::c);
}

inline
double RFCavity::getdA(const int & i,
323
                       const std::vector<double> & t,
gsell's avatar
gsell committed
324 325 326 327 328
                       const double & dz,
                       const double & frequency,
                       const std::vector<double> & F) const {
    double dt = t[i] - t[i-1];
    return dz / (frequency * frequency * dt * dt) *
329 330
        (frequency * dt * (F[i] * cos(frequency * t[i]) - F[i-1] * cos(frequency * t[i-1])) -
         (F[i] - F[i-1]) * (sin(frequency * t[i]) - sin(frequency * t[i-1])));
gsell's avatar
gsell committed
331 332 333 334
}

inline
double RFCavity::getdB(const int & i,
335
                       const std::vector<double> & t,
gsell's avatar
gsell committed
336 337 338 339
                       const double & dz,
                       const double & frequency,
                       const std::vector<double> & F) const {
    double dt = t[i] - t[i-1];
340 341 342 343 344
    return dz / (frequency * frequency * dt * dt) *
        (frequency * dt * (F[i] * sin(frequency * t[i]) - F[i-1] * sin(frequency * t[i-1])) +
         (F[i] - F[i-1]) * (cos(frequency * t[i]) - cos(frequency * t[i-1])));
}

gsell's avatar
gsell committed
345 346
inline
void RFCavity::setDesignEnergy(const double& ekin, bool)
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432
{
    designEnergy_m = ekin;
}

inline
double RFCavity::getDesignEnergy() const
{
    return designEnergy_m;
}

inline
void RFCavity::dropFieldmaps() {
    fieldmap_m = NULL;
}

inline
void RFCavity::setFieldMapFN(std::string fn) {
    filename_m = fn;
}

inline
std::string RFCavity::getFieldMapFN() const {
    return filename_m;
}

inline
void RFCavity::setAmplitudem(double vPeak) {
    scale_m = vPeak;
}

inline
double RFCavity::getAmplitudem() const {
    return scale_m;
}

inline
void RFCavity::setAmplitudeError(double vPeakError) {
    scaleError_m = vPeakError;
}

inline
double RFCavity::getAmplitudeError() const {
    return scaleError_m;
}

inline
void RFCavity::setFrequency(double freq) {
    frequency_m = freq;
}

inline
void RFCavity::setFrequencym(double freq) {
    frequency_m = freq;
}

inline
double RFCavity::getFrequencym() const {
    return frequency_m;
}

inline
void RFCavity::setPhasem(double phase) {
    phase_m = phase;
}

inline
double RFCavity::getPhasem() const {
    return phase_m;
}

inline
double RFCavity::getPhasem(double t) const {
    return phase_m + t * frequency_m;
}

inline
void RFCavity::setPhaseError(double phaseError) {
    phaseError_m = phaseError;
}

inline
double RFCavity::getPhaseError() const {
    return phaseError_m;
}

inline
433
void RFCavity::setCavityType(std::string /*type*/) {
434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504

}

inline
std::string RFCavity::getCavityType() const {
    return "SW";
}

inline
void RFCavity::setFast(bool fast) {
    fast_m = fast;
}

inline
bool RFCavity::getFast() const {
    return fast_m;
}

inline
void RFCavity::setAutophaseVeto(bool veto) {
    autophaseVeto_m = veto;
}

inline
bool RFCavity::getAutophaseVeto() const {
    return autophaseVeto_m;
}

inline
void RFCavity::setAmplitudeModel(std::shared_ptr<AbstractTimeDependence> amplitude_td) {
  amplitude_td_m = amplitude_td;
}

inline
void RFCavity::setAmplitudeModelName(std::string name) {
    amplitude_name_m=name;
}

inline
std::string RFCavity::getAmplitudeModelName() {
    return amplitude_name_m;
}

inline
void RFCavity::setPhaseModel(std::shared_ptr<AbstractTimeDependence> phase_td) {
  phase_td_m = phase_td;
}

inline
void RFCavity::setPhaseModelName(std::string name) {
    phase_name_m=name;
}

inline
std::string RFCavity::getPhaseModelName() {
    return phase_name_m;
}

inline
void RFCavity::setFrequencyModel(std::shared_ptr<AbstractTimeDependence> frequency_td) {
  frequency_td_m = frequency_td;
}

inline
void RFCavity::setFrequencyModelName(std::string name) {
  frequency_name_m=name;
}

inline
std::string RFCavity::getFrequencyModelName() {
    return frequency_name_m;
gsell's avatar
gsell committed
505 506
}

507 508 509 510 511 512 513 514 515 516 517 518 519

inline
CoordinateSystemTrafo RFCavity::getEdgeToBegin() const
{
    CoordinateSystemTrafo ret(Vector_t(0, 0, startField_m),
                              Quaternion(1, 0, 0, 0));

    return ret;
}

inline
CoordinateSystemTrafo RFCavity::getEdgeToEnd() const
{
kraus's avatar
kraus committed
520
    CoordinateSystemTrafo ret(Vector_t(0, 0, startField_m + length_m),
521 522 523 524 525 526
                              Quaternion(1, 0, 0, 0));

    return ret;
}


527
#endif // CLASSIC_RFCavity_HH