/* complex.c   implement complex arithmetic operations                 */
/*             also elementary functions and vector, matrix operations */

#include "complex.h"
#include <stdio.h>
#include <math.h>
#include "udrnrt.h"
#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))

complex cmplx(double  a, double  b)  /* make a complex number */
{
  complex c;
  c.re = a;
  c.im = b;
  return c;
}

double  real(complex a) /* get real part */
{
  return a.re;
}

double  imag(complex a) /* get imaginary part */
{
  return a.im;
}

complex cxadd(complex a, complex b)  /* add two complex numbers */
{
  complex c;
  c.re = a.re+b.re;
  c.im = a.im+b.im;
  return c;
}

complex cxsub(complex a, complex b)  /* subtract b from a, complex */
{
  complex c;
  c.re = a.re-b.re;
  c.im = a.im-b.im;
  return c;
}

complex cxmul(complex a, complex b)  /* multiply two complex numbers */
{
  complex c;
  c.re = a.re*b.re - a.im*b.im;
  c.im = a.re*b.im + a.im*b.re;
  return c;
}

complex cdmul(double   a, complex b)  /* multiply double times complex */
{
  complex c;
  c.re = a*b.re;
  c.im = a*b.im;
  return c;
}

complex cmuld(complex a, double   b)  /* multiply complex times double */
{
  complex c;
  c.re = a.re*b;
  c.im = a.im*b;
  return c;
}

complex cxdiv(complex a, complex b)  /* divide complex by complex */
{
  complex c;
  double  r;
  r = b.re*b.re+b.im*b.im;
  c.re = (a.re*b.re+a.im*b.im)/r;
  c.im = (a.im*b.re-a.re*b.im)/r;
  return c;
}

complex cdivd(complex a, double   b)  /* divide complex by double */
{
  complex c;
  c.re = a.re/b;
  c.im = a.im/b;
  return c;
}

double   cxabs(complex a)             /* return magnitude of complex munber */
{
  return sqrt(a.re*a.re+a.im*a.im);
}

double   cxarg(complex a)             /* return argument of complex munber */
{
  return atan2(a.im,a.re);
}

complex cxneg(complex a) /* negate a complex */
{
  complex c;
  c.re = -a.re;
  c.im = -a.im;
  return c;
}

complex cxconj(complex a) /* conjugate a complex */
{
  complex c;
  c.re = a.re;
  c.im = -a.im;
  return c;
}

double sumabs(complex a) /* abs re + abs im */
{
  return abs(a.re)+abs(a.im);
}

/* complex elementary functions */

complex cxsqrt(complex a) /* square root of a complex number*/
{
  complex c;
  double r, ra, sqre, sqim;
  
  r = cxabs(a);
  sqre = sqrt((r+a.re)/2.0);
  sqim = sqrt((r-a.re)/2.0);
  if(a.im<0.0) sqim = -sqim;
  return cmplx(sqre, sqim);
}

complex cxexp(complex a) /* e^a, exp of a complex number */
{
  complex v;
  double erea;

  erea = exp(a.re);
  v = cmplx(erea*cos(a.im),erea*sin(a.im));
  return v;
}

complex cxlog(complex a) /* log base e of a complex number */
{
  complex v;
  double vmod, varg;

  vmod = cxabs(a);
  varg = cxarg(a);
  v = cmplx(log(vmod), varg);
  return v;
}

complex cxsin(complex a) /* sine of a complex number */
{
  complex v;

  v = cmplx(sin(a.re)*cosh(a.im), cos(a.re)*sinh(a.im));
  return v;
}

complex cxcos(complex a) /* cosine of a complex number */
{
  complex v;

  v = cmplx(cos(a.re)*cosh(a.im), sin(a.re)*sinh(a.im));
  return v;
}

complex cxtan(complex a) /* tangent of a complex number */
{
  complex v;

  v = cxdiv(cxsin(a),cxcos(a));
  return v;
}

complex cxatan(complex a) /* arctangent of a complex number */
{
  complex t, v;

  t = cxatanh(a);
  v = cmplx(t.im, -t.re);
  return v;
}

complex cxsinh(complex a) /* hyperbolic sine of a complex number */
{
  complex v;

  v = cmplx(sinh(a.re)*cos(a.im), cosh(a.re)*sin(a.im));
  return v;
}

complex cxcosh(complex a) /* hyperbolic cosine of a complex number */
{
  complex v;

  v = cmplx(cosh(a.re)*cos(a.im), sinh(a.re)*sin(a.im));
  return v;
}

complex cxtanh(complex a) /* hyperbolic tangent of a complex number */
{
  complex v;

  v = cxdiv(cxsinh(a),cxcosh(a));
  return v;
}

complex cxatanh(complex a) /* hyperbolic arctangent of a complex number */
{
  complex one, v;
  complex lap1, lam1;

  one = cmplx(1.0,0.0);
  lap1 = cxlog(cxadd(one,a));
  lam1 = cxlog(cxsub(one,a));
  v = cdivd(cxsub(lap1,lam1),2.0);
  return v;
}

/* complex vecotor, matrix operations */

