/********************************************* * 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