Index: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QBesIKJY.cc
===================================================================
--- trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QBesIKJY.cc	(revision 1007)
+++ trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QBesIKJY.cc	(revision 1055)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4QBesIKJY.cc,v 1.2 2006/06/29 20:06:45 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4QBesIKJY.cc,v 1.3 2009/02/23 09:49:24 mkossov Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 //      ---------------- G4QBesIKJY ----------------
@@ -32,5 +32,6 @@
 //      class  for Bessel I0/I1 and K0/K1 functions in CHIPS Model
 // -------------------------------------------------------------------
-
+// Short description: Bessel functions class (can be substituted)
+// -------------------------------------------------------------------
 //#define debug
 //#define pdebug
@@ -43,69 +44,69 @@
   ex=false;
   switch (type)
-		{
-		  case BessI0:
+  {
+    case BessI0:
       nu=0;
       ik=true;
       ij=true;
-						break;
-		  case BessI1:
+      break;
+    case BessI1:
       nu=1;
       ik=true;
       ij=true;
-						break;
-		  case EBessI0:
+      break;
+    case EBessI0:
       nu=0;
       ex=true;
       ik=true;
       ij=true;
-						break;
-		  case EBessI1:
+      break;
+    case EBessI1:
       nu=1;
       ex=true;
       ik=true;
       ij=true;
-						break;
-		  case BessJ0:
+      break;
+    case BessJ0:
       nu=0;
       ik=true;
       ij=false;
-						break;
-		  case BessJ1:
+      break;
+    case BessJ1:
       nu=1;
       ik=true;
       ij=false;
-						break;
-		  case BessK0:
+      break;
+    case BessK0:
       nu=0;
       ik=false;
       ij=true;
-						break;
-		  case BessK1:
+      break;
+    case BessK1:
       nu=1;
       ik=false;
       ij=true;
-						break;
-		  case EBessK0:
+      break;
+    case EBessK0:
       nu=0;
       ex=true;
       ik=false;
       ij=true;
-						break;
-		  case EBessK1:
+      break;
+    case EBessK1:
       nu=1;
       ex=true;
       ik=false;
       ij=true;
-						break;
-		  case BessY0:
+      break;
+    case BessY0:
       nu=0;
       ik=false;
       ij=false;
-						break;
-		  case BessY1:
+      break;
+    case BessY1:
       nu=1;
       ik=false;
       ij=false;
-						break;
+      break;
   }
 }
@@ -126,7 +127,7 @@
   static const G4double EPS = 1.E-15;
   static const G4double Z1  = 1.;
-		static const G4double HF  = Z1/2;
-		static const G4double R8  = HF/4;
-		static const G4double R32 = R8/4;
+  static const G4double HF  = Z1/2;
+  static const G4double R8  = HF/4;
+  static const G4double R32 = R8/4;
   static const G4double PI  = 3.14159265358979324;
   static const G4double CE  = 0.57721566490153286;
