ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Rossegger.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Rossegger.cc
1 #include "Rossegger.h"
2 
3 #include <TFile.h>
4 #include <TMath.h>
5 #include <TTree.h>
6 
7 #include <boost/math/special_functions.hpp> //covers all the special functions.
8 
9 #include <algorithm> // for max
10 #include <cassert> // for assert
11 #include <cmath>
12 #include <cstdlib> // for exit, abs
13 #include <fstream>
14 #include <iostream>
15 #include <sstream>
16 #include <string>
17 
18 // the limu and kimu terms, that i need to think about a little while longer...
19 extern "C"
20 {
21  void dkia_(int *IFAC, double *X, double *A, double *DKI, double *DKID, int *IERRO);
22  void dlia_(int *IFAC, double *X, double *A, double *DLI, double *DLID, int *IERRO);
23 }
24 //
25 
26 //Bessel Function J_n(x):
27 #define jn(order, x) boost::math::cyl_bessel_j(order, x)
28 //Bessel (Neumann) Function Y_n(x):
29 #define yn(order, x) boost::math::cyl_neumann(order, x)
30 //Modified Bessel Function of the first kind I_n(x):
31 #define in(order, x) boost::math::cyl_bessel_i(order, x)
32 //Modified Bessel Function of the second kind K_n(x):
33 #define kn(order, x) boost::math::cyl_bessel_k(order, x)
34 #define limu(im_order, x) Rossegger::Limu(im_order, x)
35 #define kimu(im_order, x) Rossegger::Kimu(im_order, x)
36 
37 /*
38  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.
39  */
40 
41 Rossegger::Rossegger(double InnerRadius, double OuterRadius, double Rdo_Z, double precision)
42 {
43  a = InnerRadius;
44  b = OuterRadius;
45  L = Rdo_Z;
46 
48 
49  verbosity = 0;
50  pi = M_PI;
51 
53 
54  //load the greens functions:
55  char zeroesfilename[200];
56  sprintf(zeroesfilename, "rosseger_zeroes_eps%1.0E_a%2.2f_b%2.2f_L%2.2f.root", epsilon, a, b, L);
57  TFile *fileptr = TFile::Open(zeroesfilename, "READ");
58  if (!fileptr)
59  { //generate the lookuptable
62  SaveZeroes(zeroesfilename);
63  }
64  else
65  { //load it from a file
66  fileptr->Close();
67  LoadZeroes(zeroesfilename);
68  }
69 
71 
72  std::cout << "Rossegger object initialized as follows:" << std::endl;
73  std::cout << " Inner Radius = " << a << " cm." << std::endl;
74  std::cout << " Outer Radius = " << b << " cm." << std::endl;
75  std::cout << " Half Length = " << L << " cm." << std::endl;
76 
77  if (!CheckZeroes(0.01))
78  {
79  std::cout << "CheckZeroes(0.01) returned false, exiting" << std::endl;
80  exit(1);
81  }
82  return;
83 }
84 
85 double Rossegger::FindNextZero(double xstart, double epsilon, int order, double (Rossegger::*func)(int, double))
86 {
87  double previous = (this->*func)(order, xstart);
88  double x = xstart + epsilon;
89  double value = previous;
90 
91  while (!((value == 0) || (value < 0 && previous > 0) || (value > 0 && previous < 0)))
92  {
93  // Rossegger equation 5.12
94  value = (this->*func)(order, x);
95  if (value == 0) std::cout << "hit it exactly! Go buy a lottery ticket!" << std::endl;
96  if ((value == 0) || (value < 0 && previous > 0) || (value > 0 && previous < 0))
97  {
98  //when we go from one sign to the other, we have bracketed the zero
99  //the following is mathematically equivalent to finding the delta
100  //between the x of our previous value and the point where we cross the axis
101  //then returning x0=x_old+delta
102  double slope = (value - previous) / epsilon;
103  double intercept = value - slope * x;
104  double x0 = -intercept / slope;
105  if (verbosity > 1) std::cout << " " << x0 << "," << std::endl;
106  double n0 = (this->*func)(order, x - epsilon);
107  double n1 = (this->*func)(order, x + epsilon);
108  if ((n0 < 0 && n1 < 0) || (n0 > 0 && n1 > 0))
109  {
110  printf("neighbors on both sides have the same sign. Check your function resolution!\n");
111  }
112 
113  return x0;
114  }
115  previous = value;
116  x += epsilon;
117  }
118  std::cout << "logic break!\n";
119  assert(1 == 2);
120  return 0;
121 }
122 
123 void Rossegger::FindBetamn(double epsilon)
124 {
125  std::cout << "Now filling the Beta[m][n] Array..." << std::endl;
126  if (verbosity > 5) std::cout << "numberOfOrders= " << NumberOfOrders << std::endl;
127  for (int m = 0; m < NumberOfOrders; m++)
128  {
129  if (verbosity) std::cout << "Filling Beta[" << m << "][n]..." << std::endl;
130 
131  double x = epsilon;
132  for (int n = 0; n < NumberOfOrders; n++)
133  { // !!! Off by one from Rossegger convention !!!
134  x = FindNextZero(x, epsilon, m, &Rossegger::Rmn_for_zeroes);
135  Betamn[m][n] = x / b;
136  x += epsilon;
137  }
138  }
139 
140  // Now fill in the N2mn array...
141  for (int m = 0; m < NumberOfOrders; m++)
142  {
143  for (int n = 0; n < NumberOfOrders; n++)
144  {
145  // Rossegger Equation 5.17
146  // N^2_mn = 2/(pi * beta)^2 [ {Jm(beta a)/Jm(beta b)}^2 - 1 ]
147  N2mn[m][n] = 2 / (pi * pi * Betamn[m][n] * Betamn[m][n]);
148  //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);
149  double jna_over_jnb = jn(m, Betamn[m][n] * a) / jn(m, Betamn[m][n] * b);
150  N2mn[m][n] *= (jna_over_jnb * jna_over_jnb - 1.0);
151  //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
152  if (verbosity > 1) std::cout << "m: " << m << " n: " << n << " N2[m][n]: " << N2mn[m][n];
153  double integral = 0.0;
154  double step = 0.01;
155  if (verbosity > 1)
156  {
157  for (double r = a; r < b; r += step)
158  {
159  double rmnval = Rmn(m, n, r);
160 
161  integral += rmnval * rmnval * r * step; //Rmn(m,n,r)*Rmn(m,n,r)*r*step;
162  }
163  std::cout << " Int: " << integral << std::endl;
164  }
165  //N2mn[m][n] = integral;
166  }
167  }
168 
169  std::cout << "Done." << std::endl;
170 }
171 
172 void Rossegger::FindMunk(double epsilon)
173 {
174  std::cout << "Now filling the Mu[n][k] Array..." << std::endl;
175  if (verbosity > 5) std::cout << "numberOfOrders= " << NumberOfOrders << std::endl;
176  // We're looking for the zeroes of Rossegger eqn. 5.46:
177  // R_nk(mu_nk;a,b)=Limu(Beta_n*a)Kimu(Beta_n*b)-Kimu(Beta_n*a)Limu(Beta_n*b)=0
178  // since a and b are fixed, R_nk is a function solely of mu_nk and n.
179  // for each 'n' we wish to find the a set of k mu_n's that result in R_nk=0
180 
181  //could add an option here to load the munks from file if the dimensions match.
182 
183  for (int n = 0; n < NumberOfOrders; n++) // !!! Off by one from Rossegger convention !!!
184  {
185  if (verbosity) std::cout << "Filling Mu[" << n << "][k]..." << std::endl;
186  double x = epsilon;
187  for (int k = 0; k < NumberOfOrders; k++)
188  {
189  x = FindNextZero(x, epsilon, n, &Rossegger::Rnk_for_zeroes);
190  Munk[n][k] = x;
191  x += epsilon;
192  if (verbosity > 0)
193  {
194  printf("Mu[%d][%d]=%E\n", n, k, Munk[n][k]);
195  printf("adjacent values are Rnk[mu-epsilon]=%E\tRnk[mu+epsilon]=%E\n",
196  Rnk_for_zeroes(n, x - epsilon), Rnk_for_zeroes(n, x + epsilon));
197  if (verbosity > 100) printf("values of argument to limu and kimu are %f and %f\n",
198  (n + 1) * pi / L * a, (n + 1) * pi / L * b);
199  }
200  }
201  }
202 
203  // Now fill in the N2nk array...
204  for (int n = 0; n < NumberOfOrders; n++)
205  {
206  for (int k = 0; k < NumberOfOrders; k++)
207  {
208  // Rossegger Equation 5.48
209  // Integral of R_nk(r)*R_ns(r) dr/r= delta_ks*N2nk
210  // note that unlike N2mn, there is no convenient shortcut here.
211  double integral = 0.0;
212  double step = 0.001;
213 
214  for (double r = a; r < b; r += step)
215  {
216  double rnkval = Rnk_(n, k, r); //must used un-optimized. We don't have those values yet...
217 
218  integral += rnkval * rnkval / r * step;
219  }
220  if (verbosity > 1)
221  {
222  std::cout << " Int: " << integral << std::endl;
223  }
224  N2nk[n][k] = integral;
225  if (verbosity > 1) std::cout << "n: " << n << " k: " << k << " N2nk[n][k]: " << N2nk[n][k];
226  }
227  }
228 
229  std::cout << "Done." << std::endl;
230  return;
231 }
232 
233 bool Rossegger::CheckZeroes(double epsilon)
234 {
235  //confirm that the tabulated zeroes are all zeroes of their respective functions:
236  double result;
237  for (int m = 0; m < NumberOfOrders; m++)
238  {
239  for (int n = 0; n < NumberOfOrders; n++)
240  { // !!! Off by one from Rossegger convention !!!
241  result = Rmn_for_zeroes(m, Betamn[m][n] * b);
242  if (abs(result) > epsilon)
243  {
244  printf("(m=%d,n=%d) Jm(x)Ym(lx)-Jm(lx)Ym(x) = %f for x=b*%f\n", m, n, result, Betamn[m][n]);
245  return false;
246  }
247  }
248  }
249 
250  // R_nk(mu_nk;a,b)=Limu(Beta_n*a)Kimu(Beta_n*b)-Kimu(Beta_n*a)Limu(Beta_n*b)=0
251  for (int n = 0; n < NumberOfOrders; n++)
252  {
253  for (int k = 0; k < NumberOfOrders; k++)
254  { // !!! Off by one from Rossegger convention !!!
255  result = Rnk_for_zeroes(n, Munk[n][k]);
256  if (abs(result) > epsilon * 100)
257  {
258  printf("(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",
259  n, k, result, Munk[n][k]);
260  return false;
261  }
262  }
263  }
264 
265  return true;
266 }
267 
269 { //Routine used to fill the arrays of other values used repeatedly
270  //these constants depend only on the geometry of the detector
271  printf("Precalcing %d geometric constants\n", 3 * NumberOfOrders + 5 * NumberOfOrders * NumberOfOrders);
272  for (int n = 0; n < NumberOfOrders; n++)
273  {
274  BetaN[n] = (n + 1) * pi / L; // BetaN=(n+1)*pi/L as used in eg 5.32, .46
275  BetaN_a[n] = BetaN[n] * a; // BetaN*a as used in eg 5.32, .46
276  BetaN_b[n] = BetaN[n] * b; // BetaN*b as used in eg 5.32, .46
277  for (int m = 0; m < NumberOfOrders; m++)
278  {
279  km_BetaN_a[m][n] = kn(m, BetaN_a[n]); //kn(m,BetaN*a) as used in Rossegger 5.32
280  im_BetaN_a[m][n] = in(m, BetaN_a[n]); //in(m,BetaN*a) as used in Rossegger 5.32
281  km_BetaN_b[m][n] = kn(m, BetaN_b[n]); //kn(m,BetaN*b) as used in Rossegger 5.33
282  im_BetaN_b[m][n] = in(m, BetaN_b[n]); //in(m,BetaN*b) as used in Rossegger 5.33
283  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
284  }
285  }
286  return;
287 }
289 { //Routine used to fill the arrays of other values used repeatedly
290  //these constants depend on geometry and the zeroes of special functions
291  printf("Precalcing %d geometric constants\n", 6 * NumberOfOrders * NumberOfOrders);
292 
293  for (int n = 0; n < NumberOfOrders; n++)
294  {
295  for (int m = 0; m < NumberOfOrders; m++)
296  {
297  ym_Betamn_a[m][n] = yn(m, Betamn[m][n] * a); //yn(m,Betamn[m][n]*a) as used in Rossegger 5.11
298  jm_Betamn_a[m][n] = jn(m, Betamn[m][n] * a); //jn(m,Betamn[m][n]*a) as used in Rossegger 5.11
299  sinh_Betamn_L[m][n] = sinh(Betamn[m][n] * L); //sinh(Betamn[m][n]*L) as in Rossegger 5.64
300  }
301  for (int k = 0; k < NumberOfOrders; k++)
302  {
303  liMunk_BetaN_a[n][k] = limu(Munk[n][k], BetaN_a[n]); //limu(Munk[n][k],BetaN*a) as used in Rossegger 5.45
304  kiMunk_BetaN_a[n][k] = kimu(Munk[n][k], BetaN_a[n]); //kimu(Munk[n][k],BetaN*a) as used in Rossegger 5.45
305  sinh_pi_Munk[n][k] = sinh(pi * Munk[n][k]); //sinh(pi*Munk[n][k]) as in Rossegger 5.66
306  }
307  }
308  return;
309 }
310 
311 double Rossegger::Limu(double mu, double x)
312 {
313  //defined in Rossegger eqn 5.44, also a canonical 'satisfactory companion' to Kimu.
314  //could use Griddit?
315  // Rossegger Equation 5.45
316  // Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN r) - Kimu_nk(BetaN a) Limu_nk (BetaN r)
317 
318  int IFAC = 1;
319  double A = mu;
320  double DLI = 0;
321  double DERR = 0;
322  int IERRO = 0;
323 
324  double X = x;
325  dlia_(&IFAC, &X, &A, &DLI, &DERR, &IERRO);
326  return DLI;
327 }
328 double Rossegger::Kimu(double mu, double x)
329 {
330  int IFAC = 1;
331  double A = mu;
332  double DKI = 0;
333  double DERR = 0;
334  int IERRO = 0;
335 
336  double X = x;
337  dkia_(&IFAC, &X, &A, &DKI, &DERR, &IERRO);
338  return DKI;
339 }
340 
341 double Rossegger::Rmn_for_zeroes(int m, double x)
342 {
343  double lx = a * x / b;
344  // Rossegger Equation 5.12:
345 
346  return jn(m, x) * yn(m, lx) - jn(m, lx) * yn(m, x);
347 }
348 
349 double Rossegger::Rmn(int m, int n, double r)
350 {
351  if (verbosity > 100) std::cout << "Determine Rmn(" << m << "," << n << "," << r << ") = ";
352 
353  // Check input arguments for sanity...
354  int error = 0;
355  if (m < 0 || m >= NumberOfOrders) error = 1;
356  if (n < 0 || n >= NumberOfOrders) error = 1;
357  if (r < a || r > b) error = 1;
358  if (error)
359  {
360  std::cout << "Invalid arguments Rmn(" << m << "," << n << "," << r << ")" << std::endl;
361  ;
362  return 0;
363  }
364 
365  // Calculate the function using C-libraries from boost
366  // Rossegger Equation 5.11:
367  // Rmn(r) = Ym(Beta_mn a)*Jm(Beta_mn r) - Jm(Beta_mn a)*Ym(Beta_mn r)
368  double R = 0;
369  R = ym_Betamn_a[m][n] * jn(m, Betamn[m][n] * r) - jm_Betamn_a[m][n] * yn(m, Betamn[m][n] * r);
370 
371  if (verbosity > 100) std::cout << R << std::endl;
372  return R;
373 }
374 
375 double Rossegger::Rmn_(int m, int n, double r)
376 {
377  if (verbosity > 100) std::cout << "Determine Rmn(" << m << "," << n << "," << r << ") = ";
378 
379  // Check input arguments for sanity...
380  int error = 0;
381  if (m < 0 || m >= NumberOfOrders) error = 1;
382  if (n < 0 || n >= NumberOfOrders) error = 1;
383  if (r < a || r > b) error = 1;
384  if (error)
385  {
386  std::cout << "Invalid arguments Rmn(" << m << "," << n << "," << r << ")" << std::endl;
387  ;
388  return 0;
389  }
390 
391  // Calculate the function using C-libraries from boost
392  // Rossegger Equation 5.11:
393  // Rmn(r) = Ym(Beta_mn a)*Jm(Beta_mn r) - Jm(Beta_mn a)*Ym(Beta_mn r)
394  double R = 0;
395  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);
396 
397  if (verbosity > 100) std::cout << R << std::endl;
398  return R;
399 }
400 
401 double Rossegger::Rmn1(int m, int n, double r)
402 {
403  // Check input arguments for sanity...
404  int error = 0;
405  if (m < 0 || m >= NumberOfOrders) error = 1;
406  if (n < 0 || n >= NumberOfOrders) error = 1;
407  if (r < a || r > b) error = 1;
408  if (error)
409  {
410  std::cout << "Invalid arguments Rmn1(" << m << "," << n << "," << r << ")" << std::endl;
411  ;
412  return 0;
413  }
414 
415  // Calculate using the TMath functions from root.
416  // Rossegger Equation 5.32
417  // Rmn1(r) = Km(BetaN a)Im(BetaN r) - Im(BetaN a) Km(BetaN r)
418  double R = 0;
419  R = km_BetaN_a[m][n] * in(m, BetaN[n] * r) - im_BetaN_a[m][n] * kn(m, BetaN[n] * r);
420 
421  return R;
422 }
423 
424 double Rossegger::Rmn1_(int m, int n, double r)
425 {
426  // Check input arguments for sanity...
427  int error = 0;
428  if (m < 0 || m >= NumberOfOrders) error = 1;
429  if (n < 0 || n >= NumberOfOrders) error = 1;
430  if (r < a || r > b) error = 1;
431  if (error)
432  {
433  std::cout << "Invalid arguments Rmn1(" << m << "," << n << "," << r << ")" << std::endl;
434  ;
435  return 0;
436  }
437 
438  // Calculate using the TMath functions from root.
439  // Rossegger Equation 5.32
440  // Rmn1(r) = Km(BetaN a)Im(BetaN r) - Im(BetaN a) Km(BetaN r)
441  double R = 0;
442  double BetaN_ = (n + 1) * pi / L;
443  R = kn(m, BetaN_ * a) * in(m, BetaN_ * r) - in(m, BetaN_ * a) * kn(m, BetaN_ * r);
444 
445  return R;
446 }
447 
448 double Rossegger::Rmn2(int m, int n, double r)
449 {
450  // Check input arguments for sanity...
451  int error = 0;
452  if (m < 0 || m >= NumberOfOrders) error = 1;
453  if (n < 0 || n >= NumberOfOrders) error = 1;
454  if (r < a || r > b) error = 1;
455  if (error)
456  {
457  std::cout << "Invalid arguments Rmn2(" << m << "," << n << "," << r << ")" << std::endl;
458  ;
459  return 0;
460  }
461 
462  // Calculate using the TMath functions from root.
463  // Rossegger Equation 5.33
464  // Rmn2(r) = Km(BetaN b)Im(BetaN r) - Im(BetaN b) Km(BetaN r)
465  double R = 0;
466  R = km_BetaN_b[m][n] * in(m, BetaN[n] * r) - im_BetaN_b[m][n] * kn(m, BetaN[n] * r);
467 
468  return R;
469 }
470 
471 double Rossegger::Rmn2_(int m, int n, double r)
472 {
473  // Check input arguments for sanity...
474  int error = 0;
475  if (m < 0 || m >= NumberOfOrders) error = 1;
476  if (n < 0 || n >= NumberOfOrders) error = 1;
477  if (r < a || r > b) error = 1;
478  if (error)
479  {
480  std::cout << "Invalid arguments Rmn2(" << m << "," << n << "," << r << ")" << std::endl;
481  ;
482  return 0;
483  }
484 
485  // Calculate using the TMath functions from root.
486  // Rossegger Equation 5.33
487  // Rmn2(r) = Km(BetaN b)Im(BetaN r) - Im(BetaN b) Km(BetaN r)
488  double R = 0;
489  double BetaN_ = (n + 1) * pi / L;
490  R = kn(m, BetaN_ * b) * in(m, BetaN_ * r) - in(m, BetaN_ * b) * kn(m, BetaN_ * r);
491 
492  return R;
493 }
494 
495 double Rossegger::RPrime(int m, int n, double ref, double r)
496 {
497  // Check input arguments for sanity...
498  int error = 0;
499  if (m < 0 || m >= NumberOfOrders) error = 1;
500  if (n < 0 || n >= NumberOfOrders) error = 1;
501  if (ref < a || ref > b) error = 1;
502  if (r < a || r > b) error = 1;
503  if (error)
504  {
505  std::cout << "Invalid arguments RPrime(" << m << "," << n << "," << ref << "," << r << ")" << std::endl;
506  ;
507  return 0;
508  }
509 
510  double R = 0;
511  // Calculate using the TMath functions from root.
512  // Rossegger Equation 5.65
513  // Rmn2(ref,r) = BetaN/2* [ Km(BetaN ref) {Im-1(BetaN r) + Im+1(BetaN r)}
514  // - Im(BetaN ref) {Km-1(BetaN r) + Km+1(BetaN r)} ]
515  // NOTE: K-m(z) = Km(z) and I-m(z) = Im(z)... though boost handles negative orders.
516  //
517  // with: s -> ref, t -> r,
518  double BetaN_ = BetaN[n];
519  double term1 = kn(m, BetaN_ * ref) * (in(m - 1, BetaN_ * r) + in(m + 1, BetaN_ * r));
520  double term2 = in(m, BetaN_ * ref) * (kn(m - 1, BetaN_ * r) + kn(m + 1, BetaN_ * r));
521  R = BetaN_ / 2.0 * (term1 + term2);
522 
523  return R;
524 }
525 
526 double Rossegger::RPrime_(int m, int n, double ref, double r)
527 {
528  // Check input arguments for sanity...
529  int error = 0;
530  if (m < 0 || m >= NumberOfOrders) error = 1;
531  if (n < 0 || n >= NumberOfOrders) error = 1;
532  if (ref < a || ref > b) error = 1;
533  if (r < a || r > b) error = 1;
534  if (error)
535  {
536  std::cout << "Invalid arguments RPrime(" << m << "," << n << "," << ref << "," << r << ")" << std::endl;
537  ;
538  return 0;
539  }
540 
541  double R = 0;
542  // Rossegger Equation 5.65
543  // Rmn2(ref,r) = BetaN/2* [ Km(BetaN ref) {Im-1(BetaN r) + Im+1(BetaN r)}
544  // - Im(BetaN ref) {Km-1(BetaN r) + Km+1(BetaN r)} ]
545  // NOTE: K-m(z) = Km(z) and I-m(z) = Im(z)... though boost handles negative orders.
546  //
547  // with: s -> ref, t -> r,
548  double BetaN_ = (n + 1) * pi / L;
549  double term1 = kn(m, BetaN_ * ref) * (in(m - 1, BetaN_ * r) + in(m + 1, BetaN_ * r));
550  double term2 = in(m, BetaN_ * ref) * (kn(m - 1, BetaN_ * r) + kn(m + 1, BetaN_ * r));
551  R = BetaN_ / 2.0 * (term1 + term2);
552 
553  return R;
554 }
555 
556 double Rossegger::Rnk_for_zeroes(int n, double mu)
557 {
558  //unlike Rossegger, we count 'k' and 'n' from zero.
559  if (verbosity > 10) printf("Rnk_for_zeroes called with n=%d,mu=%f\n", n, mu);
560  double betana = BetaN_a[n];
561  double betanb = BetaN_b[n];
562  // Rossegger Equation 5.46
563  // Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN b) - Kimu_nk(BetaN a) Limu_nk (BetaN b)
564 
565  return limu(mu, betana) * kimu(mu, betanb) - kimu(mu, betana) * limu(mu, betanb);
566 }
567 
568 double Rossegger::Rnk_for_zeroes_(int n, double mu)
569 {
570  //unlike Rossegger, we count 'k' and 'n' from zero.
571  if (verbosity > 10) printf("Rnk_for_zeroes called with n=%d,mu=%f\n", n, mu);
572  double BetaN_ = (n + 1) * pi / L; //this is defined in the paragraph before 5.46
573  double betana = BetaN_ * a;
574  double betanb = BetaN_ * b;
575  // Rossegger Equation 5.46
576  // Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN b) - Kimu_nk(BetaN a) Limu_nk (BetaN b)
577 
578  return limu(mu, betana) * kimu(mu, betanb) - kimu(mu, betana) * limu(mu, betanb);
579 }
580 
581 double Rossegger::Rnk(int n, int k, double r)
582 {
583  // Check input arguments for sanity...
584  int error = 0;
585  if (n < 0 || n >= NumberOfOrders) error = 1;
586  if (k < 0 || k >= NumberOfOrders) error = 1;
587  if (r < a || r > b) error = 1;
588  if (error)
589  {
590  std::cout << "Invalid arguments Rnk(" << n << "," << k << "," << r << ")" << std::endl;
591  ;
592  return 0;
593  }
594  // Rossegger Equation 5.45
595  // Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN r) - Kimu_nk(BetaN a) Limu_nk (BetaN r)
596 
597  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);
598 }
599 
600 double Rossegger::Rnk_(int n, int k, double r)
601 {
602  // Check input arguments for sanity...
603  int error = 0;
604  if (n < 0 || n >= NumberOfOrders) error = 1;
605  if (k < 0 || k >= NumberOfOrders) error = 1;
606  if (r < a || r > b) error = 1;
607  if (error)
608  {
609  std::cout << "Invalid arguments Rnk(" << n << "," << k << "," << r << ")" << std::endl;
610  ;
611  return 0;
612  }
613  double BetaN_ = (n + 1) * pi / L;
614  // Rossegger Equation 5.45
615  // Rnk(r) = Limu_nk (BetaN a) Kimu_nk (BetaN r) - Kimu_nk(BetaN a) Limu_nk (BetaN r)
616 
617  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);
618 }
619 
620 double Rossegger::Ez(double r, double phi, double z, double r1, double phi1, double z1)
621 {
622  //rcc streamlined Ez
623  int error = 0;
624  if (r < a || r > b) error = 1;
625  if (phi < 0 || phi > 2 * pi) error = 1;
626  if (z < 0 || z > L) error = 1;
627  if (r1 < a || r1 > b) error = 1;
628  if (phi1 < 0 || phi1 > 2 * pi) error = 1;
629  if (z1 < 0 || z1 > L) error = 1;
630  if (error)
631  {
632  std::cout << "Invalid arguments Ez(";
633  std::cout << r << ",";
634  std::cout << phi << ",";
635  std::cout << z << ",";
636  std::cout << r1 << ",";
637  std::cout << phi1 << ",";
638  std::cout << z1;
639  std::cout << ")" << std::endl;
640  ;
641  return 0;
642  }
643  //Rossegger Equation 5.64
644  double G = 0;
645  for (int m = 0; m < NumberOfOrders; m++)
646  {
647  if (verbosity > 10) std::cout << std::endl
648  << m;
649  for (int n = 0; n < NumberOfOrders; n++)
650  {
651  if (verbosity > 10) std::cout << " " << n;
652  double term = 1; //unitless
653  if (verbosity > 10) std::cout << " " << term;
654  term *= (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1)); //unitless
655  if (verbosity > 10) std::cout << " " << term;
656  term *= Rmn(m, n, r) * Rmn(m, n, r1) / N2mn[m][n]; //units of 1/[L]^2
657  if (verbosity > 10) std::cout << " " << term;
658  if (z < z1)
659  {
660  term *= cosh(Betamn[m][n] * z) * sinh(Betamn[m][n] * (L - z1)) / sinh_Betamn_L[m][n]; //unitless
661  }
662  else
663  {
664  term *= -cosh(Betamn[m][n] * (L - z)) * sinh(Betamn[m][n] * z1) / sinh_Betamn_L[m][n]; //unitless
665  }
666  if (verbosity > 10) std::cout << " " << term;
667  G += term;
668  if (verbosity > 10) std::cout << " " << term << " " << G << std::endl;
669  }
670  }
671 
672  G = G / (2.0 * pi);
673  if (verbosity) std::cout << "Ez = " << G << std::endl;
674 
675  return G;
676 }
677 
678 double Rossegger::Ez_(double r, double phi, double z, double r1, double phi1, double z1)
679 {
680  //if(fByFile && fabs(r-r1)>MinimumDR && fabs(z-z1)>MinimumDZ) return ByFileEZ(r,phi,z,r1,phi1,z1);
681  // Check input arguments for sanity...
682  int error = 0;
683  if (r < a || r > b) error = 1;
684  if (phi < 0 || phi > 2 * pi) error = 1;
685  if (z < 0 || z > L) error = 1;
686  if (r1 < a || r1 > b) error = 1;
687  if (phi1 < 0 || phi1 > 2 * pi) error = 1;
688  if (z1 < 0 || z1 > L) error = 1;
689  if (error)
690  {
691  std::cout << "Invalid arguments Ez(";
692  std::cout << r << ",";
693  std::cout << phi << ",";
694  std::cout << z << ",";
695  std::cout << r1 << ",";
696  std::cout << phi1 << ",";
697  std::cout << z1;
698  std::cout << ")" << std::endl;
699  ;
700  return 0;
701  }
702 
703  double G = 0;
704  for (int m = 0; m < NumberOfOrders; m++)
705  {
706  if (verbosity > 10) std::cout << std::endl
707  << m;
708  for (int n = 0; n < NumberOfOrders; n++)
709  {
710  if (verbosity > 10) std::cout << " " << n;
711  double term = 1 / (2.0 * pi);
712  if (verbosity > 10) std::cout << " " << term;
713  term *= (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1)); //unitless
714  if (verbosity > 10) std::cout << " " << term;
715  term *= Rmn(m, n, r) * Rmn(m, n, r1) / N2mn[m][n]; //units of 1/[L]^2
716  if (verbosity > 10) std::cout << " " << term;
717  if (z < z1)
718  {
719  term *= cosh(Betamn[m][n] * z) * sinh(Betamn[m][n] * (L - z1)) / sinh(Betamn[m][n] * L);
720  }
721  else
722  {
723  term *= -cosh(Betamn[m][n] * (L - z)) * sinh(Betamn[m][n] * z1) / sinh(Betamn[m][n] * L);
724  ;
725  }
726  if (verbosity > 10) std::cout << " " << term;
727  G += term;
728  if (verbosity > 10) std::cout << " " << term << " " << G << std::endl;
729  }
730  }
731  if (verbosity) std::cout << "Ez = " << G << std::endl;
732 
733  return G;
734 }
735 
736 double Rossegger::Er(double r, double phi, double z, double r1, double phi1, double z1)
737 {
738  //rcc streamlined Er
739 
740  //as in Rossegger 5.65
741  //field at r, phi, z due to unit charge at r1, phi1, z1;
742  //if(fByFile && fabs(r-r1)>MinimumDR && fabs(z-z1)>MinimumDZ) return ByFileER(r,phi,z,r1,phi1,z1);
743  // Check input arguments for sanity...
744  int error = 0;
745  if (r < a || r > b) error = 1;
746  if (phi < 0 || phi > 2 * pi) error = 1;
747  if (z < 0 || z > L) error = 1;
748  if (r1 < a || r1 > b) error = 1;
749  if (phi1 < 0 || phi1 > 2 * pi) error = 1;
750  if (z1 < 0 || z1 > L) error = 1;
751  if (error)
752  {
753  std::cout << "Invalid arguments Er(";
754  std::cout << r << ",";
755  std::cout << phi << ",";
756  std::cout << z << ",";
757  std::cout << r1 << ",";
758  std::cout << phi1 << ",";
759  std::cout << z1;
760  std::cout << ")" << std::endl;
761  ;
762  return 0;
763  }
764 
765  double part = 0;
766  double G = 0;
767  for (int m = 0; m < NumberOfOrders; m++)
768  {
769  for (int n = 0; n < NumberOfOrders; n++)
770  {
771  double term = 1;
772  part = (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1)); //unitless
773  if (verbosity > 10) printf("(2 - ((m==0)?1:0))*cos(m*(phi-phi1)); = %f\n", part);
774  term *= part;
775  part = sin(BetaN[n] * z) * sin(BetaN[n] * z1); //unitless
776  if (verbosity > 10) printf("sin(BetaN[n]*z)*sin(BetaN[n]*z1); = %f\n", part);
777  term *= part;
778 
779  if (r < r1)
780  {
781  term *= RPrime(m, n, a, r) * Rmn2(m, n, r1); //units of 1/[L]
782  }
783  else
784  {
785  term *= Rmn1(m, n, r1) * RPrime(m, n, b, r); //units of 1/[L]
786  }
787  term /= bessel_denominator[m][n]; //unitless
788  G += term;
789  }
790  }
791 
792  G = G / (L * pi); //units of 1/[L] -- net is 1/[L]^2
793  if (verbosity) std::cout << "Er = " << G << std::endl;
794 
795  return G;
796 }
797 
798 double Rossegger::Er_(double r, double phi, double z, double r1, double phi1, double z1)
799 {
800  //doesn't take advantage of precalcs.
801  //as in Rossegger 5.65
802  //field at r, phi, z due to unit charge at r1, phi1, z1;
803  //if(fByFile && fabs(r-r1)>MinimumDR && fabs(z-z1)>MinimumDZ) return ByFileER(r,phi,z,r1,phi1,z1);
804  // Check input arguments for sanity...
805  int error = 0;
806  if (r < a || r > b) error = 1;
807  if (phi < 0 || phi > 2 * pi) error = 1;
808  if (z < 0 || z > L) error = 1;
809  if (r1 < a || r1 > b) error = 1;
810  if (phi1 < 0 || phi1 > 2 * pi) error = 1;
811  if (z1 < 0 || z1 > L) error = 1;
812  if (error)
813  {
814  std::cout << "Invalid arguments Er(";
815  std::cout << r << ",";
816  std::cout << phi << ",";
817  std::cout << z << ",";
818  std::cout << r1 << ",";
819  std::cout << phi1 << ",";
820  std::cout << z1;
821  std::cout << ")" << std::endl;
822  ;
823  return 0;
824  }
825 
826  double G = 0;
827  for (int m = 0; m < NumberOfOrders; m++)
828  {
829  for (int n = 0; n < NumberOfOrders; n++)
830  {
831  double term = 1 / (L * pi);
832  term *= (2 - ((m == 0) ? 1 : 0)) * cos(m * (phi - phi1));
833  double BetaN_ = (n + 1) * pi / L;
834  term *= sin(BetaN_ * z) * sin(BetaN_ * z1);
835  if (r < r1)
836  {
837  term *= RPrime_(m, n, a, r) * Rmn2_(m, n, r1);
838  }
839  else
840  {
841  term *= Rmn1_(m, n, r1) * RPrime_(m, n, b, r);
842  }
843 
844  term /= TMath::BesselI(m, BetaN_ * a) * TMath::BesselK(m, BetaN_ * b) - TMath::BesselI(m, BetaN_ * b) * TMath::BesselK(m, BetaN_ * a);
845 
846  G += term;
847  }
848  }
849 
850  if (verbosity) std::cout << "Er = " << G << std::endl;
851 
852  return G;
853 }
854 
855 double Rossegger::Ephi(double r, double phi, double z, double r1, double phi1, double z1)
856 {
857  //rcc streamlined Ephi term
858  //compute field at rphiz from charge at r1phi1z1
859  // Check input arguments for sanity...
860  int error = 0;
861  if (r < a || r > b) error = 1;
862  if (phi < 0 || phi > 2 * pi) error = 1;
863  if (z < 0 || z > L) error = 1;
864  if (r1 < a || r1 > b) error = 1;
865  if (phi1 < 0 || phi1 > 2 * pi) error = 1;
866  if (z1 < 0 || z1 > L) error = 1;
867  if (error)
868  {
869  std::cout << "Invalid arguments Ephi(";
870  std::cout << r << ",";
871  std::cout << phi << ",";
872  std::cout << z << ",";
873  std::cout << r1 << ",";
874  std::cout << phi1 << ",";
875  std::cout << z1;
876  std::cout << ")" << std::endl;
877  ;
878  return 0;
879  }
880 
881  double G = 0;
882  //Rossegger Eqn. 5.66:
883  for (int k = 0; k < NumberOfOrders; k++)
884  {
885  for (int n = 0; n < NumberOfOrders; n++)
886  {
887  double term = 1;
888  term *= sin(BetaN[n] * z) * sin(BetaN[n] * z1); //unitless
889  term *= Rnk(n, k, r) * Rnk(n, k, r1) / N2nk[n][k]; //unitless?
890 
891  //the derivative of cosh(munk(pi-|phi-phi1|)
892  if (phi > phi1)
893  {
894  term *= -sinh(Munk[n][k] * (pi - (phi - phi1))); //unitless
895  }
896  else
897  {
898  term *= sinh(Munk[n][k] * (pi - (phi1 - phi))); //unitless
899  }
900  term *= 1 / sinh_pi_Munk[n][k]; //unitless
901  G += term;
902  }
903  }
904 
905  G = G / (L * r); //units of 1/[L]^2. r comes from the phi term in cylindrical gradient expression.
906  if (verbosity) std::cout << "Ephi = " << G << std::endl;
907 
908  return G;
909 }
910 
911 double Rossegger::Ephi_(double r, double phi, double z, double r1, double phi1, double z1)
912 {
913  //compute field at rphiz from charge at r1phi1z1
914  // Check input arguments for sanity...
915  int error = 0;
916  if (r < a || r > b) error = 1;
917  if (phi < 0 || phi > 2 * pi) error = 1;
918  if (z < 0 || z > L) error = 1;
919  if (r1 < a || r1 > b) error = 1;
920  if (phi1 < 0 || phi1 > 2 * pi) error = 1;
921  if (z1 < 0 || z1 > L) error = 1;
922  if (error)
923  {
924  std::cout << "Invalid arguments Ephi(";
925  std::cout << r << ",";
926  std::cout << phi << ",";
927  std::cout << z << ",";
928  std::cout << r1 << ",";
929  std::cout << phi1 << ",";
930  std::cout << z1;
931  std::cout << ")" << std::endl;
932  ;
933  return 0;
934  }
935 
936  verbosity = 1;
937  double G = 0;
938  //Rossegger Eqn. 5.66:
939  for (int k = 0; k < NumberOfOrders; k++) //off by one from Rossegger convention!
940  {
941  if (verbosity) std::cout << "\nk=" << k;
942  for (int n = 0; n < NumberOfOrders; n++) //off by one from Rossegger convention!
943  {
944  if (verbosity) std::cout << " n=" << n;
945  double BetaN_ = (n + 1) * pi / L;
946  double term = 1 / (L * r);
947  if (verbosity) std::cout << " 1/L=" << term;
948  term *= sin(BetaN_ * z) * sin(BetaN_ * z1);
949  if (verbosity) std::cout << " *sinsin=" << term;
950  term *= Rnk_(n, k, r) * Rnk_(n, k, r1);
951  if (verbosity) std::cout << " *rnkrnk=" << term;
952  term /= N2nk[n][k];
953  if (verbosity) std::cout << " */nnknnk=" << term;
954 
955  //the derivative of cosh(munk(pi-|phi-phi1|)
956  if (phi > phi1)
957  {
958  term *= -sinh(Munk[n][k] * (pi - (phi - phi1)));
959  //term *= Munk[n][k]*sinh(Munk[n][k]*pi*(phi1-phi));
960  //this originally has a factor of Munk in front, but that cancels with one in the denominator
961  }
962  else
963  {
964  term *= sinh(Munk[n][k] * (pi - (phi1 - phi)));
965  //term *= -Munk[n][k]*sinh(Munk[n][k]*pi*(phi-phi1));
966  //this originally has a factor of Munk in front, but that cancels with one in the denominator
967  }
968  if (verbosity) std::cout << " *sinh(mu*pi-phi-phi)=" << term;
969  term *= 1 / (sinh(pi * Munk[n][k]));
970  //term *= 1/(Munk[n][k]*sinh(pi*Munk[n][k]));
971  //this originally has a factor of Munk in front, but that cancels with one in the numerator
972  G += term;
973  if (verbosity) std::cout << " /sinh=" << term << " G=" << G << std::endl;
974  }
975  }
976  if (verbosity) std::cout << "Ephi = " << G << std::endl;
977  verbosity = 0;
978 
979  return G;
980 }
981 
982 void Rossegger::SaveZeroes(const char *destfile)
983 {
984  TFile *output = TFile::Open(destfile, "RECREATE");
985  output->cd();
986 
987  TTree *tInfo = new TTree("info", "Mu[n][k] values");
988  int ord = NumberOfOrders;
989  tInfo->Branch("order", &ord);
990  tInfo->Branch("epsilon", &epsilon);
991  tInfo->Fill();
992 
993  int n, k, m;
994  double munk, betamn;
995  double n2nk, n2mn;
996  TTree *tmunk = new TTree("munk", "Mu[n][k] values");
997  tmunk->Branch("n", &n);
998  tmunk->Branch("k", &k);
999  tmunk->Branch("munk", &munk);
1000  tmunk->Branch("n2nk", &n2nk);
1001  for (n = 0; n < ord; n++)
1002  {
1003  for (k = 0; k < ord; k++)
1004  {
1005  munk = Munk[n][k];
1006  n2nk = N2nk[n][k];
1007  tmunk->Fill();
1008  }
1009  }
1010 
1011  TTree *tbetamn = new TTree("betamn", "Beta[m][n] values");
1012  tbetamn->Branch("m", &m);
1013  tbetamn->Branch("n", &n);
1014  tbetamn->Branch("betamn", &betamn);
1015  tbetamn->Branch("n2mn", &n2mn);
1016  for (m = 0; m < ord; m++)
1017  {
1018  for (n = 0; n < ord; n++)
1019  {
1020  betamn = Betamn[m][n];
1021  n2mn = N2mn[m][n];
1022  tbetamn->Fill();
1023  }
1024  }
1025 
1026  tInfo->Write();
1027  tmunk->Write();
1028  tbetamn->Write();
1029  //output->Write();
1030  output->Close();
1031  return;
1032 }
1033 
1034 void Rossegger::LoadZeroes(const char *destfile)
1035 {
1036  TFile *f = TFile::Open(destfile, "READ");
1037  printf("reading rossegger zeroes from %s\n", destfile);
1038  TTree *tInfo = (TTree *) (f->Get("info"));
1039  int ord;
1040  tInfo->SetBranchAddress("order", &ord);
1041  tInfo->SetBranchAddress("epsilon", &epsilon);
1042  tInfo->GetEntry(0);
1043  printf("order=%d,epsilon=%f\n", ord, epsilon);
1044 
1045  int n, k, m;
1046  double munk, betamn;
1047  double n2nk, n2mn;
1048  TTree *tmunk = (TTree *) (f->Get("munk"));
1049  tmunk->SetBranchAddress("n", &n);
1050  tmunk->SetBranchAddress("k", &k);
1051  tmunk->SetBranchAddress("munk", &munk);
1052  tmunk->SetBranchAddress("n2nk", &n2nk);
1053  for (int i = 0; i < tmunk->GetEntries(); i++)
1054  {
1055  tmunk->GetEntry(i);
1056  Munk[n][k] = munk;
1057  N2nk[n][k] = n2nk;
1058  }
1059 
1060  TTree *tbetamn = (TTree *) (f->Get("betamn"));
1061  tbetamn->SetBranchAddress("m", &m);
1062  tbetamn->SetBranchAddress("n", &n);
1063  tbetamn->SetBranchAddress("betamn", &betamn);
1064  tbetamn->SetBranchAddress("n2mn", &n2mn);
1065  for (int i = 0; i < tbetamn->GetEntries(); i++)
1066  {
1067  tbetamn->GetEntry(i);
1068  Betamn[m][n] = betamn;
1069  N2mn[m][n] = n2mn;
1070  }
1071 
1072  f->Close();
1073  return;
1074 }