CMS 3D CMS Logo

rz_harm_poly.cc
Go to the documentation of this file.
1 #include <typeinfo>
2 #include "rz_harm_poly.h"
3 #include <cstdlib>
4 
5 using namespace magfieldparam;
6 
8 // //
9 // Harmonic homogeneous polynomials in cylindrical system. //
10 // //
12 
13 //_______________________________________________________________________________
14 unsigned rz_harm_poly::Cnt = 0; //Number of the "rz_harm_poly" objects
15 double rz_harm_poly::phival = -11111.; //Last phi value used
16 bool rz_harm_poly::phi_set = false; //TRUE if phi value is set
17 unsigned rz_harm_poly::MaxM = 0; //Max. M among "rz_harm_poly" objects
18 
19 unsigned rz_harm_poly::TASize = 0; //TrigArr size
20 trig_pair *rz_harm_poly::TrigArr = nullptr; //Array with angular data
21 
22 //_______________________________________________________________________________
23 rz_harm_poly::rz_harm_poly(const unsigned N) {
24  //Constructor for rz_harm_poly of length N. The polynomial P(r,z) is normalized
25  //in such a way that dP/dz(r=0,z)=z^(N-1)
26  //
27  unsigned nz = N, nr = 0, nv = 0;
28  poly2d_term v3(1. / N, nr, nz);
29 
30  data = std::vector<poly2d_term>((N + 2) / 2, v3);
31 
32  while (nz >= 2) {
33  nz -= 2;
34  nr += 2;
35  nv += 1;
36  data[nv].coeff = -data[nv - 1].coeff * (nz + 1) * (nz + 2) / (nr * nr);
37  data[nv].np[0] = nr;
38  data[nv].np[1] = nz;
39  }
40  max_pwr = N;
41  if (max_pwr > NPwr) {
42  NPwr = max_pwr;
43  rz_set = false;
44  phi_set = false;
45  }
46  L = N;
47  M = 0;
48  poly2d_base_set.insert(this);
49  ++Cnt;
50 }
51 
52 //_______________________________________________________________________________
54  if (--Cnt) {
55  if (std::abs(M) >= int(MaxM)) { //a number of objects still left
56  M = 0;
57  MaxM = GetMaxM();
58  }
59  } else { //last instance -> memory cleanup
60  if (TrigArr)
61  delete[] TrigArr;
62  TrigArr = nullptr;
63  TASize = 0;
64  MaxM = 0;
65  phival = -11111.;
66  phi_set = false;
67  }
68 }
69 
70 //_______________________________________________________________________________
72  //Return max abs(M) for all rz_harm_poly objects created
73  //
74  int M_cur, M_max = 0;
75  for (auto const &poly : poly2d_base_set) {
76  if (typeid(*poly) == typeid(rz_harm_poly)) {
77  M_cur = std::abs(static_cast<rz_harm_poly const *>(poly)->M);
78  if (M_cur > M_max)
79  M_max = M_cur;
80  }
81  }
82  return M_max;
83 }
84 
85 //_______________________________________________________________________________
86 void rz_harm_poly::SetPhi(const double phi) {
87  //Set value of the angle argument, adjust the TrigArr size if neccessary
88  //and fill TrigArr if the phi value is changed
89  //
90  if (MaxM >= TASize) {
91  SetTrigArrSize(MaxM + 1);
92  FillTrigArr(phi);
93  } else if (phi != phival)
94  FillTrigArr(phi);
95  phival = phi;
96  phi_set = true;
97 }
98 
99 //_______________________________________________________________________________
100 void rz_harm_poly::SetTrigArrSize(const unsigned N) {
101  //Increase TrigArr size if neccessary
102  //
103  if (N <= TASize)
104  return;
105  if (TrigArr)
106  delete[] TrigArr;
107  TrigArr = new trig_pair[N];
108  (*TrigArr) = trig_pair(1., 0.);
109  TASize = N;
110  phi_set = false;
111 }
112 
113 //_______________________________________________________________________________
114 void rz_harm_poly::FillTrigArr(const double phi) {
115  //Fill TrigArr with trig_pair(jp*phi)
116  if (!TrigArr)
117  return;
118  trig_pair tp(phi);
119  TrigArr[1] = tp;
120  for (unsigned jp = 2; jp <= MaxM; ++jp)
121  TrigArr[jp] = TrigArr[jp - 1].Add(tp);
122 }
123 
124 //_______________________________________________________________________________
125 void rz_harm_poly::PrintTrigArr(std::ostream &out, const std::streamsize prec) {
126  out << "TrigArr: TASize = " << TASize << "\tMaxM = " << MaxM << std::endl;
127  if (TrigArr) {
128  if (MaxM < TASize) {
129  unsigned jm;
130  std::streamsize old_prec = out.precision(), wdt = prec + 7;
131  out.precision(prec);
132  out << "M: ";
133  for (jm = 0; jm <= MaxM; ++jm) {
134  out << std::setw(wdt) << std::left << jm;
135  }
136  out << "|\nCos_M: ";
137  for (jm = 0; jm <= MaxM; ++jm) {
138  out << std::setw(wdt) << std::left << TrigArr[jm].CosPhi;
139  }
140  out << "|\nSin_M: ";
141  for (jm = 0; jm <= MaxM; ++jm) {
142  out << std::setw(wdt) << std::left << TrigArr[jm].SinPhi;
143  }
144  out << "|" << std::endl;
145  out.precision(old_prec);
146  } else {
147  out << "\tTrigArr size is not adjusted." << std::endl;
148  }
149  } else {
150  out << "\tTrigArr is not allocated." << std::endl;
151  }
152 }
153 
154 //_______________________________________________________________________________
156  //Return a polynomial with increased M
157  //
158  rz_harm_poly p_out;
159  p_out.data.reserve(2 * L);
160  unsigned it;
161  poly2d_term term;
162 
163  //In 2 passes (for-cycles) to get terms in z-descending order
164  for (it = 0; it < data.size(); ++it) {
165  term = data[it];
166  if (term.np[0]) {
167  term.coeff *= int(term.np[0]) - M;
168  --term.np[0];
169  ++term.np[1];
170  p_out.data.push_back(term);
171  }
172  }
173  for (it = 0; it < data.size(); ++it) {
174  term = data[it];
175  if (term.np[1]) {
176  term.coeff *= -(int)term.np[1];
177  --term.np[1];
178  ++term.np[0];
179  p_out.data.push_back(term);
180  }
181  }
182  p_out.Collect();
183  if (!p_out.data.empty()) {
184  p_out.L = L;
185  p_out.M = M + 1;
186  if (std::abs(p_out.M) > int(MaxM))
187  MaxM = std::abs(p_out.M);
188  p_out.Scale(1. / sqrt(double((L - M) * (L + M + 1))));
189  }
190  return p_out;
191 }
192 
193 //_______________________________________________________________________________
195  //Return a polynomial with decreased M
196  //
197  rz_harm_poly p_out;
198  p_out.data.reserve(2 * L);
199  unsigned it;
200  poly2d_term term;
201 
202  //In 2 passes (for-cycles) to get terms in z-descending order
203  for (it = 0; it < data.size(); ++it) {
204  term = data[it];
205  if (term.np[0]) {
206  term.coeff *= -int(term.np[0]) - M;
207  --term.np[0];
208  ++term.np[1];
209  p_out.data.push_back(term);
210  }
211  }
212  for (it = 0; it < data.size(); ++it) {
213  term = data[it];
214  if (term.np[1]) {
215  term.coeff *= term.np[1];
216  --term.np[1];
217  ++term.np[0];
218  p_out.data.push_back(term);
219  }
220  }
221  p_out.Collect();
222  if (!p_out.data.empty()) {
223  p_out.L = L;
224  p_out.M = M - 1;
225  if (std::abs(p_out.M) > int(MaxM))
226  MaxM = std::abs(p_out.M);
227  p_out.Scale(1. / sqrt(double((L + M) * (L - M + 1))));
228  }
229  return p_out;
230 }
static void FillTrigArr(const double phi)
std::vector< poly2d_term > data
Definition: poly2d_base.h:69
static std::set< poly2d_base * > poly2d_base_set
Definition: poly2d_base.h:60
static void SetPhi(const double phi)
Definition: rz_harm_poly.cc:86
void Scale(const double C)
Definition: poly2d_base.cc:307
static trig_pair * TrigArr
Definition: rz_harm_poly.h:49
T sqrt(T t)
Definition: SSEVec.h:19
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Definition: poly.h:11
#define N
Definition: blowfish.cc:9
static unsigned NPwr
Definition: poly2d_base.h:55
static void PrintTrigArr(std::ostream &out=std::cout, const std::streamsize prec=5)
static void SetTrigArrSize(const unsigned N)