Ignore:
Timestamp:
May 28, 2009, 4:26:57 PM (15 years ago)
Author:
garnier
Message:

maj sur la beta de geant 4.9.3

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QBesIKJY.cc

    r1007 r1055  
    2525//
    2626//
    27 // $Id: G4QBesIKJY.cc,v 1.2 2006/06/29 20:06:45 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02 $
     27// $Id: G4QBesIKJY.cc,v 1.3 2009/02/23 09:49:24 mkossov Exp $
     28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2929//
    3030//      ---------------- G4QBesIKJY ----------------
    … …  
    3232//      class  for Bessel I0/I1 and K0/K1 functions in CHIPS Model
    3333// -------------------------------------------------------------------
    34 
     34// Short description: Bessel functions class (can be substituted)
     35// -------------------------------------------------------------------
    3536//#define debug
    3637//#define pdebug
    … …  
    4344  ex=false;
    4445  switch (type)
    45                 {
    46                   case BessI0:
     46  {
     47    case BessI0:
    4748      nu=0;
    4849      ik=true;
    4950      ij=true;
    50                                                 break;
    51                   case BessI1:
     51      break;
     52    case BessI1:
    5253      nu=1;
    5354      ik=true;
    5455      ij=true;
    55                                                 break;
    56                   case EBessI0:
     56      break;
     57    case EBessI0:
    5758      nu=0;
    5859      ex=true;
    5960      ik=true;
    6061      ij=true;
    61                                                 break;
    62                   case EBessI1:
     62      break;
     63    case EBessI1:
    6364      nu=1;
    6465      ex=true;
    6566      ik=true;
    6667      ij=true;
    67                                                 break;
    68                   case BessJ0:
     68      break;
     69    case BessJ0:
    6970      nu=0;
    7071      ik=true;
    7172      ij=false;
    72                                                 break;
    73                   case BessJ1:
     73      break;
     74    case BessJ1:
    7475      nu=1;
    7576      ik=true;
    7677      ij=false;
    77                                                 break;
    78                   case BessK0:
     78      break;
     79    case BessK0:
    7980      nu=0;
    8081      ik=false;
    8182      ij=true;
    82                                                 break;
    83                   case BessK1:
     83      break;
     84    case BessK1:
    8485      nu=1;
    8586      ik=false;
    8687      ij=true;
    87                                                 break;
    88                   case EBessK0:
     88      break;
     89    case EBessK0:
    8990      nu=0;
    9091      ex=true;
    9192      ik=false;
    9293      ij=true;
    93                                                 break;
    94                   case EBessK1:
     94      break;
     95    case EBessK1:
    9596      nu=1;
    9697      ex=true;
    9798      ik=false;
    9899      ij=true;
    99                                                 break;
    100                   case BessY0:
     100      break;
     101    case BessY0:
    101102      nu=0;
    102103      ik=false;
    103104      ij=false;
    104                                                 break;
    105                   case BessY1:
     105      break;
     106    case BessY1:
    106107      nu=1;
    107108      ik=false;
    108109      ij=false;
    109                                                 break;
     110      break;
    110111  }
    111112}
    … …  
    126127  static const G4double EPS = 1.E-15;
    127128  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;
     129  static const G4double HF  = Z1/2;
     130  static const G4double R8  = HF/4;
     131  static const G4double R32 = R8/4;
    131132  static const G4double PI  = 3.14159265358979324;
    132133  static const G4double CE  = 0.57721566490153286;
    … …  
    139140  static const G4double CI0[npi]={+1.00829205458740032E0,
    140141             +.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};
     142             +.00000513587726878E0,+.00000056816965808E0,+.00000008513091223E0,
     143             +.00000001238425364E0,+.00000000029801672E0,-.00000000078956698E0,
     144             -.00000000033127128E0,-.00000000004497339E0,+.00000000001799790E0,
     145             +.00000000000965748E0,+.00000000000038604E0,-.00000000000104039E0,
     146             -.00000000000023950E0,+.00000000000009554E0,+.00000000000004443E0,
     147             -.00000000000000859E0,-.00000000000000709E0,+.00000000000000087E0,
     148             +.00000000000000112E0,-.00000000000000012E0,-.00000000000000018E0};
    148149
    149150  static const G4double CI1[npi]={+.975800602326285926E0,
    … …  
    215216
    216217  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};
     218          -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0,
     219          +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0,
     220          -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0,
     221          +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0,
     222          -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0};
    222223
    223224  static const G4complex CF[npj]={
    … …  
    244245  // -------------------------------------------------------------------------------------
    245246  G4double H=0.;                     // Prototype of the result
    246                 if (ij)                            // I/K Bessel functions
     247  if (ij)                            // I/K Bessel functions
    247248  {
    248249    if (ik)                          // I0/I1/EI0/EI1 Bessel functions (symmetric)
    … …  
    263264        G4double A1=1.+DY/(XLI*V3);
    264265        G4double A2=1.+Y*(4.+(DY+Y)/(XLD*XL))/(W1*V3);
    265                                     G4double B0=1.;
    266                                     G4double    B1=1.-Y/XLI;
     266        G4double B0=1.;
     267        G4double B1=1.-Y/XLI;
    267268        G4double B2=1.-Y*(1.-Y/(XLD+XLD))/W1;
    268269        G4int    V1=3-XL;
    269270        G4double V2=V3+V3;
    270271        G4double C=0.;
    271                                     for (G4int N=3; N<=30; N++)
     272        for (G4int N=3; N<=30; N++)
    272273        {
    273274          G4double C0=C;
    … …  
    279280          G4int    W5=W4-1;
    280281          G4int    W6=W5-1;
    281                                                              V1=V1+1;
     282                   V1=V1+1;
    282283                   V2=V2+1;
    283284                   V3=V3+1;
    284                                                     G4double U1=FN*W4;
     285          G4double U1=FN*W4;
    285286          G4double E=V3/(U1*W3);
    286                                                     G4double U2=E*Y;
     287          G4double U2=E*Y;
    287288          G4double F1=1.+Y*V1/(U1*W1);
    288289          G4double F2=(1.+Y*V2/(V3*W2*W5))*U2;
    … …  
    290291          G4double A=F1*A2+F2*A1+F3*A0;
    291292          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;
     293                   C=A/B;
     294          if (std::abs(C0-C) < EPS*std::abs(C)) break;
     295          A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
    295296        }
    296297        H=C;
    … …  
    301302      {
    302303        G4double P=16./V-1.;
    303                                     G4double ALFA=P+P;
    304                                     G4double B1=0.;
     304        G4double ALFA=P+P;
     305        G4double B1=0.;
    305306        G4double B2=0.;
    306                                                   for (G4int I=npil; I>=0; I--)
     307        for (G4int I=npil; I>=0; I--)
    307308        {
    308309          if (!nu) CJ=CI0[I];
    309310          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                                   }
     311          G4double B0=CJ+ALFA*B1-B2;
     312                   B2=B1;
     313                   B1=B0;
     314        }
     315        H=std::sqrt(RPI2/V)*(B1-P*B2);
     316        if (nu && X < 0.) H=-H;
     317        if (!ex) H*=std::exp(V);
     318      }
    318319    }
    319320    else                             // K0/K1/EK0/EK1 Bessel functions
    320                   {
     321    {
    321322#ifdef debug
    322323      G4cout<<"G4BesIKJY: >>>>>>>>>>>>>> K is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
    … …  
    339340        G4double Q=HF;
    340341        G4double C=1.;
    341                                                   G4double D=B*B;
     342        G4double D=B*B;
    342343        G4double BK1=P;
    343344        for (G4int N=1; N<=nat1; N++)  // @@ "nat" can be increased
    … …  
    349350                   C*=D/FN;
    350351          G4double G=C*(P-FN*F);
    351                         G4double R=C*F;
     352          G4double R=C*F;
    352353                   BK=BK+R;
    353354                   BK1=BK1+G;
    354                                                                   if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break;
     355          if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break;
    355356        }
    356357        if (nu==1) H=BK1/B;
    … …  
    358359        if (ex) H*=std::exp(X);
    359360      }
    360                                   else if (X < 5.)
     361      else if (X < 5.)
    361362      {
    362363#ifdef debug
    … …  
    368369        G4double XN=DNUS+DNUS;
    369370        G4double A=9.-XN;
    370                                                   G4double B=25.-XN;
     371        G4double B=25.-XN;
    371372        G4double C=768*X*X;
    372                                                   G4double HX=16*X;
    373                                                   G4double C0=HX+HX+HX;;
     373        G4double HX=16*X;
     374        G4double C0=HX+HX+HX;;
    374375        G4double A0=1.;
    375376        G4double A1=(HX+7.+XN)/A;
    … …  
    378379        G4double B1=(HX+9.-XN)/A;
    379380        G4double B2=(C+C0*B)/(A*B)+1.;
    380                                                            C=0.;
    381                                                   for (G4int N=3; N<=nat2; N++)
     381                 C=0.;
     382        for (G4int N=3; N<=nat2; N++)
    382383        {
    383384          C0=C;
    … …  
    389390          G4double FN3=FN1/(FN2-3.);
    390391          G4double FN4=12*FN*FN-(1.-XN);
    391                                                                   G4double FN5=16*FN1*X;
     392          G4double FN5=16*FN1*X;
    392393          G4double RAN=1./(FNP*FNP-XN);
    393                                                                   G4double F1=FN3*(FN4-20*FN)+FN5;
     394          G4double F1=FN3*(FN4-20*FN)+FN5;
    394395          G4double F2=28*FN-FN4-8.+FN5;
    395396          G4double F3=FN3*(FNM*FNM-XN);
    396                                                                                         A=(F1*A2+F2*A1+F3*A0)*RAN;
     397                   A=(F1*A2+F2*A1+F3*A0)*RAN;
    397398                   B=(F1*B2+F2*B1+F3*B0)*RAN;
    398399                   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);
     400          if (std::abs(C0-C) < EPS*std::abs(C)) break;
     401          A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
     402        }
     403        H=C/std::sqrt(RPIH*X);
    403404        if (!ex) H*=std::exp(-X);
    404                                   }
     405      }
    405406      else
    406407      {
    … …  
    410411        G4double P=10./X-1.;
    411412        G4double ALFA=P+P;
    412                                                         G4double B1=0.;
    413                                                         G4double B2=0.;
     413        G4double B1=0.;
     414        G4double B2=0.;
    414415#ifdef debug
    415416        G4cout<<"G4BesIKJY: >>> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
    416417#endif
    417                                                   for (G4int I=npkl; I>=0; I--)
     418        for (G4int I=npkl; I>=0; I--)
    418419        {
    419420          if (!nu) CK=CK0[I];
    … …  
    422423                   B2=B1;
    423424                   B1=B0;
    424                                                   }
     425        }
    425426        H=std::sqrt(PIH/X)*(B1-P*B2);
    426427        if (!ex) H*=std::exp(-X);
    427428      }
    428                   }
     429    }
    429430  }
    430431  else
    431                 {
     432  {
    432433    if (!ik && X < 0.)
    433434    {
    … …  
    455456          B1=0.;
    456457          B2=0.;
    457                                                                                 for (G4int JT=npkl; JT>=0; JT--)
     458          for (G4int JT=npkl; JT>=0; JT--)
    458459          {
    459460            G4double B0=CB[JT]+ALFA*B1-B2;
    … …  
    463464          H=RPIH*(CE+std::log(HF*X))*H+B1-P*B2;
    464465        }
    465                                                 }
     466      }
    466467      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--)
     468      {
     469        G4double P=10./V-1.;
     470        G4double ALFA=P+P;
     471        G4complex CB1(0.,0.);
     472        G4complex CB2(0.,0.);
     473        for (G4int IT=npjl; IT>=0; IT--)
    473474        {
    474475          G4complex CB0=CC[IT]+ALFA*CB1-CB2;
    475476                    CB2=CB1;
    476477                    CB1=CB0;
    477                                                                 }
     478        }
    478479        CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI4))*(CB1-P*CB2);
    479                                                                 if (ik) H=real(CB1);
     480        if (ik) H=real(CB1);
    480481        else    H=real(-CI*CB1);
    481                                                 }
     482      }
    482483    }
    483484    else                          // J1/Y1 Bessel functions
    … …  
    494495        {
    495496          G4double B0=CD[IT]+ALFA*B1-B2;
    496                                                                                          B2=B1;
     497                   B2=B1;
    497498                   B1=B0;
    498499        }
    … …  
    510511          H=RPIH*((CE+std::log(HF*X))*H-1./X)+Y*(B1-B2);
    511512        }
    512                                                 }
     513      }
    513514      else
    514                                                 {
     515      {
    515516        G4double P=10./V-1.;
    516517        G4double ALFA=P+P;
    … …  
    520521        {
    521522          G4complex CB0=CF[IT]+ALFA*CB1-CB2;
    522                                                                                           CB2=CB1;
     523                    CB2=CB1;
    523524                    CB1=CB0;
    524525        }
    … …  
    526527        if (ik) H=real(CB1);
    527528        else    H=real(-CI*CB1);
    528                                                 }
     529      }
    529530      if (X < 0.) H=-H;
    530531    }
Note: See TracChangeset for help on using the changeset viewer.