rz_harm_poly.cc

Go to the documentation of this file.

#include <typeinfo>
#include "rz_harm_poly.h"
#include <cstdlib>

using namespace magfieldparam;

//                                                                             //
//  Harmonic homogeneous polynomials in cylindrical system.                    //
//                                                                             //

//_______________________________________________________________________________
unsigned rz_harm_poly::Cnt    = 0;      //Number of the "rz_harm_poly" objects
double   rz_harm_poly::phival = -11111.;//Last phi value used
bool     rz_harm_poly::phi_set = false; //TRUE if phi value is set
unsigned rz_harm_poly::MaxM   = 0;      //Max. M among "rz_harm_poly" objects

unsigned   rz_harm_poly::TASize  = 0;   //TrigArr size
trig_pair *rz_harm_poly::TrigArr = nullptr;   //Array with angular data

//_______________________________________________________________________________
rz_harm_poly::rz_harm_poly(const unsigned N)
{
//Constructor for rz_harm_poly of length N. The polynomial P(r,z) is normalized
//in such a way that dP/dz(r=0,z)=z^(N-1)
//
   unsigned nz = N, nr = 0, nv = 0;
   poly2d_term v3(1./N, nr, nz);

   data = std::vector<poly2d_term>((N + 2) / 2, v3);

   while (nz >= 2) {
      nz -= 2;
      nr += 2;
      nv += 1;
      data[nv].coeff = -data[nv-1].coeff*(nz+1)*(nz+2)/(nr*nr);
      data[nv].np[0] = nr;
      data[nv].np[1] = nz;
   }
   max_pwr = N;
   if (max_pwr > NPwr) {
      NPwr = max_pwr;
      rz_set = false;
      phi_set = false;
   }
   L = N;
   M = 0;
   poly2d_base_set.insert(this);
   ++Cnt;
}

//_______________________________________________________________________________
rz_harm_poly::~rz_harm_poly()
{
   if (--Cnt) {
     if (std::abs(M) >= int(MaxM)) { //a number of objects still left
         M = 0;
         MaxM = GetMaxM();
      }
   } else { //last instance -> memory cleanup
      if (TrigArr) delete [] TrigArr;
      TrigArr = nullptr;
      TASize = 0;
      MaxM = 0;
      phival = -11111.;
      phi_set = false;
   }
}

//_______________________________________________________________________________
int rz_harm_poly::GetMaxM()
{
//Return max abs(M) for all rz_harm_poly objects created
//
   int M_cur, M_max = 0;
   for (auto const& poly:  poly2d_base_set) {
      if (typeid(*poly) == typeid(rz_harm_poly)) {
         M_cur = std::abs(static_cast<rz_harm_poly const*>(poly)->M);
         if (M_cur > M_max) M_max = M_cur;
      }
   }
   return M_max;
}

//_______________________________________________________________________________
void rz_harm_poly::SetPhi(const double phi)
{
//Set value of the angle argument, adjust the TrigArr size if neccessary
//and fill TrigArr if the phi value is changed
//
   if (MaxM >= TASize) { SetTrigArrSize(MaxM+1); FillTrigArr(phi);}
   else if (phi != phival) FillTrigArr(phi);
   phival = phi;
   phi_set = true;
}

//_______________________________________________________________________________
void rz_harm_poly::SetTrigArrSize(const unsigned N)
{
//Increase TrigArr size if neccessary
//
   if (N <= TASize) return;
   if (TrigArr) delete [] TrigArr;
   TrigArr = new trig_pair [N];
   (*TrigArr) = trig_pair(1., 0.);
   TASize = N;
   phi_set = false;
}

//_______________________________________________________________________________
void rz_harm_poly::FillTrigArr(const double phi)
{
//Fill TrigArr with trig_pair(jp*phi)
   if (!TrigArr) return;
   trig_pair tp(phi);
   TrigArr[1] = tp;
   for (unsigned jp = 2; jp <= MaxM; ++jp) TrigArr[jp] = TrigArr[jp-1].Add(tp);
}

//_______________________________________________________________________________
void rz_harm_poly::PrintTrigArr(std::ostream &out, const std::streamsize prec)
{
   out << "TrigArr: TASize = " << TASize
       << "\tMaxM = " << MaxM << std::endl;
   if (TrigArr) {
      if (MaxM < TASize) {
         unsigned jm;
         std::streamsize old_prec = out.precision(), wdt = prec+7;
         out.precision(prec);
         out << "M:     ";
         for (jm = 0; jm <= MaxM; ++jm) {
            out << std::setw(wdt) << std::left << jm;
         }
         out << "|\nCos_M: ";
         for (jm = 0; jm <= MaxM; ++jm) {
            out << std::setw(wdt) << std::left << TrigArr[jm].CosPhi;
         }
         out << "|\nSin_M: ";
         for (jm = 0; jm <= MaxM; ++jm) {
            out << std::setw(wdt) << std::left << TrigArr[jm].SinPhi;
         }
         out << "|" << std::endl;
         out.precision(old_prec);
      } else {
         out << "\tTrigArr size is not adjusted." << std::endl;
      }
   } else {
      out << "\tTrigArr is not allocated." << std::endl;
   }
}

//_______________________________________________________________________________
rz_harm_poly rz_harm_poly::LadderUp()
{
//Return a polynomial with increased M
//
   rz_harm_poly p_out; p_out.data.reserve(2*L);
   unsigned it;
   poly2d_term term;

   //In 2 passes (for-cycles) to get terms in z-descending order
   for(it = 0; it < data.size(); ++it) {
      term = data[it];
      if (term.np[0]) {
         term.coeff *= int(term.np[0]) - M;
         --term.np[0];
         ++term.np[1];
         p_out.data.push_back(term);
      }
   }
   for(it = 0; it < data.size(); ++it) {
      term = data[it];
      if (term.np[1]) {
         term.coeff *= -(int)term.np[1];
         --term.np[1];
         ++term.np[0];
         p_out.data.push_back(term);
      }
   }
   p_out.Collect();
   if (!p_out.data.empty()) {
      p_out.L = L;
      p_out.M = M+1;
      if (std::abs(p_out.M) > int(MaxM)) MaxM = std::abs(p_out.M);
      p_out.Scale(1./sqrt(double((L-M)*(L+M+1))));
   }
   return p_out;
}

//_______________________________________________________________________________
rz_harm_poly rz_harm_poly::LadderDwn()
{
//Return a polynomial with decreased M
//
   rz_harm_poly p_out; p_out.data.reserve(2*L);
   unsigned it;
   poly2d_term term;

   //In 2 passes (for-cycles) to get terms in z-descending order
   for(it = 0; it < data.size(); ++it) {
      term = data[it];
      if (term.np[0]) {
         term.coeff *= -int(term.np[0]) - M;
         --term.np[0];
         ++term.np[1];
         p_out.data.push_back(term);
      }
   }
   for(it = 0; it < data.size(); ++it) {
      term = data[it];
      if (term.np[1]) {
         term.coeff *= term.np[1];
         --term.np[1];
         ++term.np[0];
         p_out.data.push_back(term);
      }
   }
   p_out.Collect();
   if (!p_out.data.empty()) {
      p_out.L = L;
      p_out.M = M-1;
      if (std::abs(p_out.M) > int(MaxM)) MaxM = std::abs(p_out.M);
      p_out.Scale(1./sqrt(double((L+M)*(L-M+1))));
   }
   return p_out;
}