1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
// complex.cmd
// Surface Evolver command file.
// Programmer: Ken Brakke, brakke@susqu.edu, http://www.susqu.edu/brakke
/* Some complex functions of complex arguments, mostly for use with i
Weierstrass representation of minimal surfaces.
Contents:
re_sqrt, im_sqrt
re_sin, im_sin
re_arcsin, im_arcsin
re_incompleteEllipticF, im_incompleteEllipticF
Each complex function is implemented as two real functions, since
Evolver doesn't have a complex or array return type for functions,
and using global variables for returns would not be thread safe.
Arguments:
rex,imx: real and imaginary parts of the complex argument
But see re_incompleteEllipticF for its further arguments.
*/
// Sqrt, cut on negative real axis
function real re_sqrt(real rex, real imx)
{
return sqrt((rex + sqrt(rex^2+imx^2))/2);
}
function real im_sqrt(real rex, real imx)
{
if imx > 0 then
return sqrt((-rex + sqrt(rex^2+imx^2))/2)
else
return -sqrt((-rex + sqrt(rex^2+imx^2))/2);
}
// Sin
function real re_sin(real rex, real imx)
{
return sin(rex)*cosh(imx);
}
function real im_sin(real rex, real imx)
{
return cos(rex)*sinh(imx);
}
// Arcsin; multiple valued - (x+2*pi*k,y),(pi-x+2*pi*k,-y)
function real re_arcsin(real rex, real imx)
{
local term;
term := rex^2 + imx^2 + 1;
if rex < 0 then
return -asin(sqrt((term-sqrt(term^2-4*rex^2))/2))
else
return asin(sqrt((term-sqrt(term^2-4*rex^2))/2));
}
function real im_arcsin(real rex, real imx)
{
local term,mag,sign;
mag := rex^2 + imx^2;
term := mag - 1;
if imx < 0 then sign := -1
else sign := 1;
if mag < 0.5 then // alternate forms to avoid catastrophic cancellation
return sign*asinh(sqrt(2*imx^2/(sqrt(term^2 + 4*imx^2) - term)))
else
return sign*asinh(sqrt((term+sqrt(term^2+4*imx^2))/2));
}
// incompleteEllipticF, Abramowitz and Stegun 17.4.11
// Have to beware branch points at +/- arcsin(1/sqrt(m)) + 2*pi*k,
// so nbr_value input argument for picking the
// proper branch by continuity; if 0, then principle branch
// is picked. Branch values differ by 2*ellipticK(m).
// nbr_test is boolean for whether to apply the continuity test.
function real re_incompleteEllipticF(real rex, real imx, real m_param,
real nbr_value,integer nbr_test)
{
local term,cotsq,lambda;
if abs(sin(rex)) < 1e-12 then newvalue := 0.0
else
{
rrex := rex;
ss := floor(rex/pi + .5);
rex -= ss*pi;
term := 1/tan(rex)^2 + m_param*sinh(imx)^2/sin(rex)^2 - (1-m_param);
cotsq := (term + sqrt(term^2 + 4*(1-m_param)/tan(rex)^2))/2;
if cotsq == 0 then lambda := pi/2
else lambda := atan(1/sqrt(cotsq));
if rex < 0 then lambda := -lambda;
newvalue := incompleteEllipticF(lambda+ss*pi,m_param);
};
if nbr_test then
{ period := 2*ellipticK(m_param);
tt := floor((newvalue-nbr_value)/period + .5);
newvalue -= tt*period;
alt_tt := floor((-newvalue-nbr_value)/period + .5);
alt_newvalue := -newvalue - alt_tt*period;
if ( abs(alt_newvalue - nbr_value) < abs(newvalue - nbr_value) ) then
newvalue := alt_newvalue;
};
return newvalue;
}
function real im_incompleteEllipticF(real rex, real imx, real m_param,
real nbr_value, integer nbr_test)
{
local term,cotsq,mu;
term := cos(rex)^2 + m_param*sinh(imx)^2 - (1-m_param)*sin(rex)^2;
if term > 0 then
mtansqp1 := (term + sqrt(term^2 + 4*(1-m_param)*sin(rex)^2*cos(rex)^2))/
2/cos(rex)^2
else
mtansqp1 := 2*(1-m_param)*sin(rex)^2/
(sqrt(term^2 + 4*(1-m_param)*sin(rex)^2*cos(rex)^2) - term);
tansq := (mtansqp1 - 1)/m_param;
mu := atan(sqrt(tansq));
if imx < 0 then mu := -mu;
return incompleteEllipticF(mu,1-m_param);
}
|