@@ -139,11 +140,11 @@
   static const G4double CI0[npi]={+1.00829205458740032E0,
              +.00845122624920943E0,+.00012700630777567E0,+.00007247591099959E0,
-						       +.00000513587726878E0,+.00000056816965808E0,+.00000008513091223E0,
-						       +.00000001238425364E0,+.00000000029801672E0,-.00000000078956698E0,
-						       -.00000000033127128E0,-.00000000004497339E0,+.00000000001799790E0,
-						       +.00000000000965748E0,+.00000000000038604E0,-.00000000000104039E0,
-						       -.00000000000023950E0,+.00000000000009554E0,+.00000000000004443E0,
-						       -.00000000000000859E0,-.00000000000000709E0,+.00000000000000087E0,
-							      +.00000000000000112E0,-.00000000000000012E0,-.00000000000000018E0};
+             +.00000513587726878E0,+.00000056816965808E0,+.00000008513091223E0,
+             +.00000001238425364E0,+.00000000029801672E0,-.00000000078956698E0,
+             -.00000000033127128E0,-.00000000004497339E0,+.00000000001799790E0,
+             +.00000000000965748E0,+.00000000000038604E0,-.00000000000104039E0,
+             -.00000000000023950E0,+.00000000000009554E0,+.00000000000004443E0,
+             -.00000000000000859E0,-.00000000000000709E0,+.00000000000000087E0,
+             +.00000000000000112E0,-.00000000000000012E0,-.00000000000000018E0};
 
   static const G4double CI1[npi]={+.975800602326285926E0,
@@ -215,9 +216,9 @@
 
   static const G4double EE[npk]={-.040172946544414076E0,-.444447147630558063E0,
-				      -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0,
-				      +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0,
-				      -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0,
-				      +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0,
-				      -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0};
+          -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0,
+          +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0,
+          -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0,
+          +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0,
+          -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0};
 
   static const G4complex CF[npj]={
@@ -244,5 +245,5 @@
   // -------------------------------------------------------------------------------------
   G4double H=0.;                     // Prototype of the result
-		if (ij)                            // I/K Bessel functions
+  if (ij)                            // I/K Bessel functions
   {
     if (ik)                          // I0/I1/EI0/EI1 Bessel functions (symmetric)
@@ -263,11 +264,11 @@
         G4double A1=1.+DY/(XLI*V3);
         G4double A2=1.+Y*(4.+(DY+Y)/(XLD*XL))/(W1*V3);
-				    G4double B0=1.;
-				    G4double	B1=1.-Y/XLI;
+        G4double B0=1.;
+        G4double B1=1.-Y/XLI;
         G4double B2=1.-Y*(1.-Y/(XLD+XLD))/W1;
         G4int    V1=3-XL;
         G4double V2=V3+V3;
         G4double C=0.;
-				    for (G4int N=3; N<=30; N++)
+        for (G4int N=3; N<=30; N++)
         {
           G4double C0=C;
@@ -279,10 +280,10 @@
           G4int    W5=W4-1;
           G4int    W6=W5-1;
-						             V1=V1+1;
+                   V1=V1+1;
                    V2=V2+1;
                    V3=V3+1;
-						    G4double U1=FN*W4;
+          G4double U1=FN*W4;
           G4double E=V3/(U1*W3);
-						    G4double U2=E*Y;
+          G4double U2=E*Y;
           G4double F1=1.+Y*V1/(U1*W1);
           G4double F2=(1.+Y*V2/(V3*W2*W5))*U2;
@@ -290,7 +291,7 @@
           G4double A=F1*A2+F2*A1+F3*A0;
           G4double B=F1*B2+F2*B1+F3*B0;
-						             C=A/B;
-						    if (std::abs(C0-C) < EPS*std::abs(C)) break;
-						    A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
+                   C=A/B;
+          if (std::abs(C0-C) < EPS*std::abs(C)) break;
+          A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
         }
         H=C;
@@ -301,22 +302,22 @@
       {
         G4double P=16./V-1.;
-				    G4double ALFA=P+P;
-				    G4double B1=0.;
+        G4double ALFA=P+P;
+        G4double B1=0.;
         G4double B2=0.;
-						  for (G4int I=npil; I>=0; I--)
+        for (G4int I=npil; I>=0; I--)
         {
           if (!nu) CJ=CI0[I];
           else     CJ=CI1[I];
-						    G4double B0=CJ+ALFA*B1-B2;
-						             B2=B1;
-						             B1=B0;
-        }
-				    H=std::sqrt(RPI2/V)*(B1-P*B2);
-				    if (nu && X < 0.) H=-H;
-				    if (!ex) H*=std::exp(V);
-				  }
+          G4double B0=CJ+ALFA*B1-B2;
+                   B2=B1;
+                   B1=B0;
+        }
+        H=std::sqrt(RPI2/V)*(B1-P*B2);
+        if (nu && X < 0.) H=-H;
+        if (!ex) H*=std::exp(V);
+      }
     }
     else                             // K0/K1/EK0/EK1 Bessel functions
-		  {
+    {
 #ifdef debug
       G4cout<<"G4BesIKJY: >>>>>>>>>>>>>> K is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
@@ -339,5 +340,5 @@
         G4double Q=HF;
         G4double C=1.;
-						  G4double D=B*B;
+        G4double D=B*B;
         G4double BK1=P;
         for (G4int N=1; N<=nat1; N++)  // @@ "nat" can be increased
@@ -349,8 +350,8 @@
                    C*=D/FN;
           G4double G=C*(P-FN*F);
-       			G4double R=C*F;
+          G4double R=C*F;
                    BK=BK+R;
                    BK1=BK1+G;
-								  if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break;
+          if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break;
         }
         if (nu==1) H=BK1/B;
@@ -358,5 +359,5 @@
         if (ex) H*=std::exp(X);
       }
