1 files changed, 49 insertions, 0 deletions
diff --git a/src/tools.jl b/src/tools.jl
new file mode 100644
index 0000000..0d3dc7f
--- /dev/null
+++ b/src/tools.jl
@@ -0,0 +1,49 @@
+## Copyright 2021 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)
+  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)
+  #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)
+  out=zeros(Float64,length(v))
+  for i in 1:length(v)
+    out[i]=Phi(v[i])
+  end
+  return out
+end
+@everywhere function dotdPhi(v)
+  out=zeros(Float64,length(v))
+  for i in 1:length(v)
+    out[i]=dPhi(v[i])
+  end
+  return out
+end