/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:12 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "zsqrtx.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	HALF	0.5e0
#define	ONE	1.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ zsqrtx(
double cin[],
double cout[])
{
	double a, a1, b, b1, r, x, y;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	double *const Cin = &cin[0] - 1;
	double *const Cout = &cout[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 2001-01-24 ZSQRTX Krogh  ZSQRT -> ZSQRTX to avoid C lib. conflicts.
	 *>> 1995-10-30 ZSQRT  Krogh  Fixed so M77CON can get S.P. for C conv.
	 *>> 1987-12-07 ZSQRT  Lawson Initial code.
	 *--Z replaces "?": ?SQRTX
	 *
	 *     Reference: Uspensky, "Theory of Equations", McGraw-Hill,
	 *     1948, Page 14-15.
	 *     Computes one of the square roots of the double precision
	 *     complex number whose real part is given in CIN(1) and
	 *     imaginary part is given in CIN(2). Result is given
	 *     similarly in COUT(1) and COUT(2).
	 *
	 *     C.L.Lawson & S.Y.Chan, JPL, June 3,1986.
	 *
	 *     ------------------------------------------------------------------
	 * */
 
 
	a = Cin[1];
	b = Cin[2];
	if (b == ZERO)
	{
		if (a == ZERO)
		{
			Cout[1] = ZERO;
			Cout[2] = ZERO;
		}
		else if (a > ZERO)
		{
			Cout[1] = sqrt( a );
			Cout[2] = ZERO;
		}
		else
		{
			Cout[1] = ZERO;
			Cout[2] = sqrt( -a );
		}
	}
	else
	{
		/*                                          Here B .ne. 0 */
		if (a == ZERO)
		{
			Cout[1] = sqrt( HALF*fabs( b ) );
			Cout[2] = Cout[1];
		}
		else
		{
			a1 = fabs( a );
			b1 = fabs( b );
			if (b1 > a1)
			{
				r = b1*sqrt( ONE + powi(a1/b1,2) );
			}
			else
			{
				r = a1*sqrt( ONE + powi(b1/a1,2) );
			}
			x = sqrt( HALF*(r + a1) );
			y = HALF*b1/x;
			if (a > ZERO)
			{
				Cout[1] = x;
				Cout[2] = y;
			}
			else
			{
				Cout[1] = y;
				Cout[2] = x;
			}
		}
		if (b < ZERO)
			Cout[2] = -Cout[2];
	}
 
	return;
} /* end of function */