1 files changed, 64 insertions, 0 deletions
diff --git a/cartoon.c b/cartoon.c
new file mode 100644
index 0000000..13e1b71
--- /dev/null
+++ b/cartoon.c
@@ -0,0 +1,64 @@
+    double xx = x[index];
+    int origin = fabs(xx) < 1e-8;
+
+    PDstandardNth2trK    = 0;
+    PDstandardNth2alp  = 0;
+    PDstandardNth22alp = origin ? PDstandardNth11alp : PDstandardNth1alp / xx;
+    PDstandardNth12alp = 0;
+    PDstandardNth23alp = 0;
+
+    PDstandardNth2beta1  = 0;
+    PDstandardNth22beta1 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
+    PDstandardNth12beta1 = 0;
+    PDstandardNth23beta1 = 0;
+    PDstandardNth2beta2  = origin ? PDstandardNth1beta1 : beta1[index] / xx;
+    PDstandardNth22beta2 = 0;
+    PDstandardNth12beta2 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
+    PDstandardNth23beta2 = origin ? PDstandardNth13beta1 : PDstandardNth3beta1 / xx;
+    PDstandardNth2beta3  = 0;
+    PDstandardNth22beta3 = origin ? PDstandardNth11beta3 : PDstandardNth1beta3 / xx;
+    PDstandardNth12beta3 = 0;
+    PDstandardNth23beta3 = 0;
+
+    PDstandardNth2gt11  = 0;
+    PDstandardNth22gt11 = origin ? PDstandardNth11gt22 : PDstandardNth1gt11 / xx - 2 * (gt11[index] - gt22[index]) / (xx * xx);
+    PDstandardNth12gt11 = 0;
+    PDstandardNth23gt11 = 0;
+    PDstandardNth2gt12  = origin ? PDstandardNth1gt11 - PDstandardNth1gt22 : (gt11[index] - gt22[index]) / xx;
+    PDstandardNth22gt12 = 0;
+    PDstandardNth12gt12 = origin ? (PDstandardNth11gt11 - PDstandardNth11gt22) / 2 : (PDstandardNth1gt11 - PDstandardNth1gt22) / xx - (gt11[index] - gt22[index]) / (xx * xx);
+    PDstandardNth23gt12 = origin ? PDstandardNth13gt11 - PDstandardNth13gt22 : (PDstandardNth3gt11  - PDstandardNth3gt22) / xx;
+    PDstandardNth2gt13  = 0;
+    PDstandardNth22gt13 = origin ? PDstandardNth11gt13 / 2 : PDstandardNth1gt13 / xx - gt13[index] / (xx * xx);
+    PDstandardNth12gt13 = 0;
+    PDstandardNth23gt13 = 0;
+    PDstandardNth2gt22  = 0;
+    PDstandardNth22gt22 = origin ? PDstandardNth11gt11 : PDstandardNth1gt22 / xx + 2 * (gt11[index] - gt22[index]) / (xx * xx);
+    PDstandardNth12gt22 = 0;
+    PDstandardNth23gt22 = 0;
+    PDstandardNth2gt23  = origin ? PDstandardNth1gt13 : gt13[index] / xx;
+    PDstandardNth22gt23 = 0;
+    PDstandardNth12gt23 = origin ? PDstandardNth11gt13 / 2 : PDstandardNth1gt13 / xx - gt13[index] / (xx * xx);
+    PDstandardNth23gt23 = origin ? PDstandardNth13gt13 : PDstandardNth3gt13 / xx;
+    PDstandardNth2gt33  = 0;
+    PDstandardNth22gt33 = origin ? PDstandardNth11gt33 : PDstandardNth1gt33 / xx;
+    PDstandardNth12gt33 = 0;
+    PDstandardNth23gt33 = 0;
+
+    PDstandardNth2phi  = 0;
+    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
+    PDstandardNth12phi = 0;
+    PDstandardNth23phi = 0;
+    PDstandardNth2trK  = 0;
+    PDstandardNth2Xt1  = 0;
+    PDstandardNth2Xt2  = origin ? PDstandardNth1Xt1 : Xt1[index] / xx;
+    PDstandardNth2Xt3  = 0;
+
+    PDstandardNth2At11 = origin ? -2 * PDstandardNth1At12 : -2 * At12[index] / xx;
+    PDstandardNth2At12 = origin ? PDstandardNth1At11 - PDstandardNth1At22 : (At11[index] - At22[index]) / xx;
+    PDstandardNth2At13 = origin ? -PDstandardNth1At23 : - At23[index] / xx;
+    PDstandardNth2At22 = origin ? 2 * PDstandardNth1At12 : 2 * At12[index] / xx;
+    PDstandardNth2At23 = origin ? PDstandardNth1At13 : At13[index] / xx;
+    PDstandardNth2At33 = 0;
+    if (fabs(y[index]) > 1e-8)
+        continue;