/*******************************************
 * Copyright (C) 1990 by Christopher R. Wren and 
 * the Massachusetts Institute of Technology
 *
 * COMMERCIAL DISCLAIMER:
 *
 * This software is intended for research purposes only and contains
 * unlicensed software, in original and modified forms, from various
 * sources including the Free Software Foundation.  These programs
 * shall not be used, rewritten, or adapted as the basis of a
 * commercial software or hardware product without first obtaining
 * appropriate licenses from the owners of the original software
 * copyrights.
 *
 * In the case of copyrighted software fragments obtained from the
 * Free Software Foundation: you can redistribute it and/or modify it
 * under the terms of the GNU Library General Public License as
 * published by the Free Software Foundation.  This software is
 * distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General
 * Public License for more details.  You should have received a copy
 * of the GNU Library General Public License and the GNU General
 * Public License along with this software; if not, write to the Free
 * Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
 * 
 * Developed  by Christopher R. Wren at the 
 * Media Laboratory, MIT, Cambridge, Massachusetts.
 * 
 *******************************************/

#include <complex.h>

fcomplex Cadd(a,b)
fcomplex a,b;
{ 
  fcomplex c;
  c.r=a.r+b.r;
  c.i=a.i+b.i;
  return c;
}

fcomplex Csub(a,b)
fcomplex a,b;
{
  fcomplex c;
  c.r=a.r-b.r;
  c.i=a.i-b.i;
  return c;
}

fcomplex Cmul(a,b)
fcomplex a,b;
{
  fcomplex c;
  c.r=a.r*b.r-a.i*b.i;
  c.i=a.i*b.r+a.r*b.i;
  return c;
}

fcomplex Cdiv(a,b)
fcomplex a,b;
{
  fcomplex c;
  float r,den;
  if (fabs(b.r) >= fabs(b.i)) {
    r=b.i/b.r;
    den=b.r+r*b.i;
    c.r=(a.r+r*a.i)/den;
    c.i=(a.i-r*a.r)/den;
  } else {
    r=b.r/b.i;
    den=b.i+r*b.r;
    c.r=(a.r*r+a.i)/den;
    c.i=(a.i*r-a.r)/den;
  }
  return c;
}

fcomplex Complex(re,im)
float re,im;
{
  fcomplex c;
  c.r=re;
  c.i=im;
  return c;
}

float Cabs(z)
fcomplex z;
{
  float x,y,ans,temp;
  
  x=fabs(z.r);
  y=fabs(z.i);
  if (x == 0.0)
    ans=y;
  else if (y == 0.0)
    ans = x;
  else if (x>y) {
    temp=y/x;
    ans=x*sqrt(1.0+temp*temp);
  } else {
    temp=x/y;
    ans=y*sqrt(1.0+temp*temp);
  }
return ans;
}

fcomplex Conjg(z)
fcomplex z;
{
  fcomplex c;
  c.r=z.r;
  c.i= -z.i;
  return c;
}

fcomplex Csqrt(z)
fcomplex z;
{
  fcomplex c;
  float x,y,w,r;
  if ((z.r == 0.0) && (z.i == 0.0)) {
    c.r=c.i=0.0;
    return c;
  } else {
    x = fabs(z.r);
    y=fabs(z.i);
    if (x>= y) {
      r=y/x;
      w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
    } else {
      r=x/y;
      w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
    }
    if (z.r >= 0.0) {
      c.r=w;
      c.i=z.i/(2.0*w);
    } else {
      c.i=(z.i >= 0) ? w : -w;
      c.r=z.i/(2.0*c.i);
    }
    return c;
  }
}

fcomplex RCmul(x,a)
float x;
fcomplex a;
{
  fcomplex c;
  c.r=x*a.r;
  c.i=x*a.i;
  return c;
}

	     
fcomplex Cpow(z,i)
fcomplex z;
int i;
{
  fcomplex x;
  int j=1;
  x=z;
  while (i > j) {
    x=Cmul(x,z);
    j++;
  }
  return x;
}