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

Last change on this file since 967 was 962, checked in by garnier, 15 years ago

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