Finished Levenburg-Marquardt for Stokes lines.

1 files changed, 23 insertions, 0 deletions

diff --git a/p-spin.hpp b/p-spin.hpp
index 50da209..15d2525 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -121,6 +121,29 @@ Tensor<Scalar, 3> ddzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
 }
 
 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()});
   Tensor<Scalar, 1> dHT = Eigen::TensorMap<Tensor<const Scalar, 1>>(dH.data(), {dH.size()});