X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=HBTAN%2Fvolya_complex.h;h=b96e69cebdda6b45f75128395b92f7a167c98b6a;hb=a9e18d1cddb941e12ceed61b5ecf29a61b8c3a81;hp=76603a716c48be9529526004f1ba3d6709405dbb;hpb=dd82cadcabd7d7498b1652c0136ecaf406c3c50c;p=u%2Fmrichter%2FAliRoot.git

diff --git a/HBTAN/volya_complex.h b/HBTAN/volya_complex.h
index 76603a716c4..b96e69cebdd 100644
--- a/HBTAN/volya_complex.h
+++ b/HBTAN/volya_complex.h
@@ -22,19 +22,21 @@ Version 1.0.1
 #define PI 3.141592653589793238462643
 #endif
 
-#include <iostream.h>
-#include <math.h>
+
+#include <TMath.h>
+#include <Riostream.h>
+//#include <math.h>
 #include <stdlib.h>
 
 /////////////////////////// Helper Error Function
 
-void ComplexError (char *msg)
+void ComplexError (const char *msg)
 {
 	cout << endl <<  msg << endl;
 	exit (1);
 }
 
-void ComplexWarning (char *msg)
+void ComplexWarning (const char *msg)
 {
 	cout << endl <<  msg << endl;
 }
@@ -165,7 +167,6 @@ public:
 
 ////////////////////////////////////////////////
 
-inline double sqr(double a)  {	return a*a; }
 
 inline Complex Complex::operator + () const
 {
@@ -201,7 +202,7 @@ inline Complex Complex::operator ! ()       //  complex conjugate  //
 
 double operator += (double &a, const Complex &b)
 {
-	if (fabs(b.Im) < small_epsilon) a+=b.Re;
+	if (TMath::Abs(b.Im) < small_epsilon) a+=b.Re;
 	else
 		ComplexError ("Error in double+=Complex: Complex is not double number");
 	return a;
@@ -209,7 +210,7 @@ double operator += (double &a, const Complex &b)
 
 double operator -= (double &a, const Complex &b)
 {
-	if (fabs(b.Im) < small_epsilon) a-=b.Re;
+	if (TMath::Abs(b.Im) < small_epsilon) a-=b.Re;
 	else
 		ComplexError ("Error in double -=Complex: Complex is not double number");
 	return a;
@@ -217,7 +218,7 @@ double operator -= (double &a, const Complex &b)
 
 double operator *= (double &a, const Complex &b)
 {
-	if (fabs(b.Im) < small_epsilon) a*=b.Re;
+	if (TMath::Abs(b.Im) < small_epsilon) a*=b.Re;
 	else
 		ComplexError ("Error in double *=Complex: Complex is not double number");
 	return a;
@@ -225,7 +226,7 @@ double operator *= (double &a, const Complex &b)
 
 double operator /= (double &a, const Complex &b)
 {
-	if (fabs(b.Im) < small_epsilon) a/=b.Re;
+	if (TMath::Abs(b.Im) < small_epsilon) a/=b.Re;
 	else
 		ComplexError ("Error in double /=Complex: Complex is not double number");
 	return a;
@@ -275,7 +276,7 @@ inline Complex Complex::operator *= (double a)
 Complex Complex::operator /= (const Complex &a)
 {
     double t1, t2, temp;
-    if (fabs(a.Re) >= fabs(a.Im))
+    if (TMath::Abs(a.Re) >= TMath::Abs(a.Im))
     {
       t1   = a.Im / a.Re;
       t2   = a.Re + a.Im * t1;
@@ -358,7 +359,7 @@ inline Complex operator / (const Complex &a, double b)
 inline Complex operator / (const Complex &a, const Complex &b)
 {
     double t1, t2;
-    if (fabs(b.Re) >= fabs(b.Im))
+    if (TMath::Abs(b.Re) >= TMath::Abs(b.Im))
     {
       t1   = b.Im / b.Re;
       t2   = b.Re + b.Im * t1;
@@ -395,12 +396,12 @@ istream& operator >> (istream &stream, Complex &a)
 
 inline int operator == (const Complex &a, double b)
 {
-	return (a.Re==b) && (fabs(a.Im)<small_epsilon);
+	return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
 }
 
 inline int operator == (double b, const Complex &a)
 {
-	return (a.Re==b) && (fabs(a.Im)<small_epsilon);
+	return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
 }
 
 inline int operator == (const Complex &a, const Complex &b)
@@ -410,12 +411,12 @@ inline int operator == (const Complex &a, const Complex &b)
 
 inline int operator != (const Complex &a, double b)
 {
-	return (a.Re!=b) || (fabs(a.Im)>small_epsilon);
+	return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
 }
 
 inline int operator != (double b, const Complex &a)
 {
-	return (a.Re!=b) || (fabs(a.Im)>small_epsilon);
+	return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
 }
 
 inline int operator != (const Complex &a, const Complex &b)
@@ -498,12 +499,12 @@ double  imag (const Complex &a)
 }
 double abs (const Complex &a)
 {
-	if (a.Im == 0) return fabs(a.Re);
-	if (a.Re == 0) return fabs(a.Im);
-    double R = fabs(a.Re), I = fabs(a.Im);
+	if (a.Im == 0) return TMath::Abs(a.Re);
+	if (a.Re == 0) return TMath::Abs(a.Im);
+    double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im);
 	return (R >= I) ?
-           (R * sqrt (1 + sqr(a.Im/a.Re))) :
-           (I * sqrt (1 + sqr(a.Re/a.Im))) ;
+           (R * sqrt (1 + a.Im*a.Im/a.Re/a.Re)) :
+           (I * sqrt (1 + a.Re*a.Re/a.Im/a.Im)) ;
 }
 
 double  Arg (const Complex &a)
@@ -527,12 +528,12 @@ inline Complex sqr (const Complex &a) { return a*a; }
 
 Complex sqrt (const Complex &a, int flag)
 {
-	if ((a.Re>=0) && (fabs(a.Im)<small_epsilon))
+	if ((a.Re>=0) && (TMath::Abs(a.Im)<small_epsilon))
        return flag ? -sqrt(a.Re) : sqrt(a.Re);
-    double R = fabs(a.Re), I = fabs(a.Im);
+    double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im);
 	double w = (R >= I) ?
-           sqrt (R/2 * (  1 + sqrt (1 + sqr(a.Im/a.Re)))):
-           sqrt (I/2 * (R/I + sqrt (1 + sqr(a.Re/a.Im))));
+           sqrt (R/2 * (  1 + sqrt (1 + a.Im*a.Im/a.Re/a.Re))):
+           sqrt (I/2 * (R/I + sqrt (1 + a.Re*a.Re/a.Im/a.Im)));
     Complex c;
     if (a.Re >= 0)
     {
@@ -760,8 +761,8 @@ void Solve3 (Complex* z, const Complex &a2, const Complex &a1, const Complex &a0
   t = sqrt((r^2) - (q^3));
 	a = ((!r * t) >=0.0) ? -((r+t)^(1.0/3)) : -((r-t)^(1.0/3));
   b = ((a == zero) ? zero : (q/a));
-  z[0] = -(a+b)/2 - a2/3 + ImUnit*sqrt(3)*(a-b)/2;
-  z[1] = -(a+b)/2 - a2/3 - ImUnit*sqrt(3)*(a-b)/2;
+  z[0] = -(a+b)/2 - a2/3 + ImUnit*sqrt(3.0)*(a-b)/2;
+  z[1] = -(a+b)/2 - a2/3 - ImUnit*sqrt(3.0)*(a-b)/2;
   z[2] = a + b - a2/3;
 }