1#include "f2c.h"
2
3#ifdef KR_headers
4double sqrt(), f__cabs();
5VOID z_sqrt(r, z) doublecomplex *r, *z;
6#else
7#undef abs
8#include "math.h"
9extern double f__cabs(double, double);
10void z_sqrt(doublecomplex *r, doublecomplex *z)
11#endif
12{
13	double mag, zi = z->i, zr = z->r;
14
15	if( (mag = f__cabs(zr, zi)) == 0.)
16		r->r = r->i = 0.;
17	else if(zr > 0)
18		{
19		r->r = sqrt(0.5 * (mag + zr) );
20		r->i = zi / r->r / 2;
21		}
22	else
23		{
24		r->i = sqrt(0.5 * (mag - zr) );
25		if(zi < 0)
26			r->i = - r->i;
27		r->r = zi / r->i / 2;
28		}
29	}
30