text: #ifndef rz_harm_poly_h
  #define rz_harm_poly_h
  
  #include <cmath>
  #include "poly2d_base.h"
  
  /////////////////////////////////////////////////////////////////////////////////
  //                                                                             //
  //  Pair (Cos(phi),Sin(Phi)). Intended for internal use by rz_harm_poly.       //
  //                                                                             //
  /////////////////////////////////////////////////////////////////////////////////
  
  //_______________________________________________________________________________
  namespace magfieldparam {
    struct trig_pair {
      double CosPhi;
      double SinPhi;
  
      trig_pair() : CosPhi(1.), SinPhi(0.) {}
      trig_pair(const trig_pair &tp) : CosPhi(tp.CosPhi), SinPhi(tp.SinPhi) {}
      trig_pair(const double C, const double S) : CosPhi(C), SinPhi(S) {}
      trig_pair(const double phi) : CosPhi(cos(phi)), SinPhi(sin(phi)) {}
  
      //Return trig_pair fo angle increased by angle of tp.
      trig_pair Add(const trig_pair &tp) {
        return trig_pair(this->CosPhi * tp.CosPhi - this->SinPhi * tp.SinPhi,
                         this->SinPhi * tp.CosPhi + this->CosPhi * tp.SinPhi);
      }
    };
  
    /////////////////////////////////////////////////////////////////////////////////
    //                                                                             //
    //  Harmonic homogeneous polynomial in cylindrical system.                     //
    //                                                                             //
    /////////////////////////////////////////////////////////////////////////////////
  
    //_______________________________________________________________________________
    class rz_harm_poly : public poly2d_base {
    private:
      unsigned L;
      int M;
  
      static unsigned Cnt;   //Number of the "rz_harm_poly" objects
      static double phival;  //Last phi value used
      static bool phi_set;   //TRUE if phi value is set
      static unsigned MaxM;  //Max. M among "rz_harm_poly" objects
  
      static unsigned TASize;     //TrigArr size
      static trig_pair *TrigArr;  //Array with angular data
  
      static void SetTrigArrSize(const unsigned N);
      static void FillTrigArr(const double phi);
  
      void PrintLM(std::ostream &out = std::cout) {
        out << "L=" << std::setw(3) << std::left << L << ", M=" << std::setw(3) << std::left << M << "; ";
      }
  
    public:
      static int GetMaxM();  //return Max. M for the class
      static unsigned ParentCount() { return poly2d_base::Count(); }
      static unsigned Count() { return Cnt; }
      static void SetPhi(const double phi);
      static void SetPoint(const double r, const double z, const double phi) {
        poly2d_base::SetPoint(r, z);
        SetPhi(phi);
      }
  
      rz_harm_poly() : poly2d_base(), L(0), M(0) { ++Cnt; }
      rz_harm_poly(const poly2d_base &S) : poly2d_base(S), L(0), M(0) { ++Cnt; }
      rz_harm_poly(const rz_harm_poly &S) : poly2d_base(S), L(S.L), M(S.M) { ++Cnt; }
      rz_harm_poly(const unsigned N);
      ~rz_harm_poly() override;
  
      bool IsPhiSet() { return phi_set; }
  
      rz_harm_poly GetDiff(int nvar) {
        rz_harm_poly R(*this);
        R.Diff(nvar);
        return R;
      }
      rz_harm_poly GetInt(int nvar) {
        rz_harm_poly R(*this);
        R.Int(nvar);
        return R;
      }
      rz_harm_poly GetIncPow(int nvar) {
        rz_harm_poly R(*this);
        R.IncPow(nvar);
        return R;
      }
      rz_harm_poly GetDecPow(int nvar) {
        rz_harm_poly R(*this);
        R.DecPow(nvar);
        return R;
      }
  
      rz_harm_poly LadderUp();
      rz_harm_poly LadderDwn();
  
      unsigned GetL() { return L; }
      int GetM() { return M; }
  
      //Next functions return value of angular terms.
      //No check is made, wheither the TrigArr is initialized.
      //User can check if IsPhiSet() == true
      double GetCos() { return TrigArr[M].CosPhi; }
      double GetSin() { return TrigArr[M].SinPhi; }
  
      void CheatL(const unsigned newL) { L = newL; }
      void Print(std::ostream &out = std::cout, const std::streamsize prec = 5) {
        PrintLM(out);
        poly2d_base::Print(out, prec);
      }
  
      static void PrintTrigArr(std::ostream &out = std::cout, const std::streamsize prec = 5);
  
    };  //class rz_harm_poly
  }  // namespace magfieldparam
  
  #endif
  