1 files changed, 75 insertions, 0 deletions
diff --git a/src/hh_root_double.c b/src/hh_root_double.c
new file mode 100644
index 0000000..5474bcb
--- /dev/null
+++ b/src/hh_root_double.c
@@ -0,0 +1,75 @@
+/*
+Copyright 2016 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.
+*/
+
+#include "hh_root_double.h"
+#include <math.h>
+#include "hh_integral_double.h"
+
+// wrapper for the integration function, used for the Newton scheme
+int hh_integration_wrapper_double(long double* out, long double in, void* args){
+  hh_args_integration_double* argument=(hh_args_integration_double*)args;
+  long double val;
+  int ret;
+
+  argument->params.W=in;
+
+  ret=hh_integrate_double(&val, argument->params, argument->abcissa, argument->weights);
+
+  if(ret<0){
+    return(ret);
+  }
+
+  // repeat with -sinphi
+  argument->params.sinphi=-argument->params.sinphi;
+
+  ret=hh_integrate_double(out, argument->params, argument->abcissa, argument->weights);
+
+  // reset sinphi
+  argument->params.sinphi=-argument->params.sinphi;
+
+  *out=argument->params.W+3*sqrtl(3)*argument->params.omega*argument->params.t2*argument->params.sinphi-argument->params.lambda*sqrtl(3)/6*(*out+val);
+
+  return(ret);
+}
+
+// wrapper for the derivative of the integration function, used for the Newton scheme
+int hh_d_integration_wrapper_double(long double* out, long double in, void* args){
+  hh_args_integration_double* argument=(hh_args_integration_double*)args;
+  long double val;
+  int ret;
+
+  argument->params.W=in;
+
+  ret=hh_d_integrate_double(&val, argument->params, argument->abcissa, argument->weights);
+
+  if(ret<0){
+    return(ret);
+  }
+
+  // repeat with -sinphi
+  argument->params.sinphi=-argument->params.sinphi;
+  argument->params.phi=-argument->params.phi;
+
+  ret=hh_d_integrate_double(out, argument->params, argument->abcissa, argument->weights);
+
+  // reset sinphi
+  argument->params.sinphi=-argument->params.sinphi;
+  argument->params.phi=-argument->params.phi;
+
+  *out=1-argument->params.lambda*sqrtl(3)/6*(*out+val);
+
+  return(ret);
+}