//**********************************************************************************// // Copyright (C) 2009-2015 Ovidio Pena // // // // This file is part of scattnlay // // // // This program is free software: you can redistribute it and/or modify // // it under the terms of the GNU General Public License as published by // // the Free Software Foundation, either version 3 of the License, or // // (at your option) any later version. // // // // This program 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 General Public License for more details. // // // // The only additional remark is that we expect that all publications // // describing work using this software, or all commercial products // // using it, cite the following reference: // // [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // // a multilayered sphere," Computer Physics Communications, // // vol. 180, Nov. 2009, pp. 2348-2354. // // // // You should have received a copy of the GNU General Public License // // along with this program. If not, see . // //**********************************************************************************// #include #include #include #include "ucomplex.h" complex Cadd(complex z1, complex z2) { complex zr; zr.r = z1.r + z2.r; zr.i = z1.i + z2.i; return zr; } complex Csub(complex z1, complex z2) { complex zr; zr.r = z1.r - z2.r; zr.i = z1.i - z2.i; return zr; } complex Cmul(complex z1, complex z2) { complex zr; zr.r = (z1.r*z2.r) - (z1.i*z2.i); zr.i = (z1.r*z2.i) + (z1.i*z2.r); return zr; } complex RCmul(double r, complex z) { complex zr; zr.r = z.r*r; zr.i = z.i*r; return zr; } complex Cdiv(complex z1, complex z2) { /* The following algorithm is used to properly handle denominator overflow: | a + b(d/c) b - a(d/c) | ---------- + ----------- I if |d| < |c| a + b I | c + d(d/c) c + d(d/c) ------- = | c + d I | b + a(c/d) -a + b(c/d) | ---------- + ----------- I if |d| >= |c| | d + c(c/d) d + c(c/d) */ complex zr; double tmp, denom; if(fabs(z2.r) > fabs(z2.i)) { tmp = z2.i/z2.r; denom = z2.r + z2.i*tmp; zr.r = (z1.r + z1.i*tmp)/denom; zr.i = (z1.i - z1.r*tmp)/denom; } else { tmp = z2.r/z2.i; denom = z2.i + z2.r*tmp; zr.r = (z1.i + z1.r*tmp)/denom; zr.i = (-z1.r + z1.i*tmp)/denom; } return zr; } complex Complex(double re, double im) { complex zr; zr.r = re; zr.i = im; return zr; } double Cabs(complex z) { double r; r = sqrt((z.r*z.r) + (z.i*z.i)); return r; } double Carc(complex z) { double r; r = atan2(z.i, z.r); return r; } complex Conjg(complex z) { complex zr; zr.r = z.r; zr.i = -z.i; return zr; } complex Cinv(complex z) { complex zr; double denom; denom = (z.r*z.r) + (z.i*z.i); zr.r = z.r/denom; zr.i = -z.i/denom; return zr; } complex Cexp(complex z) { // exp(a + ib) = exp(a)*exp(ib) = exp(a)*[cos(b) + i sin(b)] complex zr; double expz; expz = exp(z.r); zr.r = expz*cos(z.i); zr.i = expz*sin(z.i); return zr; } complex Clog(complex z) { // ln( p exp(i0)) = ln(p) + i0 + 2kpi complex zr; zr.r = log(Cabs(z)); zr.i = atan2(z.i, z.r); return zr; } complex Csqrt(complex z) { complex zr; double root, q; if((z.r != 0.0) ||(z.i != 0.0)) { root = sqrt(0.5*(fabs(z.r) + Cabs(z))); q = z.i/(2.0*root); if(z.r >= 0.0) { zr.r = root; zr.i = q; } else if(z.i < 0.0) { zr.r = -q; zr.i = -root; } else { zr.r = q; zr.i = root; } } else zr = z; return zr; } // complex trigonometric functions // complex cosinus complex Ccos(complex z) { // cos(x+iy) = cos(x).cos(iy) - sin(x).sin(iy) // cos(ix) = cosh(x) et sin(ix) = i.sinh(x) complex zr; zr.r = cos(z.r)*cosh(z.i); zr.i = -sin(z.r)*sinh(z.i); return zr; } // complex sinus complex Csin(complex z) { // sin(x+iy) = sin(x).cos(iy) + cos(x).sin(iy) // cos(ix) = cosh(x) et sin(ix) = i.sinh(x) complex zr; zr.r = sin(z.r)*cosh(z.i); zr.i = cos(z.r)*sinh(z.i); return zr; } // tangent complex Ctan(complex z) { return Cdiv(Csin(z), Ccos(z)); } // Inverse trigonometric functions // complex arc cosinus complex Carc_cos(complex z) { // arccos(z) = -i.arcch(z) complex temp = Carc_ch(z); return Complex(temp.i, -temp.r); } // complex arc sinus complex Carc_sin(complex z) { // arcsin(z) = -i.arcsh(i.z) complex temp = Carc_sh(Complex(-z.i, z.r)); return Complex(temp.i, -temp.r); } // complex arc tangent complex Carc_tan(complex z) { // arctg(z) = -i.arcth(i.z) complex temp = Carc_th(Complex(-z.i, z.r)); return Complex(temp.i, -temp.r); } // Hyberbolic complex functions // complex hyberbolic cosinus complex Cch(complex z) { // cosh(x+iy) = cosh(x).cosh(iy) + sinh(x).sinh(iy) // cosh(iy) = cos(y) et sinh(iy) = i.sin(y) complex zr; zr.r = cosh(z.r)*cos(z.i); zr.i = sinh(z.r)*sin(z.i); return zr; } // complex hyberbolic sinus complex Csh(complex z) { // sinh(x+iy) = sinh(x).cosh(iy) + cosh(x).sinh(iy) // cosh(iy) = cos(y) et sinh(iy) = i.sin(y) complex zr; zr.r = sinh(z.r)*cos(z.i); zr.i = cosh(z.r)*sin(z.i); return zr; } // complex hyberbolic tangent complex Cth(complex z) { // th(x) = sinh(x)/cosh(x) // cosh(x) > 1 qq x complex temp, zr; temp = Cch(z); zr = Csh(z); return Cdiv(zr, temp); } // Inverse complex hyperbolic functions // complex hyperbolic arc cosinus complex Carc_ch(complex z) { // _________ // arcch(z) = -/+ ln(z + V z.z - 1 ) complex zr; complex z2 = Cmul(z, z); zr = Clog(Cadd(z, Csqrt(Complex(z2.r - 1.0, z2.i)))); zr.r = -zr.r; zr.i = -zr.i; return zr; } // complex hyperbolic arc sinus complex Carc_sh(complex z) { // ________ // arcsh(z) = ln(z + V 1 + z.z) complex z2 = Cmul(z, z); return Clog(Cadd(z, Csqrt(Complex(z2.r + 1.0, z2.i)))); } // complex hyperbolic arc tangent complex Carc_th(complex z) { // arcth(z) = 1/2 ln((z + 1)/(1 - z)) return RCmul(0.5, Clog(Cdiv(Complex(z.r + 1.0, z.i), Complex(1.0 - z.r, -z.i)))); }