ParametrizedEngine/src/poly2d_base.cc

0001 #include "poly2d_base.h"
0002
0003 using namespace magfieldparam;
0004
0005 /////////////////////////////////////////////////////////////////////////////////
0006 //                                                                             //
0007 //  The "poly2d_term" represent a term of a polynomial of 2 variables.         //
0008 //                                                                             //
0009 /////////////////////////////////////////////////////////////////////////////////
0010
0011 //_______________________________________________________________________________
0012 void poly2d_term::Print(std::ostream &out, bool first_term) {
0013   if (first_term)
0014     out << coeff;
0015   else if (coeff > 0.)
0016     out << " + " << coeff;
0017   else
0018     out << " - " << -coeff;
0019   if (np[0] != 0) {
0020     out << "*r";
0021     if (np[0] != 1)
0022       out << "^" << np[0];
0023   }
0024   if (np[1] != 0) {
0025     out << "*z";
0026     if (np[1] != 1)
0027       out << "^" << np[1];
0028   }
0029 }
0030
0031 /////////////////////////////////////////////////////////////////////////////////
0032 //                                                                             //
0033 //  Base class that represent a polynomial of 2 variables. It isn't supposed   //
0034 //  to be used directly and provides no way of setting coefficients directly.  //
0035 //  Such methods must be defined in derived classes.                           //
0036 //                                                                             //
0037 /////////////////////////////////////////////////////////////////////////////////
0038
0039 //_______________________________________________________________________________
0040 double poly2d_base::rval = 0.;  //Last values of r and z used
0041 double poly2d_base::zval = 0.;
0042
0043 double **poly2d_base::rz_pow = nullptr;  //table with calculated r^n*z^m values
0044 unsigned poly2d_base::NTab = 0;          //rz_pow table size
0045 unsigned poly2d_base::NPwr = 0;          //max power in use by CLASS
0046
0047 bool poly2d_base::rz_set = false;
0048
0049 const double poly2d_base::MIN_COEFF = DBL_EPSILON;  //Thresh. for coeff == 0
0050 std::set<poly2d_base *> poly2d_base::poly2d_base_set;
0051
0052 //_______________________________________________________________________________
0053 poly2d_base::~poly2d_base() {
0054   poly2d_base_set.erase(poly2d_base_set.find(this));
0055   if (!poly2d_base_set.empty()) {  //some objects left
0056     if (max_pwr >= NPwr)
0057       NPwr = GetMaxPow();
0058   } else {
0059     if (rz_pow) {
0060       delete[] rz_pow[0];  //deleting the last instance -> memory cleanup
0061       delete[] rz_pow;
0062     }
0063     rz_pow = nullptr;
0064     rval = zval = 0.;
0065     NPwr = 0;
0066     NTab = 0;
0067     rz_set = false;
0068     //      poly2d_base_set.resize(0);
0069   }
0070 }
0071
0072 //_______________________________________________________________________________
0073 void poly2d_base::SetTabSize(const unsigned N) {
0074   if (N <= NTab)
0075     return;
0076   if (rz_pow) {
0077     delete[] rz_pow[0];
0078     delete[] rz_pow;
0079   }
0080   rz_pow = new double *[N];
0081   unsigned jr, dN = N * (N + 1) / 2;
0082   rz_pow[0] = new double[dN];
0083   memset(rz_pow[0], 0, dN * sizeof(double));
0084   rz_pow[0][0] = 1.;
0085   for (jr = 1, dN = N; jr < N; ++jr, --dN) {
0086     rz_pow[jr] = rz_pow[jr - 1] + dN;
0087   }
0088   rval = zval = 0.;
0089   NTab = N;
0090 }
0091
0092 //_______________________________________________________________________________
0093 void poly2d_base::FillTable(const double r, const double z) {
0094   if (!rz_pow)
0095     return;
0096   unsigned jr, jz;
0097   for (jz = 1; jz <= NPwr; ++jz)
0098     rz_pow[0][jz] = z * rz_pow[0][jz - 1];
0099   for (jr = 1; jr <= NPwr; ++jr) {
0100     for (jz = 0; jz <= (NPwr - jr); ++jz) {
0101       rz_pow[jr][jz] = r * rz_pow[jr - 1][jz];
0102     }
0103   }
0104 }
0105
0106 //_______________________________________________________________________________
0107 int poly2d_base::GetMaxPow() {
0108   int curp, maxp = 0;
0109   std::set<poly2d_base *>::iterator it;
0110
0111   for (it = poly2d_base_set.begin(); it != poly2d_base_set.end(); ++it) {
0112     curp = (*it)->max_pwr;
0113     if (curp > maxp)
0114       maxp = curp;
0115   }
0116   return maxp;
0117 }
0118
0119 //_______________________________________________________________________________
0120 void poly2d_base::AdjustTab() {
0121   NPwr = GetMaxPow();
0122   if (NPwr >= NTab)
0123     SetTabSize(NPwr + 1);
0124 }
0125
0126 //_______________________________________________________________________________
0127 void poly2d_base::PrintTab(std::ostream &out, const std::streamsize prec) {
0128   out << "poly2d_base table size NTab = " << NTab << "\tmax. power NPwr = " << NPwr << std::endl;
0129   if (rz_pow) {
0130     if (NPwr < NTab) {
0131       std::streamsize old_prec = out.precision(), wdt = prec + 7;
0132       out.precision(prec);
0133       out << "Table content:" << std::endl;
0134       unsigned jr, jz;
0135       for (jr = 0; jr <= NPwr; ++jr) {
0136         for (jz = 0; jz <= (NPwr - jr); ++jz) {
0137           out << std::setw(wdt) << std::left << rz_pow[jr][jz];
0138         }
0139         out << "|" << std::endl;
0140       }
0141       out.precision(old_prec);
0142     } else {
0143       out << "\tTable size is not adjusted." << std::endl;
0144     }
0145   } else {
0146     out << "\tTable is not allocated." << std::endl;
0147   }
0148 }
0149
0150 //_______________________________________________________________________________
0151 void poly2d_base::SetPoint(const double r, const double z) {
0152   if (!Count())
0153     return;
0154   if (NPwr >= NTab) {
0155     SetTabSize(NPwr + 1);
0156     FillTable(r, z);
0157   } else if ((r != rval) || (z != zval))
0158     FillTable(r, z);
0159   rz_set = true;
0160 }
0161
0162 //_______________________________________________________________________________
0163 double poly2d_base::Eval() {
0164   double S = 0.;
0165   for (unsigned j = 0; j < data.size(); ++j)
0166     S += data[j].coeff * rz_pow[data[j].np[0]][data[j].np[1]];
0167   return S;
0168 }
0169
0170 //_______________________________________________________________________________
0171 void poly2d_base::Collect() {
0172   if (data.empty())
0173     return;
0174
0175   unsigned j1, j2, rpow, zpow, noff = 0, jend = data.size();
0176   double C;
0177   std::vector<bool> mask(jend, false);
0178   max_pwr = 0;
0179
0180   for (j1 = 0; j1 < jend; ++j1) {
0181     if (mask[j1])
0182       continue;
0183     C = data[j1].coeff;
0184     rpow = data[j1].np[0];
0185     zpow = data[j1].np[1];
0186     for (j2 = j1 + 1; j2 < jend; ++j2) {
0187       if (mask[j2])
0188         continue;
0189       if ((rpow == data[j2].np[0]) && (zpow == data[j2].np[1])) {
0190         C += data[j2].coeff;
0191         mask[j2] = true;
0192         ++noff;
0193       }
0194     }
0195     if (fabs(C) > MIN_COEFF) {
0196       data[j1].coeff = C;
0197       if ((rpow = rpow + zpow) > max_pwr)
0198         max_pwr = rpow;
0199     } else {
0200       mask[j1] = true;
0201       ++noff;
0202     }
0203   }
0204   std::vector<poly2d_term> newdata;
0205   newdata.reserve(jend - noff);
0206   for (j1 = 0; j1 < jend; ++j1) {
0207     if (!(mask[j1]))
0208       newdata.push_back(data[j1]);
0209   }
0210   data.swap(newdata);
0211 }
0212
0213 //_______________________________________________________________________________
0214 void poly2d_base::Print(std::ostream &out, const std::streamsize prec) {
0215   if (data.empty()) {
0216     out << "\"poly2d_base\" object contains no terms." << std::endl;
0217     return;
0218   }
0219   out << data.size() << " terms; max. degree = " << max_pwr << ":" << std::endl;
0220   std::streamsize old_prec = out.precision();
0221   out.precision(prec);
0222   data[0].Print(out);
0223   for (unsigned it = 1; it < data.size(); ++it) {
0224     data[it].Print(out, false);
0225   }
0226   out << std::endl;
0227   out.precision(old_prec);
0228 }
0229
0230 //_______________________________________________________________________________
0231 void poly2d_base::Diff(int nvar) {
0232   //differentiate the polynomial by variable nvar.
0233   //
0234   poly2d_term v3;
0235   std::vector<poly2d_term> newdata;
0236   newdata.reserve(data.size());
0237   unsigned cur_pwr = 0, maxp = 0, oldp = max_pwr;
0238   for (unsigned it = 0; it < data.size(); ++it) {
0239     v3 = data[it];
0240     v3.coeff *= v3.np[nvar];
0241     if (v3.coeff != 0.) {
0242       --v3.np[nvar];
0243       newdata.push_back(v3);
0244       if ((cur_pwr = v3.np[0] + v3.np[1]) > maxp)
0245         maxp = cur_pwr;
0246     }
0247   }
0248   newdata.resize(newdata.size());
0249   max_pwr = maxp;
0250   data.swap(newdata);
0251   if (oldp >= NPwr)
0252     NPwr = GetMaxPow();
0253 }
0254
0255 //_______________________________________________________________________________
0256 void poly2d_base::Int(int nvar) {
0257   //Integrate the polynomial by variable# nvar. Doesn't remove terms
0258   //with zero coefficients; if you suspect they can appear, use Compress()
0259   //after the integration.
0260   //
0261   for (unsigned it = 0; it < data.size(); ++it) {
0262     data[it].coeff /= ++data[it].np[nvar];
0263   }
0264   ++max_pwr;
0265   if (max_pwr > NPwr)
0266     NPwr = GetMaxPow();
0267 }
0268
0269 //_______________________________________________________________________________
0270 void poly2d_base::IncPow(int nvar) {
0271   //Multiply the polynomial by variable# nvar
0272   //
0273   for (unsigned it = 0; it < data.size(); ++it) {
0274     ++data[it].np[nvar];
0275   }
0276   ++max_pwr;
0277   if (max_pwr > NPwr)
0278     NPwr = GetMaxPow();
0279 }
0280
0281 //_______________________________________________________________________________
0282 void poly2d_base::DecPow(int nvar) {
0283   //Divide the polynomial by variable# nvar. Remove terms with zero coefficients
0284   //and also terms where the initial power of nvar is equal zero
0285   //
0286   poly2d_term v3;
0287   std::vector<poly2d_term> newdata;
0288   newdata.reserve(data.size());
0289   unsigned cur_pwr = 0, maxp = 0, oldp = max_pwr;
0290   for (unsigned it = 0; it < data.size(); ++it) {
0291     v3 = data[it];
0292     if ((v3.coeff != 0.) && (v3.np[nvar] > 0)) {
0293       --v3.np[nvar];
0294       newdata.push_back(v3);
0295       if ((cur_pwr = v3.np[0] + v3.np[1]) > maxp)
0296         maxp = cur_pwr;
0297     }
0298   }
0299   newdata.resize(newdata.size());
0300   max_pwr = maxp;
0301   data.swap(newdata);
0302   if (oldp >= NPwr)
0303     NPwr = GetMaxPow();
0304 }
0305
0306 //_______________________________________________________________________________
0307 void poly2d_base::Scale(const double C) {
0308   //Multiply the polynomial by a constant.
0309   //
0310   if (C != 0.) {
0311     for (unsigned it = 0; it < data.size(); ++it) {
0312       data[it].coeff *= C;
0313     }
0314   } else
0315     data.resize(0);
0316 }