ParametrizedEngine/src/rz_harm_poly.cc

0001 #include <typeinfo>
0002 #include "rz_harm_poly.h"
0003 #include <cstdlib>
0004
0005 using namespace magfieldparam;
0006
0007 /////////////////////////////////////////////////////////////////////////////////
0008 //                                                                             //
0009 //  Harmonic homogeneous polynomials in cylindrical system.                    //
0010 //                                                                             //
0011 /////////////////////////////////////////////////////////////////////////////////
0012
0013 //_______________________________________________________________________________
0014 unsigned rz_harm_poly::Cnt = 0;         //Number of the "rz_harm_poly" objects
0015 double rz_harm_poly::phival = -11111.;  //Last phi value used
0016 bool rz_harm_poly::phi_set = false;     //TRUE if phi value is set
0017 unsigned rz_harm_poly::MaxM = 0;        //Max. M among "rz_harm_poly" objects
0018
0019 unsigned rz_harm_poly::TASize = 0;           //TrigArr size
0020 trig_pair *rz_harm_poly::TrigArr = nullptr;  //Array with angular data
0021
0022 //_______________________________________________________________________________
0023 rz_harm_poly::rz_harm_poly(const unsigned N) {
0024   //Constructor for rz_harm_poly of length N. The polynomial P(r,z) is normalized
0025   //in such a way that dP/dz(r=0,z)=z^(N-1)
0026   //
0027   unsigned nz = N, nr = 0, nv = 0;
0028   poly2d_term v3(1. / N, nr, nz);
0029
0030   data = std::vector<poly2d_term>((N + 2) / 2, v3);
0031
0032   while (nz >= 2) {
0033     nz -= 2;
0034     nr += 2;
0035     nv += 1;
0036     data[nv].coeff = -data[nv - 1].coeff * (nz + 1) * (nz + 2) / (nr * nr);
0037     data[nv].np[0] = nr;
0038     data[nv].np[1] = nz;
0039   }
0040   max_pwr = N;
0041   if (max_pwr > NPwr) {
0042     NPwr = max_pwr;
0043     rz_set = false;
0044     phi_set = false;
0045   }
0046   L = N;
0047   M = 0;
0048   poly2d_base_set.insert(this);
0049   ++Cnt;
0050 }
0051
0052 //_______________________________________________________________________________
0053 rz_harm_poly::~rz_harm_poly() {
0054   if (--Cnt) {
0055     if (std::abs(M) >= int(MaxM)) {  //a number of objects still left
0056       M = 0;
0057       MaxM = GetMaxM();
0058     }
0059   } else {  //last instance -> memory cleanup
0060     if (TrigArr)
0061       delete[] TrigArr;
0062     TrigArr = nullptr;
0063     TASize = 0;
0064     MaxM = 0;
0065     phival = -11111.;
0066     phi_set = false;
0067   }
0068 }
0069
0070 //_______________________________________________________________________________
0071 int rz_harm_poly::GetMaxM() {
0072   //Return max abs(M) for all rz_harm_poly objects created
0073   //
0074   int M_cur, M_max = 0;
0075   for (auto const &poly : poly2d_base_set) {
0076     if (typeid(*poly) == typeid(rz_harm_poly)) {
0077       M_cur = std::abs(static_cast<rz_harm_poly const *>(poly)->M);
0078       if (M_cur > M_max)
0079         M_max = M_cur;
0080     }
0081   }
0082   return M_max;
0083 }
0084
0085 //_______________________________________________________________________________
0086 void rz_harm_poly::SetPhi(const double phi) {
0087   //Set value of the angle argument, adjust the TrigArr size if neccessary
0088   //and fill TrigArr if the phi value is changed
0089   //
0090   if (MaxM >= TASize) {
0091     SetTrigArrSize(MaxM + 1);
0092     FillTrigArr(phi);
0093   } else if (phi != phival)
0094     FillTrigArr(phi);
0095   phival = phi;
0096   phi_set = true;
0097 }
0098
0099 //_______________________________________________________________________________
0100 void rz_harm_poly::SetTrigArrSize(const unsigned N) {
0101   //Increase TrigArr size if neccessary
0102   //
0103   if (N <= TASize)
0104     return;
0105   if (TrigArr)
0106     delete[] TrigArr;
0107   TrigArr = new trig_pair[N];
0108   (*TrigArr) = trig_pair(1., 0.);
0109   TASize = N;
0110   phi_set = false;
0111 }
0112
0113 //_______________________________________________________________________________
0114 void rz_harm_poly::FillTrigArr(const double phi) {
0115   //Fill TrigArr with trig_pair(jp*phi)
0116   if (!TrigArr)
0117     return;
0118   trig_pair tp(phi);
0119   TrigArr[1] = tp;
0120   for (unsigned jp = 2; jp <= MaxM; ++jp)
0121     TrigArr[jp] = TrigArr[jp - 1].Add(tp);
0122 }
0123
0124 //_______________________________________________________________________________
0125 void rz_harm_poly::PrintTrigArr(std::ostream &out, const std::streamsize prec) {
0126   out << "TrigArr: TASize = " << TASize << "\tMaxM = " << MaxM << std::endl;
0127   if (TrigArr) {
0128     if (MaxM < TASize) {
0129       unsigned jm;
0130       std::streamsize old_prec = out.precision(), wdt = prec + 7;
0131       out.precision(prec);
0132       out << "M:     ";
0133       for (jm = 0; jm <= MaxM; ++jm) {
0134         out << std::setw(wdt) << std::left << jm;
0135       }
0136       out << "|\nCos_M: ";
0137       for (jm = 0; jm <= MaxM; ++jm) {
0138         out << std::setw(wdt) << std::left << TrigArr[jm].CosPhi;
0139       }
0140       out << "|\nSin_M: ";
0141       for (jm = 0; jm <= MaxM; ++jm) {
0142         out << std::setw(wdt) << std::left << TrigArr[jm].SinPhi;
0143       }
0144       out << "|" << std::endl;
0145       out.precision(old_prec);
0146     } else {
0147       out << "\tTrigArr size is not adjusted." << std::endl;
0148     }
0149   } else {
0150     out << "\tTrigArr is not allocated." << std::endl;
0151   }
0152 }
0153
0154 //_______________________________________________________________________________
0155 rz_harm_poly rz_harm_poly::LadderUp() {
0156   //Return a polynomial with increased M
0157   //
0158   rz_harm_poly p_out;
0159   p_out.data.reserve(2 * L);
0160   unsigned it;
0161   poly2d_term term;
0162
0163   //In 2 passes (for-cycles) to get terms in z-descending order
0164   for (it = 0; it < data.size(); ++it) {
0165     term = data[it];
0166     if (term.np[0]) {
0167       term.coeff *= int(term.np[0]) - M;
0168       --term.np[0];
0169       ++term.np[1];
0170       p_out.data.push_back(term);
0171     }
0172   }
0173   for (it = 0; it < data.size(); ++it) {
0174     term = data[it];
0175     if (term.np[1]) {
0176       term.coeff *= -(int)term.np[1];
0177       --term.np[1];
0178       ++term.np[0];
0179       p_out.data.push_back(term);
0180     }
0181   }
0182   p_out.Collect();
0183   if (!p_out.data.empty()) {
0184     p_out.L = L;
0185     p_out.M = M + 1;
0186     if (std::abs(p_out.M) > int(MaxM))
0187       MaxM = std::abs(p_out.M);
0188     p_out.Scale(1. / sqrt(double((L - M) * (L + M + 1))));
0189   }
0190   return p_out;
0191 }
0192
0193 //_______________________________________________________________________________
0194 rz_harm_poly rz_harm_poly::LadderDwn() {
0195   //Return a polynomial with decreased M
0196   //
0197   rz_harm_poly p_out;
0198   p_out.data.reserve(2 * L);
0199   unsigned it;
0200   poly2d_term term;
0201
0202   //In 2 passes (for-cycles) to get terms in z-descending order
0203   for (it = 0; it < data.size(); ++it) {
0204     term = data[it];
0205     if (term.np[0]) {
0206       term.coeff *= -int(term.np[0]) - M;
0207       --term.np[0];
0208       ++term.np[1];
0209       p_out.data.push_back(term);
0210     }
0211   }
0212   for (it = 0; it < data.size(); ++it) {
0213     term = data[it];
0214     if (term.np[1]) {
0215       term.coeff *= term.np[1];
0216       --term.np[1];
0217       ++term.np[0];
0218       p_out.data.push_back(term);
0219     }
0220   }
0221   p_out.Collect();
0222   if (!p_out.data.empty()) {
0223     p_out.L = L;
0224     p_out.M = M - 1;
0225     if (std::abs(p_out.M) > int(MaxM))
0226       MaxM = std::abs(p_out.M);
0227     p_out.Scale(1. / sqrt(double((L + M) * (L - M + 1))));
0228   }
0229   return p_out;
0230 }