From 09bd9058e00b609d8cee2baa8023710ba74e3deb Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 6 Nov 2021 22:03:34 +0100
Subject: Finished Levenburg-Marquardt for Stokes lines.

---
 p-spin.hpp | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

(limited to 'p-spin.hpp')

diff --git a/p-spin.hpp b/p-spin.hpp
index 50da209..15d2525 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -120,6 +120,29 @@ Tensor<Scalar, 3> ddzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
          + zT.conjugate().contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * (2.0 * dH.dot(z) / pow(z², 3));
 }
 
+template <class Scalar>
+Tensor<Scalar, 3> dcdzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) {
+  Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
+  Tensor<Scalar, 1> dHT = Eigen::TensorMap<Tensor<const Scalar, 1>>(dH.data(), {dH.size()});
+  Tensor<Scalar, 2> ddHT = Eigen::TensorMap<Tensor<const Scalar, 2>>(ddH.data(), {dH.size(), dH.size()});
+
+  Matrix<Scalar> I = Matrix<Real>::Identity(z.size(), z.size());
+  Tensor<Scalar, 2> IT = Eigen::TensorMap<Tensor<const Scalar, 2>>(I.data(), {z.size(), z.size()});
+
+  Eigen::array<Eigen::IndexPair<int>, 0> ei = {};
+
+  Scalar z² = z.squaredNorm();
+
+  return ddHT.conjugate().contract(zT.conjugate(), ei) / z²
+        +IT.contract(dHT.conjugate(), ei).shuffle((std::array<int, 3>){0,2,1}) / z²
+        -zT.contract(dHT.conjugate(), ei).contract(zT.conjugate(), ei) / pow(z², 2)
+        -ddHT.conjugate().contract(zT, ip00).contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) / pow(z², 2)
+        -IT.contract(zT.conjugate(), ei) * (dH.dot(z) / pow(z², 2))
+        -IT.contract(zT.conjugate(), ei).shuffle((std::array<int, 3>){0,2,1}) * (dH.dot(z) / pow(z², 2))
+        +zT.contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * (2.0 * dH.dot(z) / pow(z², 3))
+    ;
+}
+
 template <class Scalar>
 Tensor<Scalar, 3> ddzDotConjugate(const Vector<Scalar>& z, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH, const Tensor<Scalar, 3>& dddH) {
   Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
-- 
cgit v1.2.3-54-g00ecf