/*
 *     27/Oct/1994
 *     programed by Akihiro.Sato
 * 
 *     Complex funcions
 *
 */

#include <stdio.h>
#include <math.h>
#include "complexmath.h"


void MakeComplex( real, imag, rz)  /* rz = real + i * imag */
double real,imag;
COMPLEX *rz;
{
  rz -> x = real;
  rz -> y = imag;
}  /* end. MakeComplex */


void ZeroSetofComplex( rz)  /* rz = 0.0 + i * 0.0 */
COMPLEX *rz;
{
  MakeComplex( 0.0, 0.0, rz);
}  /* end. ZeroSetofComplex */


void ESetofComplex( rz)  /* rz = 1.0 + i * 0.0 */
COMPLEX *rz;
{
  MakeComplex( 1.0, 0.0, rz);
}  /* end. ESetofComplex */


void ISetofComplex( rz)  /* rz = 0.0 + i * 1.0 */
COMPLEX *rz;
{
  MakeComplex( 0.0, 1.0, rz);
}  /* end. ISetofComplex */


void InfSetofComplex( rz)  /* rz = Infinity */
COMPLEX *rz;
{
  MakeComplex( INFINITY, 0.0, rz);
}  /* end. InfSetofComplex */

void NaNSetofComplex( rz)  /* rz = NaN + i * NaN */
COMPLEX *rz;
{
  MakeComplex( quiet_nan(0), quiet_nan(0), rz);
}  /* end. NaNSetofComplex */


void SetofComplex( z, rz)  /* rz = z */
COMPLEX z,*rz;
{
  rz -> x = z.x;
  rz -> y = z.y;
}  /* end. SetofComplex */


void ScalarTimesofComplex( c, z, rz)  /* rz = c * z */
double c;
COMPLEX z,*rz;
{
  rz -> x = c * z.x;
  rz -> y = c * z.y;
}  /* end. ScalarTimesofComplex */


void Negz( z, rz)   /* rz = -z */
COMPLEX z,*rz;
{
  if( z.x != 0.0) rz->x = -z.x;
  if( z.y != 0.0) rz->y = -z.y;
}   /* end. Negz */ 


double Rez( z) 
COMPLEX z;
{
  return( z.x);
}   /* end. Rez */


double Imz( z)
COMPLEX z;
{
  return( z.y);
}   /* end. Imz */


double Magz( z)  /* Magz(z)=|rz| */
COMPLEX z;
{
  
  double abx,aby,temp;
  
  if( z.x == INFINITY || z.y == INFINITY) return( INFINITY);

  abx = fabs( z.x); aby = fabs( z.y);

  if( z.x == 0.0 && z.y == 0.0) return( 0.0); 
  if( z.x == 0.0 && z.y !=0.0) return( aby); 
  if( z.x != 0.0 && z.y == 0.0) return( abx); 
  if( abx!= 0.0 && aby != 0.0){ 
    if( abx >= aby){ 
      temp = aby / abx;
      return( abx*sqrt( 1.0 + temp*temp));
    } 
    else{ 
      temp = abx / aby; 
      return( aby*sqrt( temp*temp + 1.0));
    } 
  } 
} /* end. Magz */


double Magz2( z)
COMPLEX z;
{
  
  double abx,aby,temp; 
  if( z.x == INFINITY || z.y == INFINITY) return( INFINITY);

  abx = fabs( z.x); aby = fabs( z.y);

  if( abx == 0.0 && aby == 0.0) return( 0.0);
  if( abx == 0.0 && aby != 0.0) return( z.y*z.y); 
  if( abx != 0.0 && aby == 0.0) return( z.x*z.x);

  if( abx != 0.0 && aby != 0.0){ 
    if( abx >= aby){ 
      temp = z.y / z.x;
      return( z.x*z.x*( 1.0 + temp*temp));
    } 
    else{ 
      temp = z.x / z.y; 
      return( z.y*z.y*( temp*temp + 1.0));
    } 
  }
} /* end. Magz2 */

