From 61ee450ddf82c6206ac3e55dd5e064d5002e57fe Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 23 Nov 2021 16:53:17 +0100
Subject: Fixed some problems with changing the Real type, and went back to
 geometric increases in fit.

---
 stokes.hpp | 23 +++++++++++------------
 1 file changed, 11 insertions(+), 12 deletions(-)

(limited to 'stokes.hpp')

diff --git a/stokes.hpp b/stokes.hpp
index 76b950f..bcb6608 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -34,7 +34,7 @@ Real dJacobi(unsigned α, unsigned β, unsigned n, Real x) {
     f += pow(y, n - s) * pow(z, s - 1) / (tgamma(s - 1 + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s + 1));
   }
 
-  return 0.5 * tgamma(n + α + 1) * tgamma(n + β + 1) * f;
+  return tgamma(n + α + 1) * tgamma(n + β + 1) * f / 2;
 }
 
 template <class Scalar, int ...ps>
@@ -44,7 +44,6 @@ public:
   std::vector<Vector<Scalar>> gs;
   Vector<Scalar> z0;
   Vector<Scalar> z1;
-  double φ;
 
   Cord(const pSpinModel<Scalar, ps...>& M, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : M(M), gs(ng, Vector<Scalar>::Zero(z2.size())) {
     Scalar H2 = M.getHamiltonian(z2);
@@ -148,7 +147,7 @@ public:
     for (unsigned i = 0; i < gs.size(); i++) {
       Real fdgn = gCoeff(i, t);
       Real dfdgn = dgCoeff(i, t);
-      Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * (
+      Vector<Scalar> dC = - 1 / ż.norm() / dz.norm() / 2 * (
           dfdgn * ż.conjugate() + fdgn * dżc * dz + fdgn * dż * dz.conjugate()
           - std::real(dz.dot(ż)) * (
               dfdgn * dz.conjugate() / dz.squaredNorm() +
@@ -163,7 +162,7 @@ public:
         Real fdgm = gCoeff(j, t);
         Real dfdgm = dgCoeff(j, t);
         Scalar CC = - std::real(dz.dot(ż)) / ż.norm() / dz.norm();
-        Vector<Scalar> dCm = - 0.5 / ż.norm() / dz.norm() * (
+        Vector<Scalar> dCm = - 1 / ż.norm() / dz.norm() / 2 * (
             dfdgm * ż.conjugate() + fdgm * dżc * dz + fdgm * dż * dz.conjugate()
             - std::real(dz.dot(ż)) * (
                 dfdgm * dz.conjugate() / dz.squaredNorm() +
@@ -171,7 +170,7 @@ public:
               )
           );
 
-        Matrix<Scalar> ddC = 0.5 / ż.norm() / dz.norm() * (
+        Matrix<Scalar> ddC = 1 / ż.norm() / dz.norm() / 2 * (
           -  (
                 dfdgn * fdgm * dżc.transpose() + dfdgm * fdgn * dżc +
                 fdgn * fdgm * ddżcdz + fdgn * fdgm * ddżdz
@@ -187,18 +186,18 @@ public:
           - CC * dz.norm() / ż.norm() * fdgn * fdgm * (
                   ddżżc + ddżcż + dżc * dż.transpose() + dż * dżc.transpose()
                 )
-          - 0.5 * CC / dz.norm() / ż.norm() * (
+          - CC / dz.norm() / ż.norm() / (Real)2 * (
             dfdgn * fdgm * dz.conjugate() * (dżc * ż + dż * ż.conjugate()).transpose() +
             dfdgm * fdgn * (dżc * ż + dż * ż.conjugate()) * dz.adjoint()
             )
-          + 0.5 * CC * ż.norm() / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.adjoint()
-          + 0.5 * CC * dz.norm() / pow(ż.norm(), 3)
+          + CC * ż.norm() / (Real)2 / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.adjoint()
+          + CC * dz.norm() / (Real)2 / pow(ż.norm(), 3)
                 * fdgn * (dżc * ż + dż * ż.conjugate())
                 * fdgm * (dżc * ż + dż * ż.conjugate()).transpose()
               )
             ;
 
-        Matrix<Scalar> dcdC = 0.5 / ż.norm() / dz.norm() * (
+        Matrix<Scalar> dcdC = 1 / ż.norm() / dz.norm() / 2 * (
           - (
                 dfdgn * fdgm * dż.adjoint() + dfdgm * fdgn * dż +
                 fdgn * fdgm * (dcdżcdz + dcdżcdz.adjoint())
@@ -215,12 +214,12 @@ public:
           - CC * dz.norm() / ż.norm() * fdgn * fdgm * (
                   dcdżcż + dcdżcż.adjoint() + dżc * dżc.adjoint() + dż * dż.adjoint()
                 )
-          - 0.5 * CC / dz.norm() / ż.norm() * (
+          - CC / dz.norm() / ż.norm() / (Real)2 * (
             dfdgn * fdgm * dz.conjugate() * (dżc * ż + dż * ż.conjugate()).adjoint() +
             dfdgm * fdgn * (dżc * ż + dż * ż.conjugate()) * dz.transpose()
             )
-          + 0.5 * CC * ż.norm() / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.transpose()
-          + 0.5 * CC * dz.norm() / pow(ż.norm(), 3) * fdgn * fdgm
+          + CC * ż.norm() / (Real)2 / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.transpose()
+          + CC * dz.norm() / (Real)2 / pow(ż.norm(), 3) * fdgn * fdgm
              * (dżc * ż + dż * ż.conjugate())
              * (dżc * ż + dż * ż.conjugate()).adjoint()
              )
-- 
cgit v1.2.3-70-g09d2