-				  else if (X < 5.)
+      else if (X < 5.)
       {
 #ifdef debug
@@ -368,8 +369,8 @@
         G4double XN=DNUS+DNUS;
         G4double A=9.-XN;
-						  G4double B=25.-XN;
+        G4double B=25.-XN;
         G4double C=768*X*X;
-						  G4double HX=16*X;
-						  G4double C0=HX+HX+HX;;
+        G4double HX=16*X;
+        G4double C0=HX+HX+HX;;
         G4double A0=1.;
         G4double A1=(HX+7.+XN)/A;
@@ -378,6 +379,6 @@
         G4double B1=(HX+9.-XN)/A;
         G4double B2=(C+C0*B)/(A*B)+1.;
-						           C=0.;
-						  for (G4int N=3; N<=nat2; N++)
+                 C=0.;
+        for (G4int N=3; N<=nat2; N++)
         {
           C0=C;
@@ -389,18 +390,18 @@
           G4double FN3=FN1/(FN2-3.);
           G4double FN4=12*FN*FN-(1.-XN);
-								  G4double FN5=16*FN1*X;
+          G4double FN5=16*FN1*X;
           G4double RAN=1./(FNP*FNP-XN);
-								  G4double F1=FN3*(FN4-20*FN)+FN5;
+          G4double F1=FN3*(FN4-20*FN)+FN5;
           G4double F2=28*FN-FN4-8.+FN5;
           G4double F3=FN3*(FNM*FNM-XN);
-      				  							A=(F1*A2+F2*A1+F3*A0)*RAN;
+                   A=(F1*A2+F2*A1+F3*A0)*RAN;
                    B=(F1*B2+F2*B1+F3*B0)*RAN;
                    C=A/B;
-						    if (std::abs(C0-C) < EPS*std::abs(C)) break;
-								  A0=A1;	A1=A2;	A2=A;	B0=B1;	B1=B2;	B2=B;
-        }
-						  H=C/std::sqrt(RPIH*X);
+          if (std::abs(C0-C) < EPS*std::abs(C)) break;
+          A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
+        }
+        H=C/std::sqrt(RPIH*X);
         if (!ex) H*=std::exp(-X);
-				  }
+      }
       else
       {
@@ -410,10 +411,10 @@
         G4double P=10./X-1.;
         G4double ALFA=P+P;
-							 G4double B1=0.;
-							 G4double B2=0.;
+        G4double B1=0.;
+        G4double B2=0.;
 #ifdef debug
         G4cout<<"G4BesIKJY: >>> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
 #endif
-						  for (G4int I=npkl; I>=0; I--)
+        for (G4int I=npkl; I>=0; I--)
         {
           if (!nu) CK=CK0[I];
@@ -422,12 +423,12 @@
                    B2=B1;
                    B1=B0;
-						  }
+        }
         H=std::sqrt(PIH/X)*(B1-P*B2);
         if (!ex) H*=std::exp(-X);
       }
-		  }
+    }
   }
   else
-		{
+  {
     if (!ik && X < 0.)
     {
@@ -455,5 +456,5 @@
           B1=0.;
           B2=0.;
-										for (G4int JT=npkl; JT>=0; JT--)
+          for (G4int JT=npkl; JT>=0; JT--)
           {
             G4double B0=CB[JT]+ALFA*B1-B2;
@@ -463,21 +464,21 @@
           H=RPIH*(CE+std::log(HF*X))*H+B1-P*B2;
         }
-						}
+      }
       else
-						{
-								G4double P=10./V-1.;
-								G4double ALFA=P+P;
-								G4complex CB1(0.,0.);
-								G4complex CB2(0.,0.);
-							 for (G4int IT=npjl; IT>=0; IT--)
+      {
+        G4double P=10./V-1.;
+        G4double ALFA=P+P;
+        G4complex CB1(0.,0.);
+        G4complex CB2(0.,0.);
+        for (G4int IT=npjl; IT>=0; IT--)
         {
           G4complex CB0=CC[IT]+ALFA*CB1-CB2;
                     CB2=CB1;
                     CB1=CB0;
-								}
+        }
         CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI4))*(CB1-P*CB2);
-								if (ik) H=real(CB1);
+        if (ik) H=real(CB1);
         else    H=real(-CI*CB1);
-						}
+      }
     }
     else                          // J1/Y1 Bessel functions
@@ -494,5 +495,5 @@
         {
           G4double B0=CD[IT]+ALFA*B1-B2;
-										         B2=B1;
+                   B2=B1;
                    B1=B0;
         }
@@ -510,7 +511,7 @@
           H=RPIH*((CE+std::log(HF*X))*H-1./X)+Y*(B1-B2);
         }
-						}
+      }
       else
-						{
+      {
         G4double P=10./V-1.;
         G4double ALFA=P+P;
@@ -520,5 +521,5 @@
         {
           G4complex CB0=CF[IT]+ALFA*CB1-CB2;
-										          CB2=CB1;
+                    CB2=CB1;
                     CB1=CB0;
         }
@@ -526,5 +527,5 @@
         if (ik) H=real(CB1);
         else    H=real(-CI*CB1);
-						}
+      }
       if (X < 0.) H=-H;
     }
