 /*
  * Khoros: $Id: complex.c,v 1.3 1992/03/20 23:42:22 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: complex.c,v 1.3 1992/03/20 23:42:22 dkhoros Exp $";
#endif

 /*
  * $Log: complex.c,v $
 * Revision 1.3  1992/03/20  23:42:22  dkhoros
 * VirtualPatch5
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1990, University of New Mexico.  All rights reserved.
 * 
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *----------------------------------------------------------------------
 */

#include "vmath.h"
#include "unmcopyright.h"	 /* Copyright 1990 by UNM */

cadd (r3,i3,r1,i1,r2,i2)
float r1,i1,r2,i2;
float *r3,*i3;
{
    /* Find complex result from (r1,i1)+(r2,i2) and put the result
       in (r3,i3) */

    *r3 = r1+r2;
    *i3 = i1+i2;
}

dcadd (r3,i3,r1,i1,r2,i2)
double r1,i1,r2,i2;
double *r3,*i3;
{
    /* Find double complex result from (r1,i1)+(r2,i2) and put the result
       in (r3,i3) */

    *r3 = r1+r2;
    *i3 = i1+i2;
}

csub (r3,i3,r1,i1,r2,i2)
float r1,i1,r2,i2;
float *r3,*i3;
{
    /* Find complex result from (r1,i1)-(r2,i2) and put the result
       in (r3,i3) */

    *r3 = r1-r2;
    *i3 = i1-i2;
}

dcsub (r3,i3,r1,i1,r2,i2)
double r1,i1,r2,i2;
double *r3,*i3;
{
    /* Find double complex result from (r1,i1)-(r2,i2) and put the result
       in (r3,i3) */

    *r3 = r1-r2;
    *i3 = i1-i2;
}

cmul(r3,i3,r1,i1,r2,i2)
float r1,i1,r2,i2;
float *r3,*i3;
  {
    /* Find complex result from (r1,i1)*(r2,i2) and put the result in
       (r3,i3) */

    *r3 = r1*r2-i1*i2;
    *i3 = r1*i2+i1*r2;
  }

dcmul(r3,i3,r1,i1,r2,i2)
double r1,i1,r2,i2;
double *r3,*i3;
  {
    /* Find double complex result from (r1,i1)*(r2,i2) and put the result in
       (r3,i3) */

    *r3 = r1*r2-i1*i2;
    *i3 = r1*i2+i1*r2;
  }

cdiv(r3,i3,r1,i1,r2,i2)
float r1,i1,r2,i2;
float *r3,*i3;
  {
    /* Find complex quotient from (r1,i1)/(r2,i2) and put the result it
       (r3,i3) */
    float f;

    f = r2*r2+i2*i2;
    *r3 = (r1*r2+i1*i2)/f;
    *i3 = (i1*r2-r1*i2)/f;
  }

dcdiv(r3,i3,r1,i1,r2,i2)
double r1,i1,r2,i2;
double *r3,*i3;
  {
    /* Find double complex quotient from (r1,i1)/(r2,i2) and put the result it
       (r3,i3) */
    float f;

    f = r2*r2+i2*i2;
    *r3 = (r1*r2+i1*i2)/f;
    *i3 = (i1*r2-r1*i2)/f;
  }

double vcabs(r,i)
float r,i;
  {
    /* Find magnitude of complex number */
    double mag;

    mag = sqrt((double)(r*r + i*i));
    return(mag);
   }

double dvcabs(r,i)
double r,i;
  {
    /* Find magnitude of double precision complex number */
    double mag;

    mag = sqrt(r*r + i*i);
    return(mag);
   }

double cmagsq(r,i)
float r,i;
  {
    /* Find power of complex number */
    double mag;

    mag = (double)(r*r + i*i);
    return(mag);
   }

double dcmagsq(r,i)
double r,i;
  {
    /* Find power of double precision complex number */
    double mag;

    mag = (double)(r*r + i*i);
    return(mag);
   }


double cang(r,i)
float r,i;
  {
    /* Find phasor angle of a complex number */
    double theta;
#ifndef atan2            /* check for inline macros */
    double atan2();
#endif

    /* Check for (0,0) case since the machines are giving back
       NaN's and Inf's! */
    if (r == 0.0 && i == 0.0) return(0.0);

    theta = atan2((double)i,(double)r);
    return(theta);
  }

double dcang(r,i)
double r,i;
  {
    /* Find phasor angle of a double precision complex number */
    double theta;
#ifndef atan2            /* check for inline macros */
    double atan2();
#endif

    /* Check for (0,0) case since the machines are giving back
       NaN's and Inf's! */
    if (r == 0.0 && i == 0.0) return(0.0);

    theta = atan2(i,r);
    return(theta);
  }

/*********************************************************************
*
*	PROGRAM: cquadratic.c
*
*       PURPOSE: This routine evaluates the following complex
*		 expression:
*			
*			s*s + coeff1*s + coeff2,
*
*		 where   s = real + imag. 
*
*        INPUT : 
*
*       OUTPUT : 
*
*   CALLED FROM: 
*
*   ROUTINES CALLED: 
*
*********************************************************************/


cquad ( real, imag, coeff1, coeff2 )
float coeff1, coeff2, *real, *imag;
{
	float tempr, tempi;

	cmul (&tempr, &tempi, *real, *imag, *real, *imag);
	cmul (real, imag, *real, *imag, coeff1, 0.0);
        cadd (real, imag, *real, *imag, coeff2, 0.0);
        cadd (real, imag, *real, *imag, tempr, tempi);
}

/*********************************************************************
*
*       PROGRAM: cconj.c
*
*       PURPOSE: this routine determines the conjugate of a complex
*                number 
*
*        INPUT : (r1,i1) - input complex pair
*                (*r2,*i2) - addresses of each member of output pair
*
*       OUTPUT : 
*
*   CALLED FROM:
*
*   ROUTINES CALLED:
*
*********************************************************************/
cconj(r2,i2,r1,i1)
float r1,i1;
float *r2,*i2;
{
  *r2 = r1;
  *i2 = (- i1);
}

