From 2617c5ddf8a2b362e1ee8ebe8686622e9f38580c Mon Sep 17 00:00:00 2001
From: HarukiST <seto.haruki@qst.go.jp>
Date: Fri, 4 Oct 2024 11:42:15 +0900
Subject: [PATCH] add in bxkappa2 and bxkappa3

---
 hypnotoad/core/mesh.py | 51 ++++++++++++++++++++++++++++++++++++++----
 1 file changed, 47 insertions(+), 4 deletions(-)

diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py
index f6d79150..98f612bf 100644
--- a/hypnotoad/core/mesh.py
+++ b/hypnotoad/core/mesh.py
@@ -82,7 +82,13 @@ class MeshRegion:
             "curl(b/B)",
             doc="Expression used to calculate curvature operator 'bxcv'",
             value_type=str,
-            allowed=["curl(b/B)", "curl(b/B) with x-y derivatives", "bxkappa"],
+            allowed=[
+                "curl(b/B)",
+                "curl(b/B) with x-y derivatives",
+                "bxkappa",
+                "bxkappa2",
+                "bxkappa3",
+            ],
         ),
         curvature_smoothing=WithMeta(
             None,
@@ -1084,7 +1090,7 @@ def calcMetric(self):
                 - self.Rxy * numpy.abs(self.Bpxy) * self.I * self.tanBeta / self.hy
             )
 
-            self.J = self.hy / numpy.abs(self.Bpxy)
+            self.J = self.bpsign * self.hy / self.Bpxy
 
             self.g_11 = (
                 1.0 / (self.Rxy * self.Bpxy * self.cosBeta) ** 2
@@ -1359,10 +1365,10 @@ def curl_bOverB_zetahat(R, Z):
             self.bxcvx = self.Bxy / 2.0 * self.curl_bOverB_x
             self.bxcvy = self.Bxy / 2.0 * self.curl_bOverB_y
             self.bxcvz = self.Bxy / 2.0 * self.curl_bOverB_z
-        elif self.user_options.curvature_type == "bxkappa":
 
-            # bxkappa for cocos1-extension by H. Seto (QST)
+        elif self.user_options.curvature_type == "bxkappa":
 
+            # b0 x kappa terms for cocos1-extension by H. Seto (QST)
             curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy")
             curlxb0w = self.DDX(
                 "#Bxy*#hy/#Bpxy"
@@ -1379,6 +1385,43 @@ def curl_bOverB_zetahat(R, Z):
             self.bxcvy = bxcvv
             self.bxcvz = bxcvw - self.I * bxcvu
 
+        elif self.user_options.curvature_type == "bxkappa2":
+
+            curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy")
+            curlxb0w = self.bpsign * (
+                self.bpsign * self.Bpxy**2 / self.Bxy / self.J * self.DDX("#hy/#Bpxy")
+                + self.DDX("#Bxy")
+                - self.Btxy * self.Rxy / self.Bxy * self.DDX("#Btxy/#Rxy")
+            )
+
+            coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2)
+            bxcvu = curlxb0u
+            bxcvv = curlxb0w * coefv
+            bxcvw = curlxb0w
+
+            self.bxcvx = bxcvu
+            self.bxcvy = bxcvv
+            self.bxcvz = bxcvw - self.I * bxcvu
+
+        elif self.user_options.curvature_type == "bxkappa3":
+
+            pprime = self.meshParent.equilibrium.regions[
+                self.equilibriumRegion.name
+            ].pprime(self.psixy)
+
+            mu0 = 4e-7 * numpy.pi
+            curlxb0u = -self.Btxy * self.Rxy / (self.J * self.Bxy**2) * self.DDY("#Bxy")
+            curlxb0w = -self.bpsign * (mu0 / self.Bxy * pprime + self.DDX("#Bxy"))
+
+            coefv = -self.Btxy * self.Bpxy * self.Rxy / (self.hy * self.Bxy**2)
+            bxcvu = curlxb0u
+            bxcvv = curlxb0w * coefv
+            bxcvw = curlxb0w
+
+            self.bxcvx = bxcvu
+            self.bxcvy = bxcvv
+            self.bxcvz = bxcvw - self.I * bxcvu
+
         else:
             raise ValueError(
                 "Unrecognized option '"