/*********************************************
* Elementary operations with complex numbers *
* ------------------------------------------ *
* See demo program tcomplex.cpp.             *
*                                            *
*          C++ version by J-P Moreau, Paris. *
*                  (www.jpmoreau.fr)         *
*********************************************/
#include "complex2.h"


  // define complex number by x and y
  void AssignXY(COMPLEX *Z, double x, double y) {
    Z->a=x; Z->b=y;
    Z->r=sqrt(x*x+y*y);
    if (x==0) {
      if (y>0) Z->t=PI/2;
      else if (y==0) Z->t=-PI/2; 
	  else Z->t=0;
    }
    else {
      Z->t=atan(y/x);
      if (x<0) {
        if (y>=0) Z->t=Z->t+PI; 
		else Z->t=Z->t-PI;
      }
      if (Z->t>PI) Z->t=Z->t-2*PI;
      if (Z->t<-PI) Z->t=Z->t+2*PI;
    }
  }

  // define complex number by r and t in radians
  void AssignRT(COMPLEX *Z, double r, double t) {
    Z->r=r; Z->t=t;
    if (Z->t>PI)  Z->t=Z->t-2*PI;
    if (Z->t<-PI) Z->t=Z->t+2*PI;
    Z->a=r*cos(t);
    Z->b=r*sin(t);
  }

  // display complex number with x and y
  void DisplayComplex(COMPLEX Z) {
    printf("(%f,%f)", Z.a,Z.b);
  }

  // display complex number with radius and phase in radians
  void DisplayComplexR(COMPLEX Z) {
    printf("(%f,%f)", Z.r,Z.t);
  }

  // add two complex numbers: Z2=Z+Z1
  void AddComplex(COMPLEX *Z2, COMPLEX Z, COMPLEX Z1) {
    Z2->a=Z.a+Z1.a; Z2->b=Z.b+Z1.b;
    AssignXY(Z2,Z2->a,Z2->b);
  }

  // substract two complex numbers: Z2=Z-Z1
  void SubComplex(COMPLEX *Z2, COMPLEX Z, COMPLEX Z1) {
    Z2->a=Z.a-Z1.a; Z2->b=Z.b-Z1.b;
    AssignXY(Z2,Z2->a,Z2->b);
  }

  // change sign of a complex number
  void ChsComplex(COMPLEX *Z1, COMPLEX Z) {
    Z1->a=-Z.a; Z1->b=-Z.b;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // multiply two complex numbers: Z2=Z*Z1
  void MulComplex(COMPLEX *Z2, COMPLEX Z, COMPLEX Z1) {
    Z2->r=Z.r*Z1.r; Z2->t=Z.t+Z1.t;
    AssignRT(Z2,Z2->r,Z2->t);
  }

  // divide two complex numbers: Z2=Z/Z1
  void DivComplex(COMPLEX *Z2, COMPLEX Z, COMPLEX Z1) {
    if (Z1.r < TINY)  Z2->r = INF;
    else Z2->r = Z.r/Z1.r;
    Z2->t=Z.t-Z1.t;
    AssignRT(Z2,Z2->r,Z2->t);
  }

  // exponential complex function: Z1=Exp(Z)
  void ExpComplex(COMPLEX *Z1, COMPLEX Z) {
    double temp;
    if (exp(Z.a) > INF) temp=INF; else temp=exp(Z.a); 
    Z1->a=temp*cos(Z.b); Z1->b=temp*sin(Z.b);
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=LN(Z)
  void LnComplex(COMPLEX *Z1, COMPLEX Z) {
    if (Z.r<=0) Z1->a=-INF;
    else Z1->a=log(Z.r);
    Z1->b=Z.t;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z2=Z^Z1
  void PowerComplex(COMPLEX *Z2, COMPLEX Z, COMPLEX Z1) {
    COMPLEX temp;
    LnComplex(&temp,Z);
    MulComplex(&temp,Z1,temp);
    ExpComplex(Z2,temp);
  }

  // Z1=COS(Z)
  void CosComplex(COMPLEX *Z1, COMPLEX Z) {
    Z1->a=(exp(-Z.b)*cos(Z.a) + exp(Z.b)*cos(-Z.a))/2;
    Z1->b=(exp(-Z.b)*sin(Z.a) + exp(Z.b)*sin(-Z.a))/2;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=SIN(Z)
  void SinComplex(COMPLEX *Z1, COMPLEX Z) {
    Z1->a= (exp(-Z.b)*sin(Z.a) - exp(Z.b)*sin(-Z.a))/2;
    Z1->b=-(exp(-Z.b)*cos(Z.a) - exp(Z.b)*cos(-Z.a))/2;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=TAN(Z)
  void TanComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX z2,z3;
    SinComplex(&z2,Z);
    CosComplex(&z3,Z);
    DivComplex(Z1,z2,z3);
  }

  // Z1=CH(Z)
  void ChComplex(COMPLEX *Z1, COMPLEX Z) {
    Z1->a=(exp(Z.a)*cos(Z.b) + exp(-Z.a)*cos(-Z.b))/2;
    Z1->b=(exp(Z.a)*sin(Z.b) + exp(-Z.a)*sin(-Z.b))/2;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=SH(Z)
  void ShComplex(COMPLEX *Z1, COMPLEX Z) {
    Z1->a= (exp(Z.a)*cos(Z.b) - exp(-Z.a)*cos(-Z.b))/2;
    Z1->b=-(exp(Z.a)*sin(Z.b) - exp(-Z.a)*sin(-Z.b))/2;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=TH(Z)
  void ThComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX z2,z3;
    ShComplex(&z2,Z);
    ChComplex(&z3,Z);
    DivComplex(Z1,z2,z3);
  }

  // Z1=ARCCOS(Z)
  void ArcCosComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX temp,u;
    temp.a=1-Z.a*Z.a+Z.b*Z.b;
    temp.b=-2*Z.a*Z.b;
    AssignXY(&temp,temp.a,temp.b);
    AssignXY(&u,0.5,0.0);
    PowerComplex(&temp,temp,u);
    temp.a=Z.a-temp.b;
    temp.b=Z.b+temp.a;
    AssignXY(&temp,temp.a,temp.b);
    LnComplex(&temp,temp);
    Z1->a=temp.b; Z1->b=-temp.a;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=ARCSIN(Z)
  void ArcSinComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX temp,u;
    temp.a=1-Z.a*Z.a+Z.b*Z.b;
    temp.b=-2*Z.a*Z.b;
    AssignXY(&temp,temp.a,temp.b);
    AssignXY(&u,0.5,0.0);
    PowerComplex(&temp,temp,u);
    temp.a=temp.a-Z.b;
    temp.b=temp.b+Z.a;
    AssignXY(&temp,temp.a,temp.b);
    LnComplex(&temp,temp);
    Z1->a=temp.b; Z1->b=-temp.a;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=ARCTAN(Z)
  void ArcTanComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX z2,z3;
    z2.a=-Z.a; z2.b=1-Z.b;
    z3.a=Z.a; z3.b=1+Z.b;
    AssignXY(&z2,z2.a,z2.b);
    AssignXY(&z3,z3.a,z3.b);
    DivComplex(&z3,z2,z3);
    LnComplex(&z3,z3);
    Z1->a=z3.b/2; Z1->b=-z3.a/2;
    AssignXY(Z1,Z1->a,Z1->b);
  }

  // Z1=ARGCH(Z)
  void ArgChComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX temp,u;
    temp.a=-1+Z.a*Z.a-Z.b*Z.b;
    temp.b=2*Z.a*Z.b;
    AssignXY(&temp,temp.a,temp.b);
    AssignXY(&u,0.5,0.0);
    PowerComplex(&temp,temp,u);
    temp.a=temp.a+Z.a;
    temp.b=temp.b+Z.b;
    AssignXY(&temp,temp.a,temp.b);
    LnComplex(Z1,temp);
  }

  // Z1=ARGSH(Z)
  void ArgShComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX temp,u;
    temp.a=1+Z.a*Z.a-Z.b*Z.b;
    temp.b=2*Z.a*Z.b;
    AssignXY(&temp,temp.a,temp.b);
    AssignXY(&u,0.5,0.0);
    PowerComplex(&temp,temp,u);
    temp.a=temp.a+Z.a;
    temp.b=temp.b+Z.b;
    AssignXY(&temp,temp.a,temp.b);
    LnComplex(Z1,temp);
  }

  // Z1=ARGTH(Z)
  void ArgThComplex(COMPLEX *Z1, COMPLEX Z) {
    COMPLEX z2,z3;
    z2.a=1+Z.a; z2.b=Z.b;
    z3.a=1-Z.a; z3.b=-Z.b;
    AssignXY(&z2,z2.a,z2.b);
    AssignXY(&z3,z3.a,z3.b);
    DivComplex(&z3,z2,z3);
    LnComplex(&z3,z3);
    Z1->a=z3.a/2; Z1->b=z3.b/2;
    AssignXY(Z1,Z1->a,Z1->b);
  }


// end of file complex2.cpp