HarmBasis3DCyl.cc

Go to the documentation of this file.

//                                                                             //
//  HarmBasis3DCyl: set of basis harmonic polynomials in cylindrical CS        //
//                                                                             //

#include "HarmBasis3DCyl.h"

using namespace magfieldparam;

//_______________________________________________________________________________
HarmBasis3DCyl::HarmBasis3DCyl(const unsigned N) {
  //Construct a basis of dimension N
  //
  Dim = N;
  Len = N * (N + 2);

  L_k = new int[Len];
  M_k = new int[Len];

  P_k = new double[Len];
  Br_k = new double[Len];
  Bz_k = new double[Len];
  Bphi_k = new double[Len];

  PtB.reserve(N);
  BrB.reserve(N);
  BzB.reserve(N);
  BphiB.reserve(N);

  rz_harm_poly::IncNPwr(N);  //In order to prevent GetMaxPow() calls
  unsigned M, vLen, k = 0;
  for (unsigned L = 1; L <= N; ++L) {
    vLen = L + 1;
    harm_poly_vec Pt_vec;
    Pt_vec.reserve(vLen);
    harm_poly_vec Br_vec;
    Br_vec.reserve(vLen);
    harm_poly_vec Bz_vec;
    Bz_vec.reserve(vLen);
    harm_poly_vec Bphi_vec;
    Bphi_vec.reserve(vLen);

    Pt_vec.push_back(rz_harm_poly(L));
    Br_vec.push_back(Pt_vec[0].GetDiff(0));
    Bz_vec.push_back(Pt_vec[0].GetDiff(1));
    Bphi_vec.push_back(rz_harm_poly());
    Bphi_vec[0].CheatL(L);

    L_k[k] = L;
    M_k[k] = 0;
    ++k;

    for (M = 1; M <= L; ++M) {
      Pt_vec.push_back(Pt_vec[M - 1].LadderUp());
      Br_vec.push_back(Pt_vec[M].GetDiff(0));
      Bz_vec.push_back(Pt_vec[M].GetDiff(1));
      Bphi_vec.push_back(Pt_vec[M].GetDecPow(0));
      Bphi_vec[M].Scale(M);
      L_k[k] = L;
      M_k[k] = M;
      ++k;
      L_k[k] = L;
      M_k[k] = -M;
      ++k;
    }
    PtB.push_back(Pt_vec);
    BrB.push_back(Br_vec);
    BzB.push_back(Bz_vec);
    BphiB.push_back(Bphi_vec);
  }
}

//_______________________________________________________________________________
HarmBasis3DCyl::~HarmBasis3DCyl() {
  delete[] Bphi_k;
  delete[] Bz_k;
  delete[] Br_k;
  delete[] P_k;
  delete[] M_k;
  delete[] L_k;
}

//_______________________________________________________________________________
void HarmBasis3DCyl::EvalRZ(harm_poly_arr &B, double *val) {
  //Fills the linear array val[Len] with values of basis polynomials.
  //Private function, intended for internal use only.
  //
  unsigned M;
  double V;
  rz_harm_poly *P;
  for (unsigned L = 1, k = 0; L <= Dim; ++L, ++k) {
    (*val) = B[k][0].Eval();
    ++val;
    for (M = 1; M <= L; ++M) {
      P = &(B[k][M]);
      V = P->Eval();
      (*val) = V * P->GetCos();
      ++val;
      (*val) = V * P->GetSin();
      ++val;
    }
  }
}

//_______________________________________________________________________________
void HarmBasis3DCyl::EvalBphi() {
  //Fills the array Bphi_k[Len] with values of phi-basis polynomials.
  //
  unsigned M;
  double V;
  double *val = Bphi_k;
  rz_harm_poly *P;
  for (unsigned L = 1, k = 0; L <= Dim; ++L, ++k) {
    (*val) = 0.;
    ++val;
    for (M = 1; M <= L; ++M) {
      P = &(BphiB[k][M]);
      V = P->Eval();
      (*val) = -V * P->GetSin();
      ++val;
      (*val) = V * P->GetCos();
      ++val;
    }
  }
}

//_______________________________________________________________________________
double HarmBasis3DCyl::GetVal(double *coeff, double *basis) {
  //return value of the expansion with coefficients coeff[Len] for the basis
  //Private function, intended for internal use only.
  //
  double S = 0.;
  for (unsigned k = 0; k < Len; ++k)
    S += coeff[k] * basis[k];
  return S;
}

//_______________________________________________________________________________
void HarmBasis3DCyl::Print(harm_poly_arr &B, std::ostream &out) {
  unsigned jL, jM, wdt = 60;
  char fc1 = '-', fc0 = out.fill(fc1);
  for (jL = 0; jL < B.size(); ++jL) {
    out << std::setw(wdt) << fc1 << std::endl;
    out << "Basis subset " << jL + 1 << std::endl;
    out << std::setw(wdt) << fc1 << std::endl;
    for (jM = 0; jM < B[jL].size(); ++jM) {
      B[jL][jM].Print(out);
    }
  }
  out.fill(fc0);
}

//_______________________________________________________________________________
void HarmBasis3DCyl::Print(std::ostream &out) {
  out << "BASIS POLYNOMIALS FOR THE POTENTIAL:\n" << std::endl;
  PrintPtB(out);
  out << "\nBASIS POLYNOMIALS FOR R-COMPONENT   OF THE FIELD:\n" << std::endl;
  PrintBrB(out);
  out << "\nBASIS POLYNOMIALS FOR Z-COMPONENT   OF THE FIELD:\n" << std::endl;
  PrintBzB(out);
  out << "\nBASIS POLYNOMIALS FOR PHI-COMPONENT OF THE FIELD:\n" << std::endl;
  PrintBphiB(out);
}