sPhenix code displayed by LXR master/​coresoftware/​calibrations/​tpc/​generator/​Rossegger.cc	Source navigation Diff markup Identifier search General search
	Version: master Revert   Change

File indexing completed on 2025-08-05 08:15:40

 #include "Rossegger.h"
 
 #include <TFile.h>
 #include <TMath.h>
 #include <TTree.h>
 
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wshadow"
 #include <boost/math/special_functions.hpp>  //covers all the special functions.
 #pragma GCC diagnostic pop
 
 #include <boost/format.hpp>
 
 #include <algorithm>  // for max
 #include <cmath>
 #include <cstdlib>  // for exit, abs
 #include <fstream>
 #include <iostream>
 #include <sstream>
 #include <string>
 
 // the limu and kimu terms, that i need to think about a little while longer...
 extern "C"
 {
   void dkia_(int *IFAC, double *X, double *A, double *DKI, double *DKID, int *IERRO);
   void dlia_(int *IFAC, double *X, double *A, double *DLI, double *DLID, int *IERRO);
 }
 //
 
 // Bessel Function J_n(x):
 #define jn(order, x) boost::math::cyl_bessel_j(order, x)
 // Bessel (Neumann) Function Y_n(x):
 #define yn(order, x) boost::math::cyl_neumann(order, x)
 // Modified Bessel Function of the first kind I_n(x):
 #define in(order, x) boost::math::cyl_bessel_i(order, x)
 // Modified Bessel Function of the second kind K_n(x):
 #define kn(order, x) boost::math::cyl_bessel_k(order, x)
 #define limu(im_order, x) Rossegger::Limu(im_order, x)
 #define kimu(im_order, x) Rossegger::Kimu(im_order, x)
 
 /*
   This is a modified/renamed copy of Carlos and Tom's "Spacecharge" class, modified to use boost instead of fortran routines, and with phi terms added.
  */
 
 Rossegger::Rossegger(double InnerRadius, double OuterRadius, double Rdo_Z, double precision)
 {
   a = InnerRadius;
   b = OuterRadius;
   L = Rdo_Z;
 
   epsilon = precision;
 
   verbosity = 0;
   pi = M_PI;
 
   PrecalcFreeConstants();
 
   // load the greens functions:
   std::string zeroesfilename = boost::str(boost::format("rosseger_zeroes_eps%1.0E_a%2.2f_b%2.2f_L%2.2f.root") % epsilon % a % b % L);
   TFile *fileptr = TFile::Open(zeroesfilename.c_str(), "READ");
   if (!fileptr)
   {  // generate the lookuptable
     FindMunk(epsilon);
     FindBetamn(epsilon);
     SaveZeroes(zeroesfilename);
   }
   else
   {  // load it from a file
     fileptr->Close();
     LoadZeroes(zeroesfilename);
   }
 
   PrecalcDerivedConstants();
 
   std::cout << "Rossegger object initialized as follows:" << std::endl;
   std::cout << "  Inner Radius = " << a << " cm." << std::endl;
   std::cout << "  Outer Radius = " << b << " cm." << std::endl;
   std::cout << "  Half  Length = " << L << " cm." << std::endl;
 
   if (!CheckZeroes(0.01))
   {
     std::cout << "CheckZeroes(0.01) returned false, exiting" << std::endl;
     exit(1);
   }
   return;
 }
 
 double Rossegger::FindNextZero(double xstart, double localepsilon, int order, double (Rossegger::*func)(int, double))
 {
   double previous = (this->*func)(order, xstart);
   double x = xstart + localepsilon;
   double value = previous;
 
   while (!((value == 0) || (value < 0 && previous > 0) || (value > 0 && previous < 0)))
   {
     //  Rossegger equation 5.12
     value = (this->*func)(order, x);
     if (value == 0)
     {
       std::cout << "hit it exactly!  Go buy a lottery ticket!" << std::endl;
     }
     if ((value == 0) || (value < 0 && previous > 0) || (value > 0 && previous < 0))
     {
       // when we go from one sign to the other, we have bracketed the zero
       // the following is mathematically equivalent to finding the delta
       // between the x of our previous value and the point where we cross the axis
       // then returning x0=x_old+delta
       double slope = (value - previous) / localepsilon;
       double intercept = value - slope * x;
       double x0 = -intercept / slope;
       if (verbosity > 1)
       {
         std::cout << " " << x0 << "," << std::endl;
       }
       double n0 = (this->*func)(order, x - localepsilon);
       double n1 = (this->*func)(order, x + localepsilon);
       if ((n0 < 0 && n1 < 0) || (n0 > 0 && n1 > 0))
       {
         std::cout << "neighbors on both sides have the same sign.  Check your function resolution!" << std::endl;
       }
 
       return x0;
     }
     previous = value;
     x += localepsilon;
   }
   std::cout << "logic break!" << std::endl;
   exit(1);
 }
 
 void Rossegger::FindBetamn(double localepsilon)
 {
   std::cout << "Now filling the Beta[m][n] Array..." << std::endl;
   if (verbosity > 5)
   {
     std::cout << "numberOfOrders= " << NumberOfOrders << std::endl;
   }
   for (int m = 0; m < NumberOfOrders; m++)
   {
     if (verbosity)
     {
       std::cout << "Filling Beta[" << m << "][n]..." << std::endl;
     }
 
     double x = localepsilon;
     for (int n = 0; n < NumberOfOrders; n++)
     {  //  !!!  Off by one from Rossegger convention  !!!
       x = FindNextZero(x, localepsilon, m, &Rossegger::Rmn_for_zeroes);
       Betamn[m][n] = x / b;
       x += localepsilon;
     }
   }
 
   //  Now fill in the N2mn array...
   for (int m = 0; m < NumberOfOrders; m++)
   {
     for (int n = 0; n < NumberOfOrders; n++)
     {
       //  Rossegger Equation 5.17
       //  N^2_mn = 2/(pi * beta)^2 [ {Jm(beta a)/Jm(beta b)}^2 - 1 ]
       N2mn[m][n] = 2 / (pi * pi * Betamn[m][n] * Betamn[m][n]);
       // N2mn[m][n] *= (jn(m,Betamn[m][n]*a)*jn(m,Betamn[m][n]*a)/jn(m,Betamn[m][n]*b)/jn(m,Betamn[m][n]*b) - 1.0);
       double jna_over_jnb = jn(m, Betamn[m][n] * a) / jn(m, Betamn[m][n] * b);
       N2mn[m][n] *= (jna_over_jnb * jna_over_jnb - 1.0);
       // rcc note!  in eq 5.17, N2nm is set with betamn[m][n], but from context that looks to be a typo.  The order is mn everywhere else
       if (verbosity > 1)
       {
         std::cout << "m: " << m << " n: " << n << " N2[m][n]: " << N2mn[m][n];
       }
       double step = 0.01;
       if (verbosity > 1)
       {
         double integral = 0.0;
         for (double r = a; r < b; r += step)
         {
           double rmnval = Rmn(m, n, r);
 
           integral += rmnval * rmnval * r * step;  // Rmn(m,n,r)*Rmn(m,n,r)*r*step;
         }
         std::cout << " Int: " << integral << std::endl;
       }
       // N2mn[m][n] = integral;
     }
   }
 
   std::cout << "Done." << std::endl;
 }
 
 void Rossegger::FindMunk(double localepsilon)
 {
   std::cout << "Now filling the Mu[n][k] Array..." << std::endl;
   if (verbosity > 5)
   {
     std::cout << "numberOfOrders= " << NumberOfOrders << std::endl;
   }
   // We're looking for the zeroes of Rossegger eqn. 5.46:
   // R_nk(mu_nk;a,b)=Limu(Beta_n*a)Kimu(Beta_n*b)-Kimu(Beta_n*a)Limu(Beta_n*b)=0
   // since a and b are fixed, R_nk is a function solely of mu_nk and n.
   // for each 'n' we wish to find the a set of k mu_n's that result in R_nk=0
 
   // could add an option here to load the munks from file if the dimensions match.
 
   for (int n = 0; n < NumberOfOrders; n++)  //  !!!  Off by one from Rossegger convention  !!!
   {
     if (verbosity)
     {
       std::cout << "Filling Mu[" << n << "][k]..." << std::endl;
     }
     double x = localepsilon;
     for (int k = 0; k < NumberOfOrders; k++)
     {
       x = FindNextZero(x, localepsilon, n, &Rossegger::Rnk_for_zeroes);
       Munk[n][k] = x;
       x += localepsilon;
       if (verbosity > 0)
       {
         std::cout << boost::str(boost::format("Mu[%d][%d]=%E") % n % k % Munk[n][k]) << std::endl;
         std::cout << boost::str(boost::format("adjacent values are Rnk[mu-localepsilon]=%E\tRnk[mu+localepsilon]=%E") % (Rnk_for_zeroes(n, x - localepsilon)) % (Rnk_for_zeroes(n, x + localepsilon))) << std::endl;
         if (verbosity > 100)
         {
           std::cout << "values of argument to limu and kimu are " << ((n + 1) * pi / L * a)
                     << " and " << ((n + 1) * pi / L * b) << std::endl;
         }
       }
     }
   }
 
   //  Now fill in the N2nk array...
   for (int n = 0; n < NumberOfOrders; n++)
   {
     for (int k = 0; k < NumberOfOrders; k++)
     {
       //  Rossegger Equation 5.48
       //  Integral of R_nk(r)*R_ns(r) dr/r= delta_ks*N2nk
       //  note that unlike N2mn, there is no convenient shortcut here.
       double integral = 0.0;
       double step = 0.001;
 
       for (double r = a; r < b; r += step)
       {
         double rnkval = Rnk_(n, k, r);  // must used un-optimized.  We don't have those values yet...
 
         integral += rnkval * rnkval / r * step;
       }
       if (verbosity > 1)
       {
         std::cout << " Int: " << integral << std::endl;
       }
       N2nk[n][k] = integral;
       if (verbosity > 1)
       {
         std::cout << "n: " << n << " k: " << k << " N2nk[n][k]: " << N2nk[n][k];
       }
     }
   }
 
   std::cout << "Done." << std::endl;
   return;
 }
 
 bool Rossegger::CheckZeroes(double localepsilon)
 {
   // confirm that the tabulated zeroes are all zeroes of their respective functions:
   double result;
   for (int m = 0; m < NumberOfOrders; m++)
   {
     for (int n = 0; n < NumberOfOrders; n++)
     {  //  !!!  Off by one from Rossegger convention  !!!
       result = Rmn_for_zeroes(m, Betamn[m][n] * b);
       if (abs(result) > localepsilon)
       {
         std::cout << boost::str(boost::format("(m=%d,n=%d) Jm(x)Ym(lx)-Jm(lx)Ym(x) = %f for x=b*%f") % m % n % result % Betamn[m][n]) << std::endl;
         return false;
       }
     }
   }
 
   // R_nk(mu_nk;a,b)=Limu(Beta_n*a)Kimu(Beta_n*b)-Kimu(Beta_n*a)Limu(Beta_n*b)=0
   for (int n = 0; n < NumberOfOrders; n++)
   {
     for (int k = 0; k < NumberOfOrders; k++)
     {  //  !!!  Off by one from Rossegger convention  !!!
       result = Rnk_for_zeroes(n, Munk[n][k]);
       if (abs(result) > localepsilon * 100)
       {
         std::cout << boost::str(boost::format("(n=%d,k=%d) limu(npi*a/L)kimu(npi*b/L)-kimu(npi*a/L)kimu(npi*b/L) = %f (>eps*100) for mu=%f") % n % k % result % Munk[n][k]) << std::endl;
         return false;
       }
     }
   }
 
   return true;
 }
 
 void Rossegger::PrecalcFreeConstants()
 {  // Routine used to fill the arrays of other values used repeatedly
   // these constants depend only on the geometry of the detector
   std::cout << "Precalcing " << (3 * NumberOfOrders + 5 * NumberOfOrders * NumberOfOrders)
             << "geometric constants" << std::endl;
   for (int n = 0; n < NumberOfOrders; n++)
   {
     BetaN[n] = (n + 1) * pi / L;  //  BetaN=(n+1)*pi/L as used in eg 5.32, .46
     BetaN_a[n] = BetaN[n] * a;    //  BetaN*a as used in eg 5.32, .46
     BetaN_b[n] = BetaN[n] * b;    //  BetaN*b as used in eg 5.32, .46
     for (int m = 0; m < NumberOfOrders; m++)
     {
       km_BetaN_a[m][n] = kn(m, BetaN_a[n]);                                                                                                                      // kn(m,BetaN*a) as used in Rossegger 5.32
       im_BetaN_a[m][n] = in(m, BetaN_a[n]);                                                                                                                      // in(m,BetaN*a) as used in Rossegger 5.32
       km_BetaN_b[m][n] = kn(m, BetaN_b[n]);                                                                                                                      // kn(m,BetaN*b) as used in Rossegger 5.33
       im_BetaN_b[m][n] = in(m, BetaN_b[n]);                                                                                                                      // in(m,BetaN*b) as used in Rossegger 5.33
       bessel_denominator[m][n] = TMath::BesselI(m, BetaN_a[n]) * TMath::BesselK(m, BetaN_b[n]) - TMath::BesselI(m, BetaN_b[n]) * TMath::BesselK(m, BetaN_a[n]);  // TMath::BesselI(m,BetaN*a)*TMath::BesselK(m,BetaN*b)-TMath::BesselI(m,BetaN*b)*TMath::BesselK(m,BetaN*a) as in Rossegger 5.65
     }
   }
   return;
 }
 void Rossegger::PrecalcDerivedConstants()
 {  // Routine used to fill the arrays of other values used repeatedly
    // these constants depend on geometry and the zeroes of special functions
   std::cout << "Precalcing " << (6 * NumberOfOrders * NumberOfOrders)
             << " geometric constants" << std::endl;
 
   for (int n = 0; n < NumberOfOrders; n++)
   {
     for (int m = 0; m < NumberOfOrders; m++)
     {
       ym_Betamn_a[m][n] = yn(m, Betamn[m][n] * a);   // yn(m,Betamn[m][n]*a) as used in Rossegger 5.11
       jm_Betamn_a[m][n] = jn(m, Betamn[m][n] * a);   // jn(m,Betamn[m][n]*a) as used in Rossegger 5.11
       sinh_Betamn_L[m][n] = sinh(Betamn[m][n] * L);  // sinh(Betamn[m][n]*L)  as in Rossegger 5.64
     }
     for (int k = 0; k < NumberOfOrders; k++)
     {
       liMunk_BetaN_a[n][k] = limu(Munk[n][k], BetaN_a[n]);  // limu(Munk[n][k],BetaN*a) as used in Rossegger 5.45
       kiMunk_BetaN_a[n][k] = kimu(Munk[n][k], BetaN_a[n]);  // kimu(Munk[n][k],BetaN*a) as used in Rossegger 5.45
       sinh_pi_Munk[n][k] = sinh(pi * Munk[n][k]);           // sinh(pi*Munk[n][k]) as in Rossegger 5.66
     }
   }
   return;
 }
 
 double Rossegger::Limu(double mu, double x)
 {
   // defined in Rossegger eqn 5.44, also a canonical 'satisfactory companion' to Kimu.
   // could use Griddit?
   //   Rossegger Equation 5.45
   //        Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN r) - Kimu_nk(BetaN a) Limu_nk (BetaN r)
 
   int IFAC = 1;
   double A = mu;
   double DLI = 0;
   double DERR = 0;
   int IERRO = 0;
 
   double X = x;
   dlia_(&IFAC, &X, &A, &DLI, &DERR, &IERRO);
   return DLI;
 }
 double Rossegger::Kimu(double mu, double x)
 {
   int IFAC = 1;
   double A = mu;
   double DKI = 0;
   double DERR = 0;
   int IERRO = 0;
 
   double X = x;
   dkia_(&IFAC, &X, &A, &DKI, &DERR, &IERRO);
   return DKI;
 }
 
 double Rossegger::Rmn_for_zeroes(int m, double x)
 {
   double lx = a * x / b;
   //  Rossegger Equation 5.12:
 
   return jn(m, x) * yn(m, lx) - jn(m, lx) * yn(m, x);
 }
 
 double Rossegger::Rmn(int m, int n, double r)
 {
   if (verbosity > 100)
   {
     std::cout << "Determine Rmn(" << m << "," << n << "," << r << ") = ";
   }
 
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rmn(" << m << "," << n << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   //  Calculate the function using C-libraries from boost
   //  Rossegger Equation 5.11:
   //         Rmn(r) = Ym(Beta_mn a)*Jm(Beta_mn r) - Jm(Beta_mn a)*Ym(Beta_mn r)
   double R = 0;
   R = ym_Betamn_a[m][n] * jn(m, Betamn[m][n] * r) - jm_Betamn_a[m][n] * yn(m, Betamn[m][n] * r);
 
   if (verbosity > 100)
   {
     std::cout << R << std::endl;
   }
   return R;
 }
 
 double Rossegger::Rmn_(int m, int n, double r)
 {
   if (verbosity > 100)
   {
     std::cout << "Determine Rmn(" << m << "," << n << "," << r << ") = ";
   }
 
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rmn(" << m << "," << n << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   //  Calculate the function using C-libraries from boost
   //  Rossegger Equation 5.11:
   //         Rmn(r) = Ym(Beta_mn a)*Jm(Beta_mn r) - Jm(Beta_mn a)*Ym(Beta_mn r)
   double R = 0;
   R = yn(m, Betamn[m][n] * a) * jn(m, Betamn[m][n] * r) - jn(m, Betamn[m][n] * a) * yn(m, Betamn[m][n] * r);
 
   if (verbosity > 100)
   {
     std::cout << R << std::endl;
   }
   return R;
 }
 
 double Rossegger::Rmn1(int m, int n, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rmn1(" << m << "," << n << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   //  Calculate using the TMath functions from root.
   //  Rossegger Equation 5.32
   //         Rmn1(r) = Km(BetaN a)Im(BetaN r) - Im(BetaN a) Km(BetaN r)
   double R = 0;
   R = km_BetaN_a[m][n] * in(m, BetaN[n] * r) - im_BetaN_a[m][n] * kn(m, BetaN[n] * r);
 
   return R;
 }
 
 double Rossegger::Rmn1_(int m, int n, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rmn1(" << m << "," << n << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   //  Calculate using the TMath functions from root.
   //  Rossegger Equation 5.32
   //         Rmn1(r) = Km(BetaN a)Im(BetaN r) - Im(BetaN a) Km(BetaN r)
   double R = 0;
   double BetaN_ = (n + 1) * pi / L;
   R = kn(m, BetaN_ * a) * in(m, BetaN_ * r) - in(m, BetaN_ * a) * kn(m, BetaN_ * r);
 
   return R;
 }
 
 double Rossegger::Rmn2(int m, int n, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rmn2(" << m << "," << n << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   //  Calculate using the TMath functions from root.
   //  Rossegger Equation 5.33
   //         Rmn2(r) = Km(BetaN b)Im(BetaN r) - Im(BetaN b) Km(BetaN r)
   double R = 0;
   R = km_BetaN_b[m][n] * in(m, BetaN[n] * r) - im_BetaN_b[m][n] * kn(m, BetaN[n] * r);
 
   return R;
 }
 
 double Rossegger::Rmn2_(int m, int n, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rmn2(" << m << "," << n << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   //  Calculate using the TMath functions from root.
   //  Rossegger Equation 5.33
   //         Rmn2(r) = Km(BetaN b)Im(BetaN r) - Im(BetaN b) Km(BetaN r)
   double R = 0;
   double BetaN_ = (n + 1) * pi / L;
   R = kn(m, BetaN_ * b) * in(m, BetaN_ * r) - in(m, BetaN_ * b) * kn(m, BetaN_ * r);
 
   return R;
 }
 
 double Rossegger::RPrime(int m, int n, double ref, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (ref < a || ref > b)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments RPrime(" << m << "," << n << "," << ref << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   double R = 0;
   //  Calculate using the TMath functions from root.
   //  Rossegger Equation 5.65
   //         Rmn2(ref,r) = BetaN/2* [   Km(BetaN ref) {Im-1(BetaN r) + Im+1(BetaN r)}
   //                                  - Im(BetaN ref) {Km-1(BetaN r) + Km+1(BetaN r)}  ]
   //  NOTE:  K-m(z) = Km(z) and I-m(z) = Im(z)... though boost handles negative orders.
   //
   // with: s -> ref,  t -> r,
   double BetaN_ = BetaN[n];
   double term1 = kn(m, BetaN_ * ref) * (in(m - 1, BetaN_ * r) + in(m + 1, BetaN_ * r));
   double term2 = in(m, BetaN_ * ref) * (kn(m - 1, BetaN_ * r) + kn(m + 1, BetaN_ * r));
   R = BetaN_ / 2.0 * (term1 + term2);
 
   return R;
 }
 
 double Rossegger::RPrime_(int m, int n, double ref, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (m < 0 || m >= NumberOfOrders)
   {
     error = 1;
   }
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (ref < a || ref > b)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments RPrime(" << m << "," << n << "," << ref << "," << r << ")" << std::endl;
     ;
     return 0;
   }
 
   double R = 0;
   //  Rossegger Equation 5.65
   //         Rmn2(ref,r) = BetaN/2* [   Km(BetaN ref) {Im-1(BetaN r) + Im+1(BetaN r)}
   //                                  - Im(BetaN ref) {Km-1(BetaN r) + Km+1(BetaN r)}  ]
   //  NOTE:  K-m(z) = Km(z) and I-m(z) = Im(z)... though boost handles negative orders.
   //
   // with: s -> ref,  t -> r,
   double BetaN_ = (n + 1) * pi / L;
   double term1 = kn(m, BetaN_ * ref) * (in(m - 1, BetaN_ * r) + in(m + 1, BetaN_ * r));
   double term2 = in(m, BetaN_ * ref) * (kn(m - 1, BetaN_ * r) + kn(m + 1, BetaN_ * r));
   R = BetaN_ / 2.0 * (term1 + term2);
 
   return R;
 }
 
 double Rossegger::Rnk_for_zeroes(int n, double mu)
 {
   // unlike Rossegger, we count 'k' and 'n' from zero.
   if (verbosity > 10)
   {
     std::cout << "Rnk_for_zeroes called with n=" << n << ",mu=" << mu << std::endl;
   }
   double betana = BetaN_a[n];
   double betanb = BetaN_b[n];
   //  Rossegger Equation 5.46
   //       Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN b) - Kimu_nk(BetaN a) Limu_nk (BetaN b)
 
   return limu(mu, betana) * kimu(mu, betanb) - kimu(mu, betana) * limu(mu, betanb);
 }
 
 double Rossegger::Rnk_for_zeroes_(int n, double mu)
 {
   // unlike Rossegger, we count 'k' and 'n' from zero.
   if (verbosity > 10)
   {
     std::cout << "Rnk_for_zeroes called with n=" << n << ",mu=" << mu << std::endl;
   }
   double BetaN_ = (n + 1) * pi / L;  // this is defined in the paragraph before 5.46
   double betana = BetaN_ * a;
   double betanb = BetaN_ * b;
   //  Rossegger Equation 5.46
   //       Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN b) - Kimu_nk(BetaN a) Limu_nk (BetaN b)
 
   return limu(mu, betana) * kimu(mu, betanb) - kimu(mu, betana) * limu(mu, betanb);
 }
 
 double Rossegger::Rnk(int n, int k, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (k < 0 || k >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rnk(" << n << "," << k << "," << r << ")" << std::endl;
     ;
     return 0;
   }
   //  Rossegger Equation 5.45
   //       Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN r) - Kimu_nk(BetaN a) Limu_nk (BetaN r)
 
   return liMunk_BetaN_a[n][k] * kimu(Munk[n][k], BetaN[n] * r) - kiMunk_BetaN_a[n][k] * limu(Munk[n][k], BetaN[n] * r);
 }
 
 double Rossegger::Rnk_(int n, int k, double r)
 {
   //  Check input arguments for sanity...
   int error = 0;
   if (n < 0 || n >= NumberOfOrders)
   {
     error = 1;
   }
   if (k < 0 || k >= NumberOfOrders)
   {
     error = 1;
   }
   if (r < a || r > b)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Rnk(" << n << "," << k << "," << r << ")" << std::endl;
     ;
     return 0;
   }
   double BetaN_ = (n + 1) * pi / L;
   //  Rossegger Equation 5.45
   //       Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN r) - Kimu_nk(BetaN a) Limu_nk (BetaN r)
 
   return limu(Munk[n][k], BetaN_ * a) * kimu(Munk[n][k], BetaN_ * r) - kimu(Munk[n][k], BetaN_ * a) * limu(Munk[n][k], BetaN_ * r);
 }
 
 double Rossegger::Ez(double r, double phi, double z, double r1, double phi1, double z1)
 {
   // rcc streamlined Ez
   int error = 0;
   if (r < a || r > b)
   {
     error = 1;
   }
   if (phi < 0 || phi > 2 * pi)
   {
     error = 1;
   }
   if (z < 0 || z > L)
   {
     error = 1;
   }
   if (r1 < a || r1 > b)
   {
     error = 1;
   }
   if (phi1 < 0 || phi1 > 2 * pi)
   {
     error = 1;
   }
   if (z1 < 0 || z1 > L)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Ez(";
     std::cout << r << ",";
     std::cout << phi << ",";
     std::cout << z << ",";
     std::cout << r1 << ",";
     std::cout << phi1 << ",";
     std::cout << z1;
     std::cout << ")" << std::endl;
     ;
     return 0;
   }
   // Rossegger Equation 5.64
   double G = 0;
   for (int m = 0; m < NumberOfOrders; m++)
   {
     if (verbosity > 10)
     {
       std::cout << std::endl
                 << m;
     }
     for (int n = 0; n < NumberOfOrders; n++)
     {
       if (verbosity > 10)
       {
         std::cout << " " << n;
       }
       double term = 1;  // unitless
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       term *= (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1));  // unitless
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       term *= Rmn(m, n, r) * Rmn(m, n, r1) / N2mn[m][n];  // units of 1/[L]^2
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       if (z < z1)
       {
         term *= cosh(Betamn[m][n] * z) * sinh(Betamn[m][n] * (L - z1)) / sinh_Betamn_L[m][n];  // unitless
       }
       else
       {
         term *= -cosh(Betamn[m][n] * (L - z)) * sinh(Betamn[m][n] * z1) / sinh_Betamn_L[m][n];  // unitless
       }
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       G += term;
       if (verbosity > 10)
       {
         std::cout << " " << term << " " << G << std::endl;
       }
     }
   }
 
   G = G / (2.0 * pi);
   if (verbosity)
   {
     std::cout << "Ez = " << G << std::endl;
   }
 
   return G;
 }
 
 double Rossegger::Ez_(double r, double phi, double z, double r1, double phi1, double z1)
 {
   // if(fByFile && fabs(r-r1)>MinimumDR && fabs(z-z1)>MinimumDZ) return ByFileEZ(r,phi,z,r1,phi1,z1);
   //   Check input arguments for sanity...
   int error = 0;
   if (r < a || r > b)
   {
     error = 1;
   }
   if (phi < 0 || phi > 2 * pi)
   {
     error = 1;
   }
   if (z < 0 || z > L)
   {
     error = 1;
   }
   if (r1 < a || r1 > b)
   {
     error = 1;
   }
   if (phi1 < 0 || phi1 > 2 * pi)
   {
     error = 1;
   }
   if (z1 < 0 || z1 > L)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Ez(";
     std::cout << r << ",";
     std::cout << phi << ",";
     std::cout << z << ",";
     std::cout << r1 << ",";
     std::cout << phi1 << ",";
     std::cout << z1;
     std::cout << ")" << std::endl;
     ;
     return 0;
   }
 
   double G = 0;
   for (int m = 0; m < NumberOfOrders; m++)
   {
     if (verbosity > 10)
     {
       std::cout << std::endl
                 << m;
     }
     for (int n = 0; n < NumberOfOrders; n++)
     {
       if (verbosity > 10)
       {
         std::cout << " " << n;
       }
       double term = 1 / (2.0 * pi);
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       term *= (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1));  // unitless
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       term *= Rmn(m, n, r) * Rmn(m, n, r1) / N2mn[m][n];  // units of 1/[L]^2
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       if (z < z1)
       {
         term *= cosh(Betamn[m][n] * z) * sinh(Betamn[m][n] * (L - z1)) / sinh(Betamn[m][n] * L);
       }
       else
       {
         term *= -cosh(Betamn[m][n] * (L - z)) * sinh(Betamn[m][n] * z1) / sinh(Betamn[m][n] * L);
         ;
       }
       if (verbosity > 10)
       {
         std::cout << " " << term;
       }
       G += term;
       if (verbosity > 10)
       {
         std::cout << " " << term << " " << G << std::endl;
       }
     }
   }
   if (verbosity)
   {
     std::cout << "Ez = " << G << std::endl;
   }
 
   return G;
 }
 
 double Rossegger::Er(double r, double phi, double z, double r1, double phi1, double z1)
 {
   // rcc streamlined Er
 
   // as in Rossegger 5.65
   // field at r, phi, z due to unit charge at r1, phi1, z1;
   // if(fByFile && fabs(r-r1)>MinimumDR && fabs(z-z1)>MinimumDZ) return ByFileER(r,phi,z,r1,phi1,z1);
   //   Check input arguments for sanity...
   int error = 0;
   if (r < a || r > b)
   {
     error = 1;
   }
   if (phi < 0 || phi > 2 * pi)
   {
     error = 1;
   }
   if (z < 0 || z > L)
   {
     error = 1;
   }
   if (r1 < a || r1 > b)
   {
     error = 1;
   }
   if (phi1 < 0 || phi1 > 2 * pi)
   {
     error = 1;
   }
   if (z1 < 0 || z1 > L)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Er(";
     std::cout << r << ",";
     std::cout << phi << ",";
     std::cout << z << ",";
     std::cout << r1 << ",";
     std::cout << phi1 << ",";
     std::cout << z1;
     std::cout << ")" << std::endl;
     ;
     return 0;
   }
 
   double part = 0;
   double G = 0;
   for (int m = 0; m < NumberOfOrders; m++)
   {
     for (int n = 0; n < NumberOfOrders; n++)
     {
       double term = 1;
       part = (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1));  // unitless
       if (verbosity > 10)
       {
         std::cout << "(2 - ((m==0)?1:0))*cos(m*(phi-phi1)); = " << part << std::endl;
       }
       term *= part;
       part = sin(BetaN[n] * z) * sin(BetaN[n] * z1);  // unitless
       if (verbosity > 10)
       {
         std::cout << "sin(BetaN[n]*z)*sin(BetaN[n]*z1); = " << part << std::endl;
       }
       term *= part;
 
       if (r < r1)
       {
         term *= RPrime(m, n, a, r) * Rmn2(m, n, r1);  // units of 1/[L]
       }
       else
       {
         term *= Rmn1(m, n, r1) * RPrime(m, n, b, r);  // units of 1/[L]
       }
       term /= bessel_denominator[m][n];  // unitless
       G += term;
     }
   }
 
   G = G / (L * pi);  // units of 1/[L] -- net is 1/[L]^2
   if (verbosity)
   {
     std::cout << "Er = " << G << std::endl;
   }
 
   return G;
 }
 
 double Rossegger::Er_(double r, double phi, double z, double r1, double phi1, double z1)
 {
   // doesn't take advantage of precalcs.
   // as in Rossegger 5.65
   // field at r, phi, z due to unit charge at r1, phi1, z1;
   // if(fByFile && fabs(r-r1)>MinimumDR && fabs(z-z1)>MinimumDZ) return ByFileER(r,phi,z,r1,phi1,z1);
   //   Check input arguments for sanity...
   int error = 0;
   if (r < a || r > b)
   {
     error = 1;
   }
   if (phi < 0 || phi > 2 * pi)
   {
     error = 1;
   }
   if (z < 0 || z > L)
   {
     error = 1;
   }
   if (r1 < a || r1 > b)
   {
     error = 1;
   }
   if (phi1 < 0 || phi1 > 2 * pi)
   {
     error = 1;
   }
   if (z1 < 0 || z1 > L)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Er(";
     std::cout << r << ",";
     std::cout << phi << ",";
     std::cout << z << ",";
     std::cout << r1 << ",";
     std::cout << phi1 << ",";
     std::cout << z1;
     std::cout << ")" << std::endl;
     ;
     return 0;
   }
 
   double G = 0;
   for (int m = 0; m < NumberOfOrders; m++)
   {
     for (int n = 0; n < NumberOfOrders; n++)
     {
       double term = 1 / (L * pi);
       term *= (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1));
       double BetaN_ = (n + 1) * pi / L;
       term *= sin(BetaN_ * z) * sin(BetaN_ * z1);
       if (r < r1)
       {
         term *= RPrime_(m, n, a, r) * Rmn2_(m, n, r1);
       }
       else
       {
         term *= Rmn1_(m, n, r1) * RPrime_(m, n, b, r);
       }
 
       term /= TMath::BesselI(m, BetaN_ * a) * TMath::BesselK(m, BetaN_ * b) - TMath::BesselI(m, BetaN_ * b) * TMath::BesselK(m, BetaN_ * a);
 
       G += term;
     }
   }
 
   if (verbosity)
   {
     std::cout << "Er = " << G << std::endl;
   }
 
   return G;
 }
 
 double Rossegger::Ephi(double r, double phi, double z, double r1, double phi1, double z1)
 {
   // rcc streamlined Ephi term
   // compute field at rphiz from charge at r1phi1z1
   //   Check input arguments for sanity...
   int error = 0;
   if (r < a || r > b)
   {
     error = 1;
   }
   if (phi < 0 || phi > 2 * pi)
   {
     error = 1;
   }
   if (z < 0 || z > L)
   {
     error = 1;
   }
   if (r1 < a || r1 > b)
   {
     error = 1;
   }
   if (phi1 < 0 || phi1 > 2 * pi)
   {
     error = 1;
   }
   if (z1 < 0 || z1 > L)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Ephi(";
     std::cout << r << ",";
     std::cout << phi << ",";
     std::cout << z << ",";
     std::cout << r1 << ",";
     std::cout << phi1 << ",";
     std::cout << z1;
     std::cout << ")" << std::endl;
     ;
     return 0;
   }
 
   double G = 0;
   // Rossegger Eqn. 5.66:
   for (int k = 0; k < NumberOfOrders; k++)
   {
     for (int n = 0; n < NumberOfOrders; n++)
     {
       double term = 1;
       term *= sin(BetaN[n] * z) * sin(BetaN[n] * z1);     // unitless
       term *= Rnk(n, k, r) * Rnk(n, k, r1) / N2nk[n][k];  // unitless?
 
       // the derivative of cosh(munk(pi-|phi-phi1|)
       if (phi > phi1)
       {
         term *= -sinh(Munk[n][k] * (pi - (phi - phi1)));  // unitless
       }
       else
       {
         term *= sinh(Munk[n][k] * (pi - (phi1 - phi)));  // unitless
       }
       term *= 1 / sinh_pi_Munk[n][k];  // unitless
       G += term;
     }
   }
 
   G = G / (L * r);  // units of 1/[L]^2.  r comes from the phi term in cylindrical gradient expression.
   if (verbosity)
   {
     std::cout << "Ephi = " << G << std::endl;
   }
 
   return G;
 }
 
 double Rossegger::Ephi_(double r, double phi, double z, double r1, double phi1, double z1)
 {
   // compute field at rphiz from charge at r1phi1z1
   //   Check input arguments for sanity...
   int error = 0;
   if (r < a || r > b)
   {
     error = 1;
   }
   if (phi < 0 || phi > 2 * pi)
   {
     error = 1;
   }
   if (z < 0 || z > L)
   {
     error = 1;
   }
   if (r1 < a || r1 > b)
   {
     error = 1;
   }
   if (phi1 < 0 || phi1 > 2 * pi)
   {
     error = 1;
   }
   if (z1 < 0 || z1 > L)
   {
     error = 1;
   }
   if (error)
   {
     std::cout << "Invalid arguments Ephi(";
     std::cout << r << ",";
     std::cout << phi << ",";
     std::cout << z << ",";
     std::cout << r1 << ",";
     std::cout << phi1 << ",";
     std::cout << z1;
     std::cout << ")" << std::endl;
     ;
     return 0;
   }
 
   verbosity = 1;
   double G = 0;
   // Rossegger Eqn. 5.66:
   for (int k = 0; k < NumberOfOrders; k++)  // off by one from Rossegger convention!
   {
     if (verbosity)
     {
       std::cout << "\nk=" << k;
     }
     for (int n = 0; n < NumberOfOrders; n++)  // off by one from Rossegger convention!
     {
       if (verbosity)
       {
         std::cout << " n=" << n;
       }
       double BetaN_ = (n + 1) * pi / L;
       double term = 1 / (L * r);
       if (verbosity)
       {
         std::cout << " 1/L=" << term;
       }
       term *= sin(BetaN_ * z) * sin(BetaN_ * z1);
       if (verbosity)
       {
         std::cout << " *sinsin=" << term;
       }
       term *= Rnk_(n, k, r) * Rnk_(n, k, r1);
       if (verbosity)
       {
         std::cout << " *rnkrnk=" << term;
       }
       term /= N2nk[n][k];
       if (verbosity)
       {
         std::cout << " */nnknnk=" << term;
       }
 
       // the derivative of cosh(munk(pi-|phi-phi1|)
       if (phi > phi1)
       {
         term *= -sinh(Munk[n][k] * (pi - (phi - phi1)));
         // term *=  Munk[n][k]*sinh(Munk[n][k]*pi*(phi1-phi));
         // this originally has a factor of Munk in front, but that cancels with one in the denominator
       }
       else
       {
         term *= sinh(Munk[n][k] * (pi - (phi1 - phi)));
         // term *=  -Munk[n][k]*sinh(Munk[n][k]*pi*(phi-phi1));
         // this originally has a factor of Munk in front, but that cancels with one in the denominator
       }
       if (verbosity)
       {
         std::cout << " *sinh(mu*pi-phi-phi)=" << term;
       }
       term *= 1 / (sinh(pi * Munk[n][k]));
       // term *= 1/(Munk[n][k]*sinh(pi*Munk[n][k]));
       // this originally has a factor of Munk in front, but that cancels with one in the numerator
       G += term;
       if (verbosity)
       {
         std::cout << "  /sinh=" << term << " G=" << G << std::endl;
       }
     }
   }
   if (verbosity)
   {
     std::cout << "Ephi = " << G << std::endl;
   }
   verbosity = 0;
 
   return G;
 }
 
 void Rossegger::SaveZeroes(const std::string &destfile)
 {
   TFile *output = TFile::Open(destfile.c_str(), "RECREATE");
   output->cd();
 
   TTree *tInfo = new TTree("info", "Mu[n][k] values");
   int ord = NumberOfOrders;
   tInfo->Branch("order", &ord);
   tInfo->Branch("epsilon", &epsilon);
   tInfo->Fill();
 
   int n, k, m;
   double munk, betamn;
   double n2nk, n2mn;
   TTree *tmunk = new TTree("munk", "Mu[n][k] values");
   tmunk->Branch("n", &n);
   tmunk->Branch("k", &k);
   tmunk->Branch("munk", &munk);
   tmunk->Branch("n2nk", &n2nk);
   for (n = 0; n < ord; n++)
   {
     for (k = 0; k < ord; k++)
     {
       munk = Munk[n][k];
       n2nk = N2nk[n][k];
       tmunk->Fill();
     }
   }
 
   TTree *tbetamn = new TTree("betamn", "Beta[m][n] values");
   tbetamn->Branch("m", &m);
   tbetamn->Branch("n", &n);
   tbetamn->Branch("betamn", &betamn);
   tbetamn->Branch("n2mn", &n2mn);
   for (m = 0; m < ord; m++)
   {
     for (n = 0; n < ord; n++)
     {
       betamn = Betamn[m][n];
       n2mn = N2mn[m][n];
       tbetamn->Fill();
     }
   }
 
   tInfo->Write();
   tmunk->Write();
   tbetamn->Write();
   // output->Write();
   output->Close();
   return;
 }
 
 void Rossegger::LoadZeroes(const std::string &destfile)
 {
   TFile *f = TFile::Open(destfile.c_str(), "READ");
   std::cout << "reading rossegger zeroes from " << destfile << std::endl;
   TTree *tInfo = (TTree *) (f->Get("info"));
   int ord;
   tInfo->SetBranchAddress("order", &ord);
   tInfo->SetBranchAddress("epsilon", &epsilon);
   tInfo->GetEntry(0);
   std::cout << "order=" << ord << ",epsilon=" << epsilon << std::endl;
 
   int n, k, m;
   double munk, betamn;
   double n2nk, n2mn;
   TTree *tmunk = (TTree *) (f->Get("munk"));
   tmunk->SetBranchAddress("n", &n);
   tmunk->SetBranchAddress("k", &k);
   tmunk->SetBranchAddress("munk", &munk);
   tmunk->SetBranchAddress("n2nk", &n2nk);
   for (int i = 0; i < tmunk->GetEntries(); i++)
   {
     tmunk->GetEntry(i);
     Munk[n][k] = munk;
     N2nk[n][k] = n2nk;
   }
 
   TTree *tbetamn = (TTree *) (f->Get("betamn"));
   tbetamn->SetBranchAddress("m", &m);
   tbetamn->SetBranchAddress("n", &n);
   tbetamn->SetBranchAddress("betamn", &betamn);
   tbetamn->SetBranchAddress("n2mn", &n2mn);
   for (int i = 0; i < tbetamn->GetEntries(); i++)
   {
     tbetamn->GetEntry(i);
     Betamn[m][n] = betamn;
     N2mn[m][n] = n2mn;
   }
 
   f->Close();
   return;
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team