void cxvadd(int n, complex A[n], complex B[n], complex C[n]) /* C = A + B */
{
  int i;

  for(i=0; i<n; i++)
    C[i] = cxadd(A[i],B[i]);
}

void cxvsub(int n, complex A[n], complex B[n], complex C[n])  /* C = A - B */
{
  int i;

  for(i=0; i<n; i++)
    C[i] = cxsub(A[i],B[i]);
}

void cxvcopy(int n, complex A[n], complex B[n])  /* copy vector A to B */
{
  int i;

  for(i=0; i<n; i++)
    B[i] = A[i];
}

double cxvnorm(int n, complex A[n]) /* vector norm of A */
{
  int i;
  double sum = 0.0;
  
  for(i=0; i<n; i++)
    sum = sum + cxabs(A[i]);
  return sqrt(sum);
}

void cxvident(int n, int i, complex A[n])  /* identity vector at i */
{
  int j;
  complex zero;

  zero = cmplx(0.0,0.0);
  for(j=0; j<n; j++)
    A[j] = zero;
  A[i] =  cmplx(1.0,0.0);
}

void cxvprint(int n, complex A[n]) /* print complex vector */
{
  int i;

  for(i=0; i<n; i++)
    printf("A[%d]=(%g,%g) \n", i, A[i].re, A[i].im);

}

void cxvran(int n, complex A[n]) /* generate random complex vector */
{
  int i;

  for(i=0; i<n; i++)
    A[i] = cmplx(udrnrt(),udrnrt());

}

void cxmadd(int n, complex A[n][n], complex B[n][n], complex C[n][n]) /*C=A+B*/
{
  int i, j;

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      C[i][j] = cxadd(A[i][j],B[i][j]);
}

void cxmsub(int n, complex A[n][n], complex B[n][n], complex C[n][n]) /*C=A-B*/
{
  int i, j;

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      C[i][j] = cxsub(A[i][j],B[i][j]);
}

void cxmmul(int n, complex A[n][n], complex B[n][n], complex C[n][n]) /*C=A*B*/
{
  int i, j, k;
  complex zero = cmplx(0.0,0.0);

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
    {
      C[i][j] = zero;
      for(k=0; k<n; k++)
        C[i][j] = cxadd(C[i][j],cxmul(A[i][k],B[k][j]));
    }
}

void cxmcopy(int n, complex A[n][n], complex B[n][n]) /* copy matrix A to B */
{
  int i, j;

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      B[i][j] = A[i][j];
}

double cxmnorm(int n, complex A[n][n]) /* matrix norm of A */
{
  int i, j;
  double sum = 0.0;
  
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      sum = sum + cxabs(A[i][j]);
  return sqrt(sum);
}

void cxmident(int n, complex A[n][n])      /* initialize identity matrix   */
{
  int i, j;
  complex zero, one;

  zero = cmplx(0.0,0.0);
  one = cmplx(1.0,0.0);
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
      A[i][j] = zero;
    A[i][i] = one;
  }
}

void cxmprint(int n, complex A[n][n])      /* print complex matrix         */
{
  int i, j;

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      printf("A[%d][%d]=(%g,%g) \n", i, j, A[i][j].re, A[i][j].im);
}

void cxmran(int n, complex A[n][n]) /* generate random complex matrix */
{
  int i, j;

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      A[i][j] = cmplx(udrnrt(),udrnrt());

}

void cxmvmul(int n, complex A[n][n], complex B[n], complex C[n])
{
  int i, j, k;
  complex zero = cmplx(0.0,0.0);

  for(i=0; i<n; i++)
  {
    C[i] = zero;
    for(j=0; j<n; j++)
      C[i] = cxadd(C[i],cxmul(A[i][j],B[j]));
  }
}

void cxsimeq(int n, complex A[n][n], complex Y[n], complex X[n])
{
  /* solve linear system of equations with complex coefficients  */
  /*       [A] * |X| = |Y|, |A| destroyed                        */
  int row[n];          /* row interchange indices */
  int hold , I_pivot;  /* pivot indices */
  complex pivot;       /* pivot element value */
  double abs_pivot;
  int i,j,k;
  complex temp[n];
  complex zero = cmplx(0.0,0.0);

  cxvcopy(n, Y, X); /* RHS becomes solution */
  /* set up row interchange vector */
  for(k=0; k<n; k++) row[k] = k;
  for(k=0; k<n; k++)
  {
    /* find largest element for povot */
    pivot = A[row[k]][k];
    abs_pivot = cxabs(pivot);
    I_pivot = k;
    for(i=k; i<n; i++)
    {
      if(cxabs(A[row[i]][k]) > abs_pivot)
      {
        I_pivot = i;
        pivot = A[row[i]][k];
        abs_pivot = cxabs(pivot);
      }
    }
    hold = row[k];
    row[k] = row[I_pivot];
    row[I_pivot] = hold;
    if(abs_pivot < 1.0E-15 ) /* singular, delete row */
    {
      for(j=k+1; j<n; j++)
        A[row[k]][j] = zero;
      X[row[k]] = zero;
      printf("redundant row (singular) %d \n", row[k]);
    }
    else
    {
      /* reduce about pivot */
      for(j=k+1; j<n; j++)
        A[row[k]][j] = cxdiv(A[row[k]][j],A[row[k]][k]);
      X[row[k]] = cxdiv(X[row[k]],A[row[k]][k]);
      /* inner reduction loop */
      for(i=0; i<n; i++)
      {
        if( i != k)
	{
          for(j=k+1; j<n; j++)
            A[row[i]][j] = cxsub(A[row[i]][j],cxmul(A[row[i]][k],A[row[k]][j]));
          X[row[i]] = cxsub(X[row[i]],cxmul(A[row[i]][k],X[row[k]]));
	}
      }
    }
  }
  /* unscramble rows */
  for(i=0; i<n; i++)
    temp[i] = X[i];
  for(i=0; i<n; i++)
    X[i] = temp[row[i]];
} /* end cxsimeq */

