From e72af82c3ed16b81cdb5043c58abbdbb3cf02102 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian@jauslin.org>
Date: Mon, 4 Oct 2021 11:12:34 -0400
Subject: Initial commit

---
 src/tools.jl | 49 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 49 insertions(+)
 create mode 100644 src/tools.jl

(limited to 'src/tools.jl')

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
-- 
cgit v1.2.3-54-g00ecf