1 #include "v3p_f2c.h"
2 
3 #ifdef KR_headers
4 extern double sqrt(), f__cabs();
5 
c_sqrt(r,z)6 VOID c_sqrt(r, z) complex *r, *z;
7 #else
8 #undef abs
9 #include "math.h"
10 #ifdef __cplusplus
11 extern "C" {
12 #endif
13 extern double f__cabs(double, double);
14 
15 #undef complex
16 #define complex v3p_netlib_complex
17 
18 void c_sqrt(complex *r, complex *z)
19 #endif
20 {
21         double mag, t;
22         double zi = z->i, zr = z->r;
23 
24         if( (mag = f__cabs(zr, zi)) == 0.)
25                 r->r = r->i = 0.;
26         else if(zr > 0)
27                 {
28                 r->r = t = sqrt(0.5 * (mag + zr) );
29                 t = zi / t;
30                 r->i = 0.5 * t;
31                 }
32         else
33                 {
34                 t = sqrt(0.5 * (mag - zr) );
35                 if(zi < 0)
36                         t = -t;
37                 r->i = t;
38                 t = zi / t;
39                 r->r = 0.5 * t;
40                 }
41         }
42 #ifdef __cplusplus
43 }
44 #endif
45