1 files changed, 144 insertions, 0 deletions
diff --git a/src/hh_root.c b/src/hh_root.c
new file mode 100644
index 0000000..6df4c82
--- /dev/null
+++ b/src/hh_root.c
@@ -0,0 +1,144 @@
+/*
+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.h"
+
+#include <mpfr.h>
+
+#include "hh_integral.h"
+
+// wrapper for the integration function, used for the Newton scheme
+int integration_wrapper(mpfr_t* out, mpfr_t in, void* extra_args){
+  mpfr_t tmp;
+  hh_params params;
+  int ret;
+
+  mpfr_init(tmp);
+
+  mpfr_set(((args_integration*)extra_args)->params.W, in, MPFR_RNDN);
+  
+  // out = I+
+  ret=hh_integrate(out, ((args_integration*)extra_args)->params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights);
+  if(ret<0){
+    mpfr_clear(tmp);
+    return(ret);
+  }
+
+  // clone params (to set sinphi=-sinphi in order to compute I_-)
+  mpfr_inits(params.t1, params.t2, params.lambda, params.W, params.sinphi, NULL);
+  mpfr_set(params.t1, ((args_integration*)extra_args)->params.t1, MPFR_RNDN);
+  mpfr_set(params.t2, ((args_integration*)extra_args)->params.t2, MPFR_RNDN);
+  mpfr_set(params.lambda, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN);
+  mpfr_set(params.W, ((args_integration*)extra_args)->params.W, MPFR_RNDN);
+  mpfr_neg(params.sinphi, ((args_integration*)extra_args)->params.sinphi, MPFR_RNDN);
+
+  // tmp = I-
+  ret=hh_integrate(&tmp, params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights);
+  if(ret<0){
+    mpfr_clear(tmp);
+    mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, NULL);
+    return(ret);
+  }
+
+  mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, NULL);
+
+  // out=I+ + I-
+  mpfr_add(*out, *out, tmp, MPFR_RNDN);
+  //// tmp free
+
+  // tmp = sqrt(3)
+  mpfr_sqrt_ui(tmp, 3, MPFR_RNDN);
+
+  // out = W-sqrt(3)*lambda/6*(I+ + I-)
+  mpfr_mul(*out, *out, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN);
+  mpfr_div_ui(*out, *out, 6, MPFR_RNDN);
+  mpfr_mul(*out, *out, tmp, MPFR_RNDN);
+  mpfr_sub(*out, ((args_integration*)extra_args)->params.W, *out, MPFR_RNDN);
+
+  // tmp = 3*sqrt(3)*t2*sin(phi)
+  mpfr_mul_ui(tmp, tmp, 3, MPFR_RNDN);
+  mpfr_mul(tmp, tmp, ((args_integration*)extra_args)->params.t2, MPFR_RNDN);
+  mpfr_mul(tmp, tmp, ((args_integration*)extra_args)->params.sinphi, MPFR_RNDN);
+
+  // W+w*3*sqrt(3)*t2*sin(phi)-sqrt(3)*lambda/6*(I+ + I-)
+  if(((args_integration*)extra_args)->params.omega==1){
+    mpfr_add(*out, *out, tmp, MPFR_RNDN);
+  }
+  else{
+    mpfr_sub(*out, *out, tmp, MPFR_RNDN);
+  }
+  //// tmp free
+
+  mpfr_clear(tmp);
+  return(0);
+}
+
+// wrapper for the derivative of the integration function, used for the Newton scheme
+int d_integration_wrapper(mpfr_t* out, mpfr_t in, void* extra_args){
+  mpfr_t tmp;
+  hh_params params;
+  int ret;
+
+  mpfr_init(tmp);
+
+  mpfr_set(((args_integration*)extra_args)->params.W, in, MPFR_RNDN);
+  
+  // out = dI+
+  ret=hh_d_integrate(out, ((args_integration*)extra_args)->params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights);
+  if(ret<0){
+    mpfr_clear(tmp);
+    return(ret);
+  }
+
+  // clone params (to set sinphi=-sinphi in order to compute dI_-)
+  mpfr_inits(params.t1, params.t2, params.lambda, params.W, params.sinphi, params.phi, NULL);
+  mpfr_set(params.t1, ((args_integration*)extra_args)->params.t1, MPFR_RNDN);
+  mpfr_set(params.t2, ((args_integration*)extra_args)->params.t2, MPFR_RNDN);
+  mpfr_set(params.lambda, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN);
+  mpfr_set(params.W, ((args_integration*)extra_args)->params.W, MPFR_RNDN);
+  params.omega=((args_integration*)extra_args)->params.omega;
+  mpfr_neg(params.sinphi, ((args_integration*)extra_args)->params.sinphi, MPFR_RNDN);
+  mpfr_neg(params.phi, ((args_integration*)extra_args)->params.phi, MPFR_RNDN);
+
+  // tmp = dI-
+  ret=hh_d_integrate(&tmp, params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights);
+  if(ret<0){
+    mpfr_clear(tmp);
+    mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, params.phi, NULL);
+    return(ret);
+  }
+
+  mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, params.phi, NULL);
+
+  // out=dI+ + dI-
+  mpfr_add(*out, *out, tmp, MPFR_RNDN);
+  //// tmp free
+
+  // tmp = sqrt(3)
+  mpfr_sqrt_ui(tmp, 3, MPFR_RNDN);
+
+  // out = 1-sqrt(3)*lambda/6*(dI+ + dI-)
+  mpfr_mul(*out, *out, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN);
+  mpfr_div_ui(*out, *out, 6, MPFR_RNDN);
+  mpfr_mul(*out, *out, tmp, MPFR_RNDN);
+  //// tmp free
+  mpfr_ui_sub(*out, 1, *out, MPFR_RNDN);
+
+  mpfr_clear(tmp);
+
+  return(0);
+}
+