diff options
| author | goodale <goodale@17b73243-c579-4c4c-a9d2-2d5706c11dac> | 2005-01-24 23:56:11 +0000 |
|---|---|---|
| committer | goodale <goodale@17b73243-c579-4c4c-a9d2-2d5706c11dac> | 2005-01-24 23:56:11 +0000 |
| commit | aea70b578fc76529f1519ce85b1b4591db01d80d (patch) | |
| tree | 571569b509f8e3fa0e25d79aea05bfba8711af97 | |
| parent | 910923177db17fbb3f9970463ff22bc61a494939 (diff) | |
Patch from Yaakoub to add a complex version of the 'pow' function.
git-svn-id: http://svn.cactuscode.org/flesh/trunk@3966 17b73243-c579-4c4c-a9d2-2d5706c11dac
| -rw-r--r-- | src/include/cctk_Complex.h | 6 | ||||
| -rw-r--r-- | src/main/Complex.c | 79 |
2 files changed, 83 insertions, 2 deletions
diff --git a/src/include/cctk_Complex.h b/src/include/cctk_Complex.h index e47b7e59..041884f3 100644 --- a/src/include/cctk_Complex.h +++ b/src/include/cctk_Complex.h @@ -31,7 +31,8 @@ cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b); \ cctk_complex CCTK_Cmplx##Sin (cctk_complex complex_number); \ cctk_complex CCTK_Cmplx##Cos (cctk_complex complex_number); \ cctk_complex CCTK_Cmplx##Exp (cctk_complex complex_number); \ -cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number); +cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number); \ +cctk_complex CCTK_Cmplx##Pow (cctk_complex complex_number, cctk_real w); /* declare complex functions for all available precisions */ @@ -63,6 +64,7 @@ DECLARE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32) #define CCTK_CmplxCos CCTK_Cmplx8Cos #define CCTK_CmplxExp CCTK_Cmplx8Exp #define CCTK_CmplxSqrt CCTK_Cmplx8Sqrt +#define CCTK_CmplxPow CCTK_Cmplx8Pow #elif CCTK_REAL_PRECISION_8 #define CCTK_Cmplx CCTK_Cmplx16 #define CCTK_CmplxReal CCTK_Cmplx16Real @@ -77,6 +79,7 @@ DECLARE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32) #define CCTK_CmplxCos CCTK_Cmplx16Cos #define CCTK_CmplxExp CCTK_Cmplx16Exp #define CCTK_CmplxSqrt CCTK_Cmplx16Sqrt +#define CCTK_CmplxPow CCTK_Cmplx16Pow #elif CCTK_REAL_PRECISION_16 #define CCTK_Cmplx CCTK_Cmplx32 #define CCTK_CmplxReal CCTK_Cmplx32Real @@ -91,6 +94,7 @@ DECLARE_CMPLX_FUNCTIONS (CCTK_Cmplx32, CCTK_REAL16, CCTK_COMPLEX32) #define CCTK_CmplxCos CCTK_Cmplx32Cos #define CCTK_CmplxExp CCTK_Cmplx32Exp #define CCTK_CmplxSqrt CCTK_Cmplx32Sqrt +#define CCTK_CmplxPow CCTK_Cmplx32Pow #endif #ifdef __cplusplus diff --git a/src/main/Complex.c b/src/main/Complex.c index 998978a3..9415c2f2 100644 --- a/src/main/Complex.c +++ b/src/main/Complex.c @@ -507,6 +507,82 @@ cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number) \ return (result); \ } + /*@@ + @routine CCTK_CmplxPow + @date + @author Yaakoub Y El Khamra + @desc + Raises a complex number to a given power + @enddesc + + @var complex_number + @vdesc The complex number + @vtype CCTK_COMPLEX + @vio in + @endvar + + @returntype CCTK_COMPLEX + @returndesc + The square root + @endreturndesc +@@*/ +#define DEFINE_CCTK_CMPLX_POW(CCTK_Cmplx, cctk_real, cctk_complex) \ +cctk_complex CCTK_Cmplx##Pow (cctk_complex complex_number, cctk_real w) \ +{ \ + cctk_real R, theta, phi; \ + cctk_complex result; \ + \ + if ( complex_number.Re > 0) \ + { \ + theta = 0.0; \ + phi = atan(complex_number.Im/complex_number.Re) + theta; \ + R = sqrt(complex_number.Re*complex_number.Re + complex_number.Im*complex_number.Im);\ + R = pow(R,w); \ + result.Re = R * cos (w * phi); \ + result.Im = R * sin (w * phi); \ + } \ + else if ( complex_number.Re < 0 && complex_number.Im >= 0 ) \ + { \ + theta = 4.0 * atan (1.0); \ + phi = atan(complex_number.Im/complex_number.Re) + theta; \ + R = sqrt(complex_number.Re*complex_number.Re + complex_number.Im*complex_number.Im);\ + R = pow(R,w); \ + result.Re = R * cos (w * phi); \ + result.Im = R * sin (w * phi); \ + } \ + else if ( complex_number.Re < 0 && complex_number.Im < 0 ) \ + { \ + theta = 4.0 * atan (1.0); \ + phi = atan(complex_number.Im/complex_number.Re) + theta; \ + R = sqrt(complex_number.Re*complex_number.Re + complex_number.Im*complex_number.Im);\ + R = pow(R,w); \ + result.Re = R * cos (w * phi); \ + result.Im = R * sin (w * phi); \ + } \ + else if ( fabs(complex_number.Re) <= 1e-20 && fabs(complex_number.Im) < 1e-20 )\ + { \ + result.Re = 0.0; \ + result.Im = 0.0; \ + } \ + else if ( fabs(complex_number.Re) <= 1e-20 && complex_number.Im < 0 ) \ + { \ + phi = 2.0 * atan (1.0); \ + R = sqrt(complex_number.Re*complex_number.Re + complex_number.Im*complex_number.Im);\ + R = pow(R,w); \ + result.Re = R * cos (w * phi); \ + result.Im = R * sin (w * phi); \ + } \ + if ( fabs(complex_number.Re) <= 1e-20 && complex_number.Im < 0 ) \ + { \ + phi = -2.0 * atan (1.0); \ + R = sqrt(complex_number.Re*complex_number.Re + complex_number.Im*complex_number.Im);\ + R = pow(R,w); \ + result.Re = R * cos (w * phi); \ + result.Im = R * sin (w * phi); \ + } \ + \ + return (result); \ +} /* macro to define a set of complex functions for a given precision */ #define DEFINE_CMPLX_FUNCTIONS(CCTK_Cmplx, cctk_real, cctk_complex) \ @@ -522,7 +598,8 @@ cctk_complex CCTK_Cmplx##Sqrt (cctk_complex complex_number) \ DEFINE_CCTK_CMPLX_SIN (CCTK_Cmplx, cctk_real, cctk_complex) \ DEFINE_CCTK_CMPLX_COS (CCTK_Cmplx, cctk_real, cctk_complex) \ DEFINE_CCTK_CMPLX_EXP (CCTK_Cmplx, cctk_real, cctk_complex) \ - DEFINE_CCTK_CMPLX_SQRT (CCTK_Cmplx, cctk_real, cctk_complex) + DEFINE_CCTK_CMPLX_SQRT (CCTK_Cmplx, cctk_real, cctk_complex) \ + DEFINE_CCTK_CMPLX_POW (CCTK_Cmplx, cctk_real, cctk_complex) \ /* define complex functions for all available precisions */ #ifdef CCTK_REAL4 |