complex cxdet(int n, complex A[n][n])
{
  complex D = cmplx(1.0,0.0);  /* determinant */
  complex B[n][n];  /* working matrix */
  int row[n];  /* row interchange indicies */
  int hold , I_pivot;  /* pivot indicies */
  complex pivot;  /* pivot element value */
  double abs_pivot;
  int i, j, k;
  complex zero = cmplx(0.0,0.0);

  cxmcopy(n,A,B);
  for(k=0; k<n; k++) row[k]= k;
  for(k=0; k<n-1; k++)
  {
    pivot = B[row[k]][k];
    abs_pivot = cxabs(pivot);
    I_pivot = k;
    for(i=k; i<n; i++){
      if( cxabs(B[row[i]][k]) > abs_pivot )
      {
        I_pivot = i;
        pivot = B[row[i]][k];
        abs_pivot = cxabs(pivot);
      }
    }
    if(I_pivot != k)
    {
      hold = row[k];
      row[k] = row[I_pivot];
      row[I_pivot] = hold;
      D = cxneg(D);
    }
    if(abs_pivot < 1.0E-15) return zero;
    D = cxmul(D,B[row[k]][k]);
    for(j=k+1; j<n; j++)
      B[row[k]][j] = cxdiv(B[row[k]][j],B[row[k]][k]);
    for(i=0; i<n; i++)
      if(i != k)
        for(j=k+1; j<n; j++)
          B[row[i]][j] = cxsub(B[row[i]][j],cxmul(B[row[i]][k],B[row[k]][j]));
  }
  return cxmul(D,B[row[n-1]][n-1]);
} /* end cxdet */

void cxinv(int n, complex A[n][n], complex AA[n][n])
{
    int row[n], col[n];
    complex temp[n];
    int hold , I_pivot , J_pivot;
    double abs_pivot;
    complex pivot;
    int i, j, k;
    complex one = cmplx(1.0,0.0);

    cxmcopy(n, A, AA);
    for (k=0; k<n; k++)
    {
      row[k] = k;
      col[k] = k;
    }
    for (k=0; k<n; k++)
    {
      pivot = AA[row[k]][col[k]] ;
      I_pivot = k;
      J_pivot = k;
      for(i=k; i<n; i++)
      {
        for(j=k; j<n; j++)
        {
          abs_pivot = cxabs(pivot) ;
          if(cxabs(AA[row[i]][col[j]]) > abs_pivot )
          {
            I_pivot = i ;
            J_pivot = j ;
            pivot = AA[row[i]][col[j]] ;
          }
        }
      }
      if(cxabs(pivot) < 1.0E-15)
      {
        printf("MATRIX is SINGULAR !!! \n");
        return;
      }
      hold = row[k];
      row[k]= row[I_pivot];
      row[I_pivot] = hold ;
      hold = col[k];
      col[k]= col[J_pivot];
      col[J_pivot] = hold ;
      AA[row[k]][col[k]] = cxdiv(one,pivot) ;
      for (j=0; j<n; j++)
        if( j != k )
          AA[row[k]][col[j]] = cxmul(AA[row[k]][col[j]],AA[row[k]][col[k]]);
      for (i=0; i<n; i++)
      {
        if( k != i )
        {
          for (j=0; j<n; j++)
          {
            if( k != j )
            {
              AA[row[i]][col[j]] = cxsub(AA[row[i]][col[j]], 
				 cxmul(AA[row[i]][col[k]],AA[row[k]][col[j]]));
            }
          }
          AA[row[i]][col[k]] = cxneg(
                                cxmul(AA[row[i]][col[k]],AA[row[k]][col[k]]));
        }
      }
    }
    /* unscramble rows */
    for (j=0; j<n; j++)
    {
      for (i=0; i<n; i++)
        temp[col[i]] = AA[row[i]][j];
      for (i=0; i<n; i++)
        AA[i][j] = temp[i] ;
    }
    /* unscramble columns */
    for (i=0; i<n; i++)
    {
      for (j=0; j<n; j++)
        temp[row[j]] = AA[i][col[j]] ;
      for (j=0; j<n; j++)
        AA[i][j] = temp[j] ;
    }
} /* end cxinv */

/* end complex.c */