source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QBesIKJY.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 18.8 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4QBesIKJY.cc,v 1.4 2009/11/10 17:13:46 mkossov Exp $
28// GEANT4 tag $Name: hadr-chips-V09-03-08 $
29//
30//      ---------------- G4QBesIKJY ----------------
31//             by Mikhail Kossov, Sept 1999.
32//      class  for Bessel I0/I1 and K0/K1 functions in CHIPS Model
33// -------------------------------------------------------------------
34// Short description: Bessel functions class (can be substituted)
35// -------------------------------------------------------------------
36
37//#define debug
38//#define pdebug
39
40#include "G4QBesIKJY.hh"
41
42// Constructor
43G4QBesIKJY::G4QBesIKJY(G4QBIType type)
44{
45  ex=false;
46  switch (type)
47  {
48    case BessI0:
49      nu=0;
50      ik=true;
51      ij=true;
52      break;
53    case BessI1:
54      nu=1;
55      ik=true;
56      ij=true;
57      break;
58    case EBessI0:
59      nu=0;
60      ex=true;
61      ik=true;
62      ij=true;
63      break;
64    case EBessI1:
65      nu=1;
66      ex=true;
67      ik=true;
68      ij=true;
69      break;
70    case BessJ0:
71      nu=0;
72      ik=true;
73      ij=false;
74      break;
75    case BessJ1:
76      nu=1;
77      ik=true;
78      ij=false;
79      break;
80    case BessK0:
81      nu=0;
82      ik=false;
83      ij=true;
84      break;
85    case BessK1:
86      nu=1;
87      ik=false;
88      ij=true;
89      break;
90    case EBessK0:
91      nu=0;
92      ex=true;
93      ik=false;
94      ij=true;
95      break;
96    case EBessK1:
97      nu=1;
98      ex=true;
99      ik=false;
100      ij=true;
101      break;
102    case BessY0:
103      nu=0;
104      ik=false;
105      ij=false;
106      break;
107    case BessY1:
108      nu=1;
109      ik=false;
110      ij=false;
111      break;
112  }
113}
114
115G4QBesIKJY::~G4QBesIKJY(){;}
116
117G4double G4QBesIKJY::operator() (G4double X) const
118{
119  static const G4int nat1 = 15;            // a # of attempts to reach the X<1 accuracy
120  static const G4int nat2 = nat1+nat1;     // a # of attempts to reach the X<5 accuracy
121  static const G4int npi = 25;
122  static const G4int npil = npi-1;
123  static const G4int npk = 17;
124  static const G4int npkl = npk-1;
125  static const G4int npj = 20;
126  static const G4int npjl = npj-1;
127  static const G4complex CI(0,1);
128  static const G4double EPS = 1.E-15;
129  static const G4double Z1  = 1.;
130  static const G4double HF  = Z1/2;
131  static const G4double R8  = HF/4;
132  static const G4double R32 = R8/4;
133  static const G4double PI  = 3.14159265358979324;
134  static const G4double CE  = 0.57721566490153286;
135  static const G4double PIH = PI/2;
136  static const G4double PI4 = PIH/2;   // PI/4
137  static const G4double PI3 = PIH+PI4; // 3*PI/4
138  static const G4double RPIH = 2./PI;
139  static const G4double RPI2 = RPIH/4;
140
141  static const G4double CI0[npi]={+1.00829205458740032E0,
142             +.00845122624920943E0,+.00012700630777567E0,+.00007247591099959E0,
143             +.00000513587726878E0,+.00000056816965808E0,+.00000008513091223E0,
144             +.00000001238425364E0,+.00000000029801672E0,-.00000000078956698E0,
145             -.00000000033127128E0,-.00000000004497339E0,+.00000000001799790E0,
146             +.00000000000965748E0,+.00000000000038604E0,-.00000000000104039E0,
147             -.00000000000023950E0,+.00000000000009554E0,+.00000000000004443E0,
148             -.00000000000000859E0,-.00000000000000709E0,+.00000000000000087E0,
149             +.00000000000000112E0,-.00000000000000012E0,-.00000000000000018E0};
150
151  static const G4double CI1[npi]={+.975800602326285926E0,
152           -.024467442963276385E0,-.000277205360763829E0,-.000009732146728020E0,
153           -.000000629724238640E0,-.000000065961142154E0,-.000000009613872919E0,
154           -.000000001401140901E0,-.000000000047563167E0,+.000000000081530681E0,
155           +.000000000035408148E0,+.000000000005102564E0,-.000000000001804409E0,
156           -.000000000001023594E0,-.000000000000052678E0,+.000000000000107094E0,
157           +.000000000000026120E0,-.000000000000009561E0,-.000000000000004713E0,
158           +.000000000000000829E0,+.000000000000000743E0,-.000000000000000080E0,
159           -.000000000000000117E0,+.000000000000000011E0,+.000000000000000019E0};
160
161  static const G4double CK0[npk]={+.988408174230825800E0,-.011310504646928281E0,
162           +.000269532612762724E0,-.000011106685196665E0,+.000000632575108500E0,
163           -.000000045047337641E0,+.000000003792996456E0,-.000000000364547179E0,
164           +.000000000039043756E0,-.000000000004579936E0,+.000000000000580811E0, 
165           -.000000000000078832E0,+.000000000000011360E0,-.000000000000001727E0,
166           +.000000000000000275E0,-.000000000000000046E0,+.000000000000000008E0};
167
168  static const G4double CK1[npk]={+1.035950858772358331E0,+.035465291243331114E0,
169            -.000468475028166889E0,+.000016185063810053E0,-.000000845172048124E0,
170            +.000000057132218103E0,-.000000004645554607E0,+.000000000435417339E0,
171            -.000000000045757297E0,+.000000000005288133E0,-.000000000000662613E0,
172            +.000000000000089048E0,-.000000000000012726E0,+.000000000000001921E0,
173            -.000000000000000305E0,+.000000000000000050E0,-.000000000000000009E0};
174
175  static const G4double CA[npk]={+.157727971474890120E0,-.008723442352852221E0,
176          +.265178613203336810E0,-.370094993872649779E0,+.158067102332097261E0,
177          -.034893769411408885E0,+.004819180069467605E0,-.000460626166206275E0,
178          +.000032460328821005E0,-.000001761946907762E0,+.000000076081635924E0,
179          -.000000002679253531E0,+.000000000078486963E0,-.000000000001943835E0,
180          +.000000000000041253E0,-.000000000000000759E0,+.000000000000000012E0};
181
182  static const G4double CB[npk]={-.021505111449657551E0,-.275118133043518791E0,
183          +.198605634702554156E0,+.234252746109021802E0,-.165635981713650413E0,
184          +.044621379540669282E0,-.006932286291523188E0,+.000719117403752303E0,
185          -.000053925079722939E0,+.000003076493288108E0,-.000000138457181230E0,
186          +.000000005051054369E0,-.000000000152582850E0,+.000000000003882867E0,
187          -.000000000000084429E0,+.000000000000001587E0,-.000000000000000026E0};
188
189  static const G4complex CC[npj]={
190    G4complex(+.998988089858965153E0,-.012331520578544144E0),
191    G4complex(-.001338428549971856E0,-.012249496281259475E0),
192    G4complex(-.000318789878061893E0,+.000096494184993423E0),
193    G4complex(+.000008511232210657E0,+.000013655570490357E0),
194    G4complex(+.000000691542349139E0,-.000000851806644426E0),
195    G4complex(-.000000090770101537E0,-.000000027244053414E0),
196    G4complex(+.000000001454928079E0,+.000000009646421338E0),
197    G4complex(+.000000000926762487E0,-.000000000683347518E0),
198    G4complex(-.000000000139166198E0,-.000000000060627380E0),
199    G4complex(+.000000000003237975E0,+.000000000021695716E0),
200    G4complex(+.000000000002535357E0,-.000000000002304899E0),
201    G4complex(-.000000000000559090E0,-.000000000000122554E0),
202    G4complex(+.000000000000041919E0,+.000000000000092314E0),
203    G4complex(+.000000000000008733E0,-.000000000000016778E0),
204    G4complex(-.000000000000003619E0,+.000000000000000754E0),
205    G4complex(+.000000000000000594E0,+.000000000000000462E0),
206    G4complex(-.000000000000000010E0,-.000000000000000159E0),
207    G4complex(-.000000000000000024E0,+.000000000000000025E0),
208    G4complex(+.000000000000000008E0,+.000000000000000000E0),
209    G4complex(-.000000000000000001E0,-.000000000000000001E0)};
210
211  static const G4double CD[npk]={+0.648358770605264921E0,-1.191801160541216873E0,
212         +1.287994098857677620E0,-0.661443934134543253E0,+0.177709117239728283E0,
213         -0.029175524806154208E0,+0.003240270182683857E0,-0.000260444389348581E0,
214         +0.000015887019239932E0,-0.000000761758780540E0,+0.000000029497070073E0,
215         -0.000000000942421298E0,+0.000000000025281237E0,-0.000000000000577740E0,
216         +0.000000000000011386E0,-0.000000000000000196E0,+0.000000000000000003E0};
217
218  static const G4double EE[npk]={-.040172946544414076E0,-.444447147630558063E0,
219          -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0,
220          +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0,
221          -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0,
222          +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0,
223          -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0};
224
225  static const G4complex CF[npj]={
226    G4complex(+1.001702234853820996E0,+.037261715000537654E0),
227     G4complex(+.002255572846561180E0,+.037145322479807690E0),
228     G4complex(+.000543216487508013E0,-.000137263238201907E0),
229     G4complex(-.000011179461895408E0,-.000019851294687597E0),
230     G4complex(-.000000946901382392E0,+.000001070014057386E0),
231     G4complex(+.000000111032677121E0,+.000000038305261714E0),
232     G4complex(-.000000001294398927E0,-.000000011628723277E0),
233     G4complex(-.000000001114905944E0,+.000000000759733092E0),
234     G4complex(+.000000000157637232E0,+.000000000075476075E0),
235     G4complex(-.000000000002830457E0,-.000000000024752781E0),
236     G4complex(-.000000000002932169E0,+.000000000002493893E0),
237     G4complex(+.000000000000617809E0,+.000000000000156198E0),
238     G4complex(-.000000000000043162E0,-.000000000000103385E0),
239     G4complex(-.000000000000010133E0,+.000000000000018129E0),
240     G4complex(+.000000000000003973E0,-.000000000000000709E0),
241     G4complex(-.000000000000000632E0,-.000000000000000520E0),
242     G4complex(+.000000000000000006E0,+.000000000000000172E0),
243     G4complex(+.000000000000000027E0,-.000000000000000026E0),
244     G4complex(-.000000000000000008E0,-.000000000000000000E0),
245     G4complex(+.000000000000000001E0,+.000000000000000001E0)};
246  // -------------------------------------------------------------------------------------
247  G4double H=0.;                     // Prototype of the result
248  if (ij)                            // I/K Bessel functions
249  {
250    if (ik)                          // I0/I1/EI0/EI1 Bessel functions (symmetric)
251    {
252      G4double V=std::abs(X);
253      G4double CJ=0.;                // Prototype of the element of the CI0/CI1 matrix
254      if (V < 8.)
255      {
256        G4double HFV=HF*V;
257        G4double Y=HFV*HFV;
258        G4int    V3=nu+1;
259        G4int    XL=V3+1;
260        G4int    XLI=XL+1;
261        G4int    XLD=XLI+1;
262        G4int    W1=XLD+1;
263        G4double A0=1.;
264        G4double DY=Y+Y;
265        G4double A1=1.+DY/(XLI*V3);
266        G4double A2=1.+Y*(4.+(DY+Y)/(XLD*XL))/(W1*V3);
267        G4double B0=1.;
268        G4double B1=1.-Y/XLI;
269        G4double B2=1.-Y*(1.-Y/(XLD+XLD))/W1;
270        G4int    V1=3-XL;
271        G4double V2=V3+V3;
272        G4double C=0.;
273        for (G4int N=3; N<=30; N++)
274        {
275          G4double C0=C;
276          G4double FN=N;
277                   W1=W1+2;
278          G4int    W2=W1-1;
279          G4int    W3=W2-1;
280          G4int    W4=W3-1;
281          G4int    W5=W4-1;
282          G4int    W6=W5-1;
283                   V1=V1+1;
284                   V2=V2+1;
285                   V3=V3+1;
286          G4double U1=FN*W4;
287          G4double E=V3/(U1*W3);
288          G4double U2=E*Y;
289          G4double F1=1.+Y*V1/(U1*W1);
290          G4double F2=(1.+Y*V2/(V3*W2*W5))*U2;
291          G4double F3=-Y*Y*U2/(W4*W5*W5*W6);
292          G4double A=F1*A2+F2*A1+F3*A0;
293          G4double B=F1*B2+F2*B1+F3*B0;
294                   C=A/B;
295          if (std::abs(C0-C) < EPS*std::abs(C)) break;
296          A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
297        }
298        H=C;
299        if (nu==1) H*=HF*X;
300        if (ex) H*=std::exp(-V);
301      }
302      else
303      {
304        G4double P=16./V-1.;
305        G4double ALFA=P+P;
306        G4double B1=0.;
307        G4double B2=0.;
308        for (G4int I=npil; I>=0; I--)
309        {
310          if (!nu) CJ=CI0[I];
311          else     CJ=CI1[I];
312          G4double B0=CJ+ALFA*B1-B2;
313                   B2=B1;
314                   B1=B0;
315        }
316        H=std::sqrt(RPI2/V)*(B1-P*B2);
317        if (nu && X < 0.) H=-H;
318        if (!ex) H*=std::exp(V);
319      }
320    }
321    else                             // K0/K1/EK0/EK1 Bessel functions
322    {
323#ifdef debug
324      G4cout<<"G4BesIKJY: >>>>>>>>>>>>>> K is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
325#endif
326      G4double CK=0.;                // Prototype of the element of the CI0/CI1 matrix
327      if (X < 0.)
328      {
329        G4cout<<"G4BesIKJY::NegativeArg in K-BessFun X="<<X<<", n="<<nu<<",E="<<ex<<G4endl;
330        return H;
331      }
332      else if (X < 1.)
333      {
334#ifdef debug
335        G4cout<<"G4BesIKJY: >>>> [ X < 1 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
336#endif
337        G4double B=HF*X;
338        G4double BK=-(std::log(B)+CE);
339        G4double F=BK;
340        G4double P=HF;
341        G4double Q=HF;
342        G4double C=1.;
343        G4double D=B*B;
344        G4double BK1=P;
345        for (G4int N=1; N<=nat1; N++)  // @@ "nat" can be increased
346        {
347          G4double FN=N;
348                   P/=FN;
349                   Q/=FN;
350                   F=(F+P+Q)/FN;
351                   C*=D/FN;
352          G4double G=C*(P-FN*F);
353          G4double R=C*F;
354                   BK=BK+R;
355                   BK1=BK1+G;
356          if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break;
357        }
358        if (nu==1) H=BK1/B;
359        else       H=BK;
360        if (ex) H*=std::exp(X);
361      }
362      else if (X < 5.)
363      {
364#ifdef debug
365        G4cout<<"G4BesIKJY: >>>> [ X < 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
366#endif
367        G4int NUS=0;              // @@ NU**2 for future NU>1 applications
368        if (nu==1) NUS=1;
369        G4double DNUS=NUS+NUS;
370        G4double XN=DNUS+DNUS;
371        G4double A=9.-XN;
372        G4double B=25.-XN;
373        G4double C=768*X*X;
374        G4double HX=16*X;
375        G4double C0=HX+HX+HX;;
376        G4double A0=1.;
377        G4double A1=(HX+7.+XN)/A;
378        G4double A2=(C+C0*(XN+23.)+XN*(XN+62.)+129.)/(A*B);
379        G4double B0=1.;
380        G4double B1=(HX+9.-XN)/A;
381        G4double B2=(C+C0*B)/(A*B)+1.;
382                 C=0.;
383        for (G4int N=3; N<=nat2; N++)
384        {
385          C0=C;
386          G4double FN=N;
387          G4double FN2=FN+FN;
388          G4double FNP=FN2+1.;
389          G4double FN1=FN2-1.;
390          G4double FNM=FN1-4.;
391          G4double FN3=FN1/(FN2-3.);
392          G4double FN4=12*FN*FN-(1.-XN);
393          G4double FN5=16*FN1*X;
394          G4double RAN=1./(FNP*FNP-XN);
395          G4double F1=FN3*(FN4-20*FN)+FN5;
396          G4double F2=28*FN-FN4-8.+FN5;
397          G4double F3=FN3*(FNM*FNM-XN);
398                   A=(F1*A2+F2*A1+F3*A0)*RAN;
399                   B=(F1*B2+F2*B1+F3*B0)*RAN;
400                   C=A/B;
401          if (std::abs(C0-C) < EPS*std::abs(C)) break;
402          A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
403        }
404        H=C/std::sqrt(RPIH*X);
405        if (!ex) H*=std::exp(-X);
406      }
407      else
408      {
409#ifdef debug
410        G4cout<<"G4BesIKJY: >>> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
411#endif
412        G4double P=10./X-1.;
413        G4double ALFA=P+P;
414        G4double B1=0.;
415        G4double B2=0.;
416#ifdef debug
417        G4cout<<"G4BesIKJY: >>> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
418#endif
419        for (G4int I=npkl; I>=0; I--)
420        {
421          if (!nu) CK=CK0[I];
422          else     CK=CK1[I];
423          G4double B0=CK+ALFA*B1-B2;
424                   B2=B1;
425                   B1=B0;
426        }
427        H=std::sqrt(PIH/X)*(B1-P*B2);
428        if (!ex) H*=std::exp(-X);
429      }
430    }
431  }
432  else
433  {
434    if (!ik && X < 0.)
435    {
436      G4cout<<"G4BesIKJY::NegativeArgument in Y BesselFunction X="<<X<<", nu="<<nu<<G4endl;
437      return H;
438    }
439    G4double V=std::abs(X);
440    if (!nu)                          // J0/Y0 Bessel functions
441    {
442      if (V < 8.)
443      {
444        G4double P=R32*V*V-1.;
445        G4double ALFA=P+P;
446        G4double B1=0.;
447        G4double B2=0.;
448        for (G4int IT=npkl; IT>=0; IT--)
449        {
450          G4double B0=CA[IT]+ALFA*B1-B2;
451                   B2=B1;
452                   B1=B0;
453        }
454        H=B1-P*B2;
455        if (!ik)
456        {
457          B1=0.;
458          B2=0.;
459          for (G4int JT=npkl; JT>=0; JT--)
460          {
461            G4double B0=CB[JT]+ALFA*B1-B2;
462                     B2=B1;
463                     B1=B0;
464          }
465          H=RPIH*(CE+std::log(HF*X))*H+B1-P*B2;
466        }
467      }
468      else
469      {
470        G4double P=10./V-1.;
471        G4double ALFA=P+P;
472        G4complex CB1(0.,0.);
473        G4complex CB2(0.,0.);
474        for (G4int IT=npjl; IT>=0; IT--)
475        {
476          G4complex CB0=CC[IT]+ALFA*CB1-CB2;
477                    CB2=CB1;
478                    CB1=CB0;
479        }
480        CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI4))*(CB1-P*CB2);
481        if (ik) H=real(CB1);
482        else    H=real(-CI*CB1);
483      }
484    }
485    else                          // J1/Y1 Bessel functions
486    {
487      if (V < 8.)
488      {
489        G4double Y=R8*V;
490        G4double Y2=Y*Y;
491        G4double P=Y2+Y2-1.;
492        G4double ALFA=P+P;
493        G4double B1=0.;
494        G4double B2=0.;
495        for (G4int IT=npkl; IT>=0; IT--)
496        {
497          G4double B0=CD[IT]+ALFA*B1-B2;
498                   B2=B1;
499                   B1=B0;
500        }
501        H=Y*(B1-P*B2);
502        if (!ik)
503        {
504          B1=0.;
505          B2=0.;
506          for (G4int JT=npkl; JT>=0; JT--)
507          {
508            G4double B0=EE[JT]+ALFA*B1-B2;
509                     B2=B1;
510                     B1=B0;
511          }
512          H=RPIH*((CE+std::log(HF*X))*H-1./X)+Y*(B1-B2);
513        }
514      }
515      else
516      {
517        G4double P=10./V-1.;
518        G4double ALFA=P+P;
519        G4complex CB1(0.,0.);
520        G4complex CB2(0.,0.);
521        for (G4int IT=npjl; IT>=0; IT--)
522        {
523          G4complex CB0=CF[IT]+ALFA*CB1-CB2;
524                    CB2=CB1;
525                    CB1=CB0;
526        }
527        CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI3))*(CB1-P*CB2);
528        if (ik) H=real(CB1);
529        else    H=real(-CI*CB1);
530      }
531      if (X < 0.) H=-H;
532    }
533  }
534  return H;
535}
Note: See TracBrowser for help on using the repository browser.