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