double Argz( z) 
COMPLEX z; 
{ 
  return( atan2( z.y, z.x));
} /* end. Argz */


void Addz( z1, z2, rz) /* rz = z1 + z2 */ 
COMPLEX z1,z2,*rz; 
{ 
  rz->x = z1.x + z2.x;
  rz->y = z1.y + z2.y; 
} /* end. Addz */


void Subz( z1, z2, rz) /* rz = z1 - z2 */ 
COMPLEX z1,z2,*rz;
{ 
  rz->x =z1.x - z2.x; 
  rz->y = z1.y - z2.y;
} /* end. Subz */


void Mulz( z1, z2, rz) /* rz = z1 * z2 */ 
COMPLEX z1,z2,*rz; 
{ 
  double abz1,abz2;

  abz1 = Magz( z1); abz2 = Magz( z2);

  if( abz1 == INFINITY){ 
    if( abz2 == INFINITY) InfSetofComplex( rz); /* Infinity * Infinity */ 
    if( abz2 == 0.0) NaNSetofComplex( rz); /* Infinity * 0 */ 
  } 
  if( abz1 == 0.0){ 
    if( abz2 == INFINITY) NaNSetofComplex( rz); /* 0 * Infinity */
    if( abz2 == 0.0) ZeroSetofComplex( rz); /* 0 * 0 */ 
  } 
  else{ 
    rz->x = (z1.x) * (z2.x) - (z1.y) * (z2.y);
    rz->y = (z1.x) * (z2.y) + (z1.y) * (z2.x);
  }
  
} /*end. Mulz */


int Divz( z1, z2, rz) /* rz = z1 / z2 */ 
     COMPLEX z1,z2,*rz; 
{ 
  double abx2,aby2,abz1,abz2,s,temp;

  abz1 = Magz( z1);  abz2 = Magz( z2);
  abx2 = fabs( z2.x);  aby2 = fabs( z2.y); 
  if( abz2 == 0.0){ 
    if( abz1 == 0.0 || abz1 == INFINITY){ 
      NaNSetofComplex( rz);
      return(0);
    } /* in case of 0/0 or Inf/0 */
    else{ 
      InfSetofComplex( rz); 
      return(1);
    } /* in case of z/0 */
  }

  if( abz2 == INFINITY){ 
    if( abz1 == 0.0 || abz1 == INFINITY){ 
      NaNSetofComplex( rz); 
      return(0); 
    } /* in case of 0/Inf or Inf/Inf */
    else{ 
      ZeroSetofComplex( rz); 
      return(1);
    } /* in case of 0/z */
  }

  if( z1.x != 0.0 && z1.y == 0.0 && z2.x != 0.0 && z2.y == 0.0){
    MakeComplex( z1.x / z2.x, 0.0, rz); 
    return(1);
  }

  if( z1.x == 0.0 && z1.y != 0.0 && z2.x == 0.0 && z2.y != 0.0){
    MakeComplex( z1.y / z2.y, 0.0, rz);
    return(1);
  }
  else{
    if( abx2 >= aby2){ 
      temp = z2.y / z2.x; 
      s = z2.x + z2.y * temp;
      rz-> x = ( z1.x + z1.y * temp ) / s;
      rz -> y = ( -z1.x * temp + z1.y ) / s;
    }
    else{
      temp = z2.x / z2.y;
      s = z2.x * temp + z2.y;
      rz -> x = ( z1.x * temp + z1.y ) / s; 
      rz -> y = ( -z1.x + z1.y * temp ) / s;
    } 
  } return(1);

} /* end. Divz */


double sign( x) /* sign returns plus plus, minus or zero. */ 
     double x;
{ 
  if( x == 0.0) return( 0.0);
  if( signbit( x) != 1) return( 1.0);
  if( signbit( x) == 1) return( -1.0);
} /* end. sign */


