Ian Jauslin
summaryrefslogtreecommitdiff
blob: 5d1678d2d29311ad450c7bd6509a447360ed8142 (plain)
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
## Copyright 2021-2023 Ian Jauslin
## 
## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License.
## You may obtain a copy of the License at
## 
##     http://www.apache.org/licenses/LICENSE-2.0
## 
## Unless required by applicable law or agreed to in writing, software
## distributed under the License is distributed on an "AS IS" BASIS,
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
## See the License for the specific language governing permissions and
## limitations under the License.

# \Phi(x)=2*(1-sqrt(1-x))/x
@everywhere function Phi(
  x::Float64
)
  if abs(x)>1e-5
    return 2*(1-sqrt(abs(1-x)))/x
  else
    return 1+x/4+x^2/8+5*x^3/64+7*x^4/128+21*x^5/512
  end
end
# \partial\Phi
@everywhere function dPhi(
  x::Float64
)
  #if abs(x-1)<1e-5
  #  @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x))
  #end
  if abs(x)>1e-5
    return 1/(sqrt(abs(1-x))*x)*(1-x>=0 ? 1 : -1)-2*(1-sqrt(abs(1-x)))/x^2
  else
    return 1/4+x/4+15*x^2/64+7*x^3/32+105*x^4/512+99*x^5/512
  end
end

# apply Phi to every element of a vector
@everywhere function dotPhi(
  v::Array{Float64,1}
)
  out=zeros(Float64,length(v))
  for i in 1:length(v)
    out[i]=Phi(v[i])
  end
  return out
end
@everywhere function dotdPhi(
  v::Array{Float64,1}
)
  out=zeros(Float64,length(v))
  for i in 1:length(v)
    out[i]=dPhi(v[i])
  end
  return out
end

@everywhere function sinc(x::Float64)
  return (x == 0 ? 1 : sin(x)/x)
end