ParametrizedEngine/src/HarmBasis3DCyl.cc

0001 /////////////////////////////////////////////////////////////////////////////////
0002 //                                                                             //
0003 //  HarmBasis3DCyl: set of basis harmonic polynomials in cylindrical CS        //
0004 //                                                                             //
0005 /////////////////////////////////////////////////////////////////////////////////
0006
0007 #include "HarmBasis3DCyl.h"
0008
0009 using namespace magfieldparam;
0010
0011 //_______________________________________________________________________________
0012 HarmBasis3DCyl::HarmBasis3DCyl(const unsigned N) {
0013   //Construct a basis of dimension N
0014   //
0015   Dim = N;
0016   Len = N * (N + 2);
0017
0018   L_k = new int[Len];
0019   M_k = new int[Len];
0020
0021   P_k = new double[Len];
0022   Br_k = new double[Len];
0023   Bz_k = new double[Len];
0024   Bphi_k = new double[Len];
0025
0026   PtB.reserve(N);
0027   BrB.reserve(N);
0028   BzB.reserve(N);
0029   BphiB.reserve(N);
0030
0031   rz_harm_poly::IncNPwr(N);  //In order to prevent GetMaxPow() calls
0032   unsigned M, vLen, k = 0;
0033   for (unsigned L = 1; L <= N; ++L) {
0034     vLen = L + 1;
0035     harm_poly_vec Pt_vec;
0036     Pt_vec.reserve(vLen);
0037     harm_poly_vec Br_vec;
0038     Br_vec.reserve(vLen);
0039     harm_poly_vec Bz_vec;
0040     Bz_vec.reserve(vLen);
0041     harm_poly_vec Bphi_vec;
0042     Bphi_vec.reserve(vLen);
0043
0044     Pt_vec.push_back(rz_harm_poly(L));
0045     Br_vec.push_back(Pt_vec[0].GetDiff(0));
0046     Bz_vec.push_back(Pt_vec[0].GetDiff(1));
0047     Bphi_vec.push_back(rz_harm_poly());
0048     Bphi_vec[0].CheatL(L);
0049
0050     L_k[k] = L;
0051     M_k[k] = 0;
0052     ++k;
0053
0054     for (M = 1; M <= L; ++M) {
0055       Pt_vec.push_back(Pt_vec[M - 1].LadderUp());
0056       Br_vec.push_back(Pt_vec[M].GetDiff(0));
0057       Bz_vec.push_back(Pt_vec[M].GetDiff(1));
0058       Bphi_vec.push_back(Pt_vec[M].GetDecPow(0));
0059       Bphi_vec[M].Scale(M);
0060       L_k[k] = L;
0061       M_k[k] = M;
0062       ++k;
0063       L_k[k] = L;
0064       M_k[k] = -M;
0065       ++k;
0066     }
0067     PtB.push_back(Pt_vec);
0068     BrB.push_back(Br_vec);
0069     BzB.push_back(Bz_vec);
0070     BphiB.push_back(Bphi_vec);
0071   }
0072 }
0073
0074 //_______________________________________________________________________________
0075 HarmBasis3DCyl::~HarmBasis3DCyl() {
0076   delete[] Bphi_k;
0077   delete[] Bz_k;
0078   delete[] Br_k;
0079   delete[] P_k;
0080   delete[] M_k;
0081   delete[] L_k;
0082 }
0083
0084 //_______________________________________________________________________________
0085 void HarmBasis3DCyl::EvalRZ(harm_poly_arr &B, double *val) {
0086   //Fills the linear array val[Len] with values of basis polynomials.
0087   //Private function, intended for internal use only.
0088   //
0089   unsigned M;
0090   double V;
0091   rz_harm_poly *P;
0092   for (unsigned L = 1, k = 0; L <= Dim; ++L, ++k) {
0093     (*val) = B[k][0].Eval();
0094     ++val;
0095     for (M = 1; M <= L; ++M) {
0096       P = &(B[k][M]);
0097       V = P->Eval();
0098       (*val) = V * P->GetCos();
0099       ++val;
0100       (*val) = V * P->GetSin();
0101       ++val;
0102     }
0103   }
0104 }
0105
0106 //_______________________________________________________________________________
0107 void HarmBasis3DCyl::EvalBphi() {
0108   //Fills the array Bphi_k[Len] with values of phi-basis polynomials.
0109   //
0110   unsigned M;
0111   double V;
0112   double *val = Bphi_k;
0113   rz_harm_poly *P;
0114   for (unsigned L = 1, k = 0; L <= Dim; ++L, ++k) {
0115     (*val) = 0.;
0116     ++val;
0117     for (M = 1; M <= L; ++M) {
0118       P = &(BphiB[k][M]);
0119       V = P->Eval();
0120       (*val) = -V * P->GetSin();
0121       ++val;
0122       (*val) = V * P->GetCos();
0123       ++val;
0124     }
0125   }
0126 }
0127
0128 //_______________________________________________________________________________
0129 double HarmBasis3DCyl::GetVal(double *coeff, double *basis) {
0130   //return value of the expansion with coefficients coeff[Len] for the basis
0131   //Private function, intended for internal use only.
0132   //
0133   double S = 0.;
0134   for (unsigned k = 0; k < Len; ++k)
0135     S += coeff[k] * basis[k];
0136   return S;
0137 }
0138
0139 //_______________________________________________________________________________
0140 void HarmBasis3DCyl::Print(harm_poly_arr &B, std::ostream &out) {
0141   unsigned jL, jM, wdt = 60;
0142   char fc1 = '-', fc0 = out.fill(fc1);
0143   for (jL = 0; jL < B.size(); ++jL) {
0144     out << std::setw(wdt) << fc1 << std::endl;
0145     out << "Basis subset " << jL + 1 << std::endl;
0146     out << std::setw(wdt) << fc1 << std::endl;
0147     for (jM = 0; jM < B[jL].size(); ++jM) {
0148       B[jL][jM].Print(out);
0149     }
0150   }
0151   out.fill(fc0);
0152 }
0153
0154 //_______________________________________________________________________________
0155 void HarmBasis3DCyl::Print(std::ostream &out) {
0156   out << "BASIS POLYNOMIALS FOR THE POTENTIAL:\n" << std::endl;
0157   PrintPtB(out);
0158   out << "\nBASIS POLYNOMIALS FOR R-COMPONENT   OF THE FIELD:\n" << std::endl;
0159   PrintBrB(out);
0160   out << "\nBASIS POLYNOMIALS FOR Z-COMPONENT   OF THE FIELD:\n" << std::endl;
0161   PrintBzB(out);
0162   out << "\nBASIS POLYNOMIALS FOR PHI-COMPONENT OF THE FIELD:\n" << std::endl;
0163   PrintBphiB(out);
0164 }