void Sqrtz( z, rz) /* rz = sqrt( z) */ 
COMPLEX z, *rz;
{
  double r,th,abz;
  
  abz = Magz( z);

  if( abz == INFINITY){
    InfSetofComplex( rz); 
    return;
  }
  if( z.x == 0.0 && z.y == 0.0){
    ZeroSetofComplex( rz);
    return;
  }
 
  if( z.y == 0.0){

    if( z.x > 0.0) MakeComplex( sqrt( z.x), 0.0, rz);
    else MakeComplex( 0.0, sqrt( fabs( z.x)), rz);

  } 
  else{
 
    r = Magz( z);
    th = Argz( z) / 2;

    rz -> x = sign( cos( th)) * sqrt( ( r + z.x)/2);
    rz -> y = sign( sin( th)) * sqrt( ( r - z.x)/2);

  }

}    /* end. Sqrtz */


int powzn( z1, z2, n, rz)    /* rz = pow( z1, z2, n) */
COMPLEX z1,z2,*rz;
int n;
{
  double zz;
  double r,rr,lnr,arg;
  double th,e,sith,coth;

  zz = Magz( z1);

  if( zz == INFINITY){
    InfSetofComplex( rz);
    return(1);
  }

  if( zz == 0.0){
    if( z2.x == 0.0 && z2.y == 0.0){
      NaNSetofComplex( rz);
      return(0);
    }
    if( z2.y != 0.0 ){
      NaNSetofComplex( rz);
      return(0);
    }
    if( z2.x != 0.0 && z2.y == 0.0){
      ZeroSetofComplex( rz);
      return(1);
    }
  }
  else{
    if( z2.x == 0.0 && z2.y == 0.0){
      ESetofComplex( rz);
      return(1);
    }
    else{

      r = zz;
      rr = pow( r, z2.x);
      lnr = log( r);
      arg = Argz( z1) + 2 * n * M_PI;

      e = exp( - z2.y * arg );
      th = z2.y * lnr + z2.x * arg;

      coth = cos( th);
      sith = sin( th);

      if( rr * e == INFINITY) InfSetofComplex( rz);
      else{
	rz -> x = rr * e *coth;
	rz -> y = rr * e *sith;
       }
      return(1);

    }
  }

}  /* end. powzn */


int Powz( z1, z2, rz)   /* rz = pow( z1, z2) */
COMPLEX z1,z2,*rz;
{

  return( powzn( z1, z2, 0, rz));

}  /* end. powz */


void sinz( z, rz)    /* rz = sin( z) */
COMPLEX z,*rz;
{
  double ch, sh;

  ch = cosh( z.y);  sh = sinh( z.y);

  if( ch == INFINITY || sh == INFINITY) InfSetofComplex( rz);
  else{
    rz -> x = sin( z.x) * cosh( z.y);
    rz -> y = cos( z.x) * sinh( z.y);
  }
}  /* end. sinz */


void cosz( z, rz)     /* rz = cos( z) */
COMPLEX z,*rz;
{
  double ch,sh;

  ch = cosh( z.y); sh = sinh( z.y);

  if( ch == INFINITY || sh == INFINITY) InfSetofComplex( rz);
  else{
    rz -> x = cos( z.x) * cosh( z.y);
    rz -> y = -sin( z.x) * sinh( z.y);
  }
}  /* end. cosz */


int tanz( z, rz)     /* rz = tan( z) */ 
COMPLEX z,*rz;
{
  COMPLEX s,c;

  if( z.x == 0.0 && z.y == 0.0){
    ZeroSetofComplex( rz);
    return(1);
  }
  if( z.x != 0.0 && z.y == 0.0){
    MakeComplex( tan( z.x), 0.0, rz);
    return(1);
  }
  if( z.x == 0.0 && z.y != 0.0){
    MakeComplex( 0.0, tanh( z.y), rz);
    return(1);
  }
  if( z.x != 0.0 && z.y != 0.0){

    sinz( z, &s);
    cosz( z, &c);

    if( Divz( s, c, rz) == 0){
      NaNSetofComplex( rz);
      return(0);
    }
    else return(1);
  }

}  /* end. tanz */


void expz( z, rz)     /* rz = exp( z) */ 
COMPLEX z,*rz;
{
  double ex;

  ex = exp( z.x);
  rz -> x = ex*cos( z.y);
  rz -> y = ex*sin( z.y);

}  /* end. expz */ 


void sinhz( z, rz)  /* rz = sinhz( z) */
COMPLEX z,*rz;
{
  double sh, ch;

  sh = sinh( z.x); ch = cosh( z.x);
  if( sh == INFINITY || ch == INFINITY) InfSetofComplex( rz);
  else{
    rz -> x = sh * cos( z.y);
    rz -> y = ch * sin( z.y);
  }
}  /* end. sinhz */


void coshz( z, rz)  /* rz = cosh( z) */
COMPLEX z,*rz;
{
  double sh, ch;

  sh = sinh( z.x); ch = cosh( z.x);
  if( sh == INFINITY || ch == INFINITY) InfSetofComplex( rz);
  else{
    rz -> x = ch * cos( z.y);
    rz -> y = sh * sin( z.y);
  }
}  /* end. coshz */


int tanhz( z, rz)  /* rz = tanh( z) */
COMPLEX z,*rz;
{
  COMPLEX sh,ch;

  if( z.x == 0.0 && z.y == 0.0){
    ZeroSetofComplex( rz);
    return(1);
  }
  if( z.x != 0.0 && z.y == 0.0){
    MakeComplex( tanh( z.x), 0.0, rz);
    return(1);
  }
  if( z.x == 0.0 && z.y != 0.0){
    MakeComplex( 0.0, tan( z.y), rz);
    return(1);
  }
  if( z.x != 0.0 && z.y != 0.0){
 
    sinhz( z, &sh);
    coshz( z, &ch);

    if( Divz( sh, ch, rz) == 0){
      NaNSetofComplex( rz);
      return(0);
    }
    else return(1);
  }

}  /* end. tanhz */

int Logz( z, rz)  /* rz = log( z) */
COMPLEX z,*rz;
{
  double zz,th;

  zz = Magz( z);

  if( zz == 0.0){
    NaNSetofComplex( rz);
    return(0);
  }
  else{
    rz -> x = log( zz);
    th = Argz( z);
    rz -> y = th;
    return(1);
  }

}  /* end. Logz */ 


int Arcsinz( z, rz) /* rz=Arcsin(z)=-i*Log{i*z+sqrt(1-z*z)} */
COMPLEX z,*rz;
{
  COMPLEX tempz1,tempz2;
  COMPLEX e,i;

  ESetofComplex( &e);
  ISetofComplex( &i);
  
  Mulz( z, z, &tempz1);
  Subz( e, tempz1, &tempz1);
  Sqrtz( tempz1, &tempz1);

  Mulz( i, z, &tempz2);
  Addz( tempz1, tempz2, &tempz1);
  if( Logz( tempz1, &tempz1) == 1){
    Mulz( i, tempz1, &tempz1);
    Negz( tempz1, rz);
    return(1);
  }
  else{
    NaNSetofComplex( rz);
    return(0);
  }
  
}   /* end. Arcsinz */


int Arccosz( z, rz)  /* rz=Arccos(z)=-i*Log{z+sqrt(z*z-1)} */  
COMPLEX z,*rz;
{
  COMPLEX temp;
  COMPLEX e,i;

  ESetofComplex( &e);
  ISetofComplex( &i);

  Mulz( z, z, &temp);
  Subz( temp, e, &temp);
  Sqrtz( temp, &temp);

  Addz( z, temp, &temp);
  
  if( Logz( temp, &temp) == 1){
    Mulz( i, temp, &temp); 
    Negz( temp, rz);
    return(1);
  }
  else{
    NaNSetofComplex( rz);
    return(0);
  }

}   /* end. Arccosz */


int Arctanz( z, rz)   /* rz=Arctan(z)=i/2*Log{(1-i*z)/(1+i*z)} */
COMPLEX z,*rz;
{
  COMPLEX temp0,temp1,temp2,temp3;
 
  COMPLEX e,i,halfi;

  if( z.x == 0.0 && z.y == 0.0){
    ZeroSetofComplex( rz);
    return(1);
  }
  if( z.x != 0.0 && z.y == 0.0){
    MakeComplex( atan( z.x), 0.0, rz);
    return(1);
  }
  if( z.x == 0.0 && z.y != 0.0){

    temp0.x = Imz( z);  temp0.y = 0.0;
    ISetofComplex( &i);
    Arctanhz( temp0, &temp0);

    Mulz( i, temp0, rz);
    return(1);
  } 
  if( z.x != 0.0 && z.y != 0.0){
    
    ESetofComplex( &e);
    ISetofComplex( &i);
    MakeComplex( 0.0, 0.5, &halfi);

    Mulz( i, z, &temp0);
    Subz( e, temp0, &temp1);
    Addz( e, temp0, &temp2);

    if ( Divz( temp1, temp2, &temp3) == 1 && Logz( temp3, &temp3) == 1)
      Mulz( halfi, temp3, rz);
    else{
      NaNSetofComplex( rz);
      return(0);
    }
  }
  return(1);

}   /* end. Arctanz */


int Arcsinhz( z, rz)    /* rz=Arcsinh(z)=Log{z+sqrt(z*z+1)} */
COMPLEX z,*rz;
{
  COMPLEX temp;
  COMPLEX e;

  ESetofComplex( &e);

  Mulz( z, z, &temp);
  Addz( temp, e, &temp);
  Sqrtz( temp, &temp);
  Addz( z, temp, &temp);

  if( Logz( temp, rz) == 1) return(1);
  else{
    NaNSetofComplex( rz);
    return(0);
  }

}  /* end. Arcsinhz */

int Arccoshz( z, rz)   /* rz=Arccosh(z)=Log(z+sqrt(z*z-1)} */
COMPLEX z,*rz;
{
  COMPLEX tempz;
  COMPLEX e;

  ESetofComplex( &e);

  Mulz( z, z, &tempz);
  Subz( tempz, e, &tempz);
  Sqrtz( tempz, &tempz);
  Addz( z, tempz, &tempz);
  if ( Logz( tempz, rz) == 1) return(1);
  else{
    NaNSetofComplex( rz);
    return(0);
  }

}   /* end. Arccoshz */

int Arctanhz( z, rz)     /* rz=Arctanh(z)=1/2*Log{(1+z)/(1-z)} */
COMPLEX z,*rz;
{
  COMPLEX temp0,temp1,temp2;
  COMPLEX e,halfe;
 
  if( z.x == 0.0){
    MakeComplex( 0.0, atan( z.y), rz);
    return(1);
  }
  else{

    if( fabs( z.x) <= 1.0 && z.y == 0.0){
      MakeComplex( atanh( z.x), 0.0, rz);
      return(1);
    } 
    else{

      ESetofComplex( &e);
      MakeComplex( 0.5, 0.0, &halfe);

      Addz( e, z, &temp0);
      Subz( e, z, &temp1);
      if( Divz( temp0, temp1, &temp2) == 1 && Logz( temp2, &temp2) == 1){
	Mulz( halfe, temp2, rz);
	return(1);
      }
      else{
	NaNSetofComplex( rz);
	return(0);
      }  
    }
  }

}   /* end. Arctanhz */
    
/* complexmath.c EOF */