1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2c Copyright (C) INRIA 3c 4c Copyright (C) 2012 - 2016 - Scilab Enterprises 5c 6c This file is hereby licensed under the terms of the GNU GPL v2.0, 7c pursuant to article 5.3.4 of the CeCILL v.2.1. 8c This file was originally licensed under the terms of the CeCILL v2.1, 9c and continues to be available under such terms. 10c For more information, see the COPYING file which you should have received 11c along with this program. 12c 13 subroutine wasin(zr, zi, ar, ai) 14* 15* PURPOSE 16* Compute the arcsin of a complex number 17* a = ar + i ai = asin(z), z = zr + i zi 18* 19* CALLING LIST / PARAMETERS 20* subroutine wasin(zr,zi,ar,ai) 21* double precision zr,zi,ar,ai 22* 23* zr,zi: real and imaginary parts of the complex number 24* ar,ai: real and imaginary parts of the result 25* ar,ai may have the same memory cases than zr et zi 26* 27* REFERENCE 28* This is a Fortran-77 translation of an algorithm by 29* T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which 30* appears in their article : 31* "Implementing the Complex Arcsine and Arccosine 32* Functions Using Exception Handling", ACM, TOMS, 33* Vol 23, No. 3, Sept 1997, p. 299-335 34* 35* with some modifications so as don't rely on ieee handle 36* trap functions. 37* 38* AUTHOR 39* Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr> 40* Thanks to Tom Fairgrieve 41* 42 implicit none 43 44* PARAMETERS 45 double precision zr, zi, ar, ai 46 47* EXTERNAL FUNCTIONS 48 double precision dlamch, logp1 49 external dlamch, logp1 50 integer isanan 51 external isanan 52 53* CONSTANTS 54 double precision LN2, HALFPI, Across, Bcross 55 parameter (LN2 = 0.6931471805599453094172321d0, 56 $ HALFPI = 1.5707963267948966192313216d0, 57 $ Across = 1.5d0, 58 $ Bcross = 0.6417d0) 59 60* LOCAL VARIABLES 61 double precision x, y, A, B, R, S, Am1, szr, szi 62 63 64* STATIC VARIABLES 65 double precision LSUP, LINF, EPSM 66 save LSUP, LINF, EPSM 67 logical first 68 save first 69 data first /.true./ 70* TEXT 71* got f.p. parameters used by the algorithm 72 if (first) then 73 LSUP = sqrt(dlamch('o'))/8.d0 74 LINF = 4.d0*sqrt(dlamch('u')) 75 EPSM = sqrt(dlamch('e')) 76 first = .false. 77 endif 78 79* avoid memory pb ... 80 x = abs(zr) 81 y = abs(zi) 82 szr = sign(1.d0,zr) 83 szi = sign(1.d0,zi) 84 85 86 if (LINF .le. min(x,y) .and. max(x,y) .le. LSUP ) then 87* we are in the safe region 88 R = sqrt((x+1.d0)**2 + y**2) 89 S = sqrt((x-1.d0)**2 + y**2) 90 A = 0.5d0*(R + S) 91 B = x/A 92 93* compute the real part 94 if ( B .le. Bcross ) then 95 ar = asin(B) 96 elseif ( x .le. 1.d0 ) then 97 ar = atan( x / sqrt( 98 $ 0.5d0*(A+x)*((y**2)/(R+(x+1.d0))+(S+(1.d0-x)))) ) 99 else 100 ar = atan( x / 101 $ (y*sqrt(0.5d0*((A+x)/(R+(x+1.d0))+(A+x)/(S+(x-1.d0)))))) 102 endif 103 104* compute the imaginary part 105 if ( A .le. Across ) then 106 if ( x .lt. 1.d0 ) then 107 Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) 108 else 109 Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) 110 endif 111 ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) 112 else 113 ai = log(A + sqrt(A**2 - 1.d0)) 114 endif 115 116 else 117* HANDLE BLOC : evaluation in the special regions ... 118 if ( y .le. EPSM*abs(x-1.d0) ) then 119 if (x .lt. 1.d0 ) then 120 ar = asin(x) 121 ai = y/sqrt((1.d0+x)*(1.d0-x)) 122 else 123 ar = HALFPI 124 if ( x .le. LSUP ) then 125 ai = logp1((x-1.d0) + sqrt((x-1.d0)*(x+1.d0))) 126 else 127 ai = LN2 + log(x) 128 endif 129 endif 130 131 elseif (y .lt. LINF) then 132 if (isanan(x).eq.1) then 133 ar = x 134 else 135 ar = HALFPI - sqrt(y) 136 endif 137 ai = sqrt(y) 138 139 elseif (EPSM*y - 1.d0 .ge. x) then 140 ar = x/y 141 ai = LN2 + log(y) 142 143 elseif (x .gt. 1.d0) then 144 ar = atan(x/y) 145 ai = LN2 + log(y) + 0.5d0*logp1((x/y)**2) 146 147 else 148 A = sqrt(1 + y**2) 149 ar = x/A 150 ai = 0.5*logp1(2.d0*y*(y+A)) 151 endif 152 endif 153 154* recover the signs 155 ar = szr * ar 156 if (y.eq.0d00 .and. szr.gt.0d00) then 157 szi = - szi 158 endif 159 ai = szi * ai 160 161 end 162 163 164