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;
}