adrn/Make-all-Hessians.ipynb

## Make-all-Hessians.ipynb
{
  "cells": [
    {
      "metadata": {
        "trusted": true,
        "ExecuteTime": {
          "start_time": "2020-10-14T20:20:14.614257Z",
          "end_time": "2020-10-14T20:20:14.667993Z"
        }
      },
      "cell_type": "code",
      "source": "import astropy.coordinates as coord\nimport astropy.units as u\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport numpy as np\nfrom scipy.special import gamma, gammainc\n\n# gala\nimport gala.coordinates as gc\nimport gala.dynamics as gd\nimport gala.potential as gp\nfrom gala.units import galactic\n\nimport sympy as sy",
      "execution_count": 5,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:52:45.754312Z",
          "start_time": "2020-10-14T18:52:45.749368Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "potential_names = []\nfor name in dir(gp):\n    pot = getattr(gp, name)\n    \n    try:\n        issubclass(pot, object)\n    except TypeError:\n        continue\n        \n    if (not issubclass(pot, gp.CPotentialBase) or \n            issubclass(pot, dict) or\n            pot == gp.CPotentialBase or \n            pot == gp.PotentialBase):\n        continue\n    \n    print(name)\n    potential_names.append(name)",
      "execution_count": 2,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "HenonHeilesPotential\nHernquistPotential\nIsochronePotential\nJaffePotential\nKeplerPotential\nLeeSutoTriaxialNFWPotential\nLogarithmicPotential\nLongMuraliBarPotential\nMiyamotoNagaiPotential\nNFWPotential\nNullPotential\nPlummerPotential\nPowerLawCutoffPotential\nSCFPotential\nSatohPotential\nStonePotential\n"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:30:44.367966Z",
          "start_time": "2020-10-14T18:30:42.578492Z"
        },
        "scrolled": true,
        "trusted": true
      },
      "cell_type": "code",
      "source": "for name in potential_names:\n    pot = getattr(gp, name)\n    \n    try:\n        Phi, v, p = pot.to_sympy()\n    except NotImplementedError:\n        continue\n        \n    Hess = sy.hessian(Phi, list(v.values()))\n\n    CSE_results = sy.cse(Hess, sy.numbered_symbols(\"tmp_\"))\n    \n    print(f\"-- {name} --\")\n    for helper in CSE_results[0]:\n        print(\"double \" + sy.ccode(helper[1], helper[0]))\n        \n    tmp = sy.MatrixSymbol('hess', pot.ndim**2, 1)\n    print(sy.ccode(CSE_results[1][0].reshape(pot.ndim**2, 1), assign_to=tmp))\n    \n    print(\"\\n\\n\")",
      "execution_count": 3,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "-- HenonHeilesPotential --\ndouble tmp_0 = 2.0*y;\ndouble tmp_1 = 2.0*x;\nhess[0] = tmp_0 + 1.0;\nhess[1] = tmp_1;\nhess[2] = tmp_1;\nhess[3] = 1.0 - tmp_0;\n\n\n\n-- HernquistPotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = tmp_0 + tmp_1 + tmp_2;\ndouble tmp_4 = sqrt(tmp_3);\ndouble tmp_5 = c + tmp_4;\ndouble tmp_6 = G*m;\ndouble tmp_7 = tmp_6/pow(tmp_5, 2);\ndouble tmp_8 = tmp_7/tmp_4;\ndouble tmp_9 = 2*tmp_6/(tmp_3*pow(tmp_5, 3));\ndouble tmp_10 = tmp_7/pow(tmp_3, 3.0/2.0);\ndouble tmp_11 = tmp_9*x;\ndouble tmp_12 = tmp_10*x;\ndouble tmp_13 = -tmp_11*y - tmp_12*y;\ndouble tmp_14 = -tmp_11*z - tmp_12*z;\ndouble tmp_15 = y*z;\ndouble tmp_16 = -tmp_10*tmp_15 - tmp_15*tmp_9;\nhess[0] = -tmp_0*tmp_10 - tmp_0*tmp_9 + tmp_8;\nhess[1] = tmp_13;\nhess[2] = tmp_14;\nhess[3] = tmp_13;\nhess[4] = -tmp_1*tmp_10 - tmp_1*tmp_9 + tmp_8;\nhess[5] = tmp_16;\nhess[6] = tmp_14;\nhess[7] = tmp_16;\nhess[8] = -tmp_10*tmp_2 - tmp_2*tmp_9 + tmp_8;\n\n\n\n-- IsochronePotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = pow(b, 2) + tmp_0 + tmp_1 + tmp_2;\ndouble tmp_4 = sqrt(tmp_3);\ndouble tmp_5 = b + tmp_4;\ndouble tmp_6 = G*m;\ndouble tmp_7 = tmp_6/pow(tmp_5, 2);\ndouble tmp_8 = tmp_7/tmp_4;\ndouble tmp_9 = 2*tmp_6/(tmp_3*pow(tmp_5, 3));\ndouble tmp_10 = tmp_7/pow(tmp_3, 3.0/2.0);\ndouble tmp_11 = tmp_9*x;\ndouble tmp_12 = tmp_10*x;\ndouble tmp_13 = -tmp_11*y - tmp_12*y;\ndouble tmp_14 = -tmp_11*z - tmp_12*z;\ndouble tmp_15 = y*z;\ndouble tmp_16 = -tmp_10*tmp_15 - tmp_15*tmp_9;\nhess[0] = -tmp_0*tmp_10 - tmp_0*tmp_9 + tmp_8;\nhess[1] = tmp_13;\nhess[2] = tmp_14;\nhess[3] = tmp_13;\nhess[4] = -tmp_1*tmp_10 - tmp_1*tmp_9 + tmp_8;\nhess[5] = tmp_16;\nhess[6] = tmp_14;\nhess[7] = tmp_16;\nhess[8] = -tmp_10*tmp_2 - tmp_2*tmp_9 + tmp_8;\n\n\n\n-- JaffePotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = tmp_0 + tmp_1 + tmp_2;\ndouble tmp_4 = 1.0/tmp_3;\ndouble tmp_5 = sqrt(tmp_3);\ndouble tmp_6 = c + tmp_5;\ndouble tmp_7 = pow(tmp_6, -2);\ndouble tmp_8 = tmp_7*x;\ndouble tmp_9 = 1.0/tmp_5;\ndouble tmp_10 = 1.0/tmp_6;\ndouble tmp_11 = tmp_10*tmp_9;\ndouble tmp_12 = G*m/c;\ndouble tmp_13 = tmp_12*(tmp_11*x - tmp_8);\ndouble tmp_14 = tmp_13*tmp_4;\ndouble tmp_15 = pow(tmp_3, -3.0/2.0);\ndouble tmp_16 = tmp_13*tmp_15*tmp_6;\ndouble tmp_17 = tmp_10*tmp_15;\ndouble tmp_18 = tmp_4*tmp_7;\ndouble tmp_19 = 2*tmp_9/pow(tmp_6, 3);\ndouble tmp_20 = tmp_11 - tmp_7;\ndouble tmp_21 = tmp_12*tmp_6;\ndouble tmp_22 = tmp_21*tmp_9;\ndouble tmp_23 = tmp_19*x;\ndouble tmp_24 = tmp_4*tmp_8;\ndouble tmp_25 = tmp_17*x;\ndouble tmp_26 = tmp_14*y - tmp_16*y + tmp_22*(tmp_23*y - tmp_24*y - tmp_25*y);\ndouble tmp_27 = tmp_4*z;\ndouble tmp_28 = tmp_13*tmp_27 - tmp_16*z + tmp_22*(tmp_23*z - tmp_24*z - tmp_25*z);\ndouble tmp_29 = tmp_7*y;\ndouble tmp_30 = tmp_11*y - tmp_29;\ndouble tmp_31 = tmp_12*tmp_30;\ndouble tmp_32 = tmp_15*tmp_21;\ndouble tmp_33 = tmp_30*tmp_32;\ndouble tmp_34 = y*z;\ndouble tmp_35 = tmp_22*(-tmp_17*tmp_34 + tmp_19*tmp_34 - tmp_27*tmp_29) + tmp_27*tmp_31 - tmp_33*z;\ndouble tmp_36 = tmp_11*z - tmp_7*z;\nhess[0] = tmp_14*x - tmp_16*x + tmp_22*(-tmp_0*tmp_17 - tmp_0*tmp_18 + tmp_0*tmp_19 + tmp_20);\nhess[1] = tmp_26;\nhess[2] = tmp_28;\nhess[3] = tmp_26;\nhess[4] = tmp_22*(-tmp_1*tmp_17 - tmp_1*tmp_18 + tmp_1*tmp_19 + tmp_20) + tmp_31*tmp_4*y - tmp_33*y;\nhess[5] = tmp_35;\nhess[6] = tmp_28;\nhess[7] = tmp_35;\nhess[8] = tmp_12*tmp_27*tmp_36 + tmp_22*(-tmp_17*tmp_2 - tmp_18*tmp_2 + tmp_19*tmp_2 + tmp_20) - tmp_32*tmp_36*z;\n\n\n\n-- KeplerPotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = tmp_0 + tmp_1 + tmp_2;\ndouble tmp_4 = G*m;\ndouble tmp_5 = tmp_4/pow(tmp_3, 3.0/2.0);\ndouble tmp_6 = 3*tmp_4/pow(tmp_3, 5.0/2.0);\ndouble tmp_7 = tmp_6*x;\ndouble tmp_8 = -tmp_7*y;\ndouble tmp_9 = -tmp_7*z;\ndouble tmp_10 = -tmp_6*y*z;\nhess[0] = -tmp_0*tmp_6 + tmp_5;\nhess[1] = tmp_8;\nhess[2] = tmp_9;\nhess[3] = tmp_8;\nhess[4] = -tmp_1*tmp_6 + tmp_5;\nhess[5] = tmp_10;\nhess[6] = tmp_9;\nhess[7] = tmp_10;\nhess[8] = -tmp_2*tmp_6 + tmp_5;\n\n\n\n-- LogarithmicPotential --\ndouble tmp_0 = pow(q1, -2);\ndouble tmp_1 = pow(v_c, 2);\ndouble tmp_2 = tmp_0*tmp_1;\ndouble tmp_3 = pow(x, 2);\ndouble tmp_4 = pow(q2, -2);\ndouble tmp_5 = pow(y, 2);\ndouble tmp_6 = pow(q3, -2);\ndouble tmp_7 = pow(z, 2);\ndouble tmp_8 = pow(r_h, 2) + tmp_0*tmp_3 + tmp_4*tmp_5 + tmp_6*tmp_7;\ndouble tmp_9 = 1.0/tmp_8;\ndouble tmp_10 = 2.0/pow(tmp_8, 2);\ndouble tmp_11 = tmp_1*tmp_10;\ndouble tmp_12 = tmp_4*y;\ndouble tmp_13 = tmp_10*tmp_2*x;\ndouble tmp_14 = -tmp_12*tmp_13;\ndouble tmp_15 = tmp_6*z;\ndouble tmp_16 = -tmp_13*tmp_15;\ndouble tmp_17 = tmp_1*tmp_9;\ndouble tmp_18 = -tmp_11*tmp_12*tmp_15;\nhess[0] = tmp_2*tmp_9 - tmp_11*tmp_3/pow(q1, 4);\nhess[1] = tmp_14;\nhess[2] = tmp_16;\nhess[3] = tmp_14;\nhess[4] = tmp_17*tmp_4 - tmp_11*tmp_5/pow(q2, 4);\nhess[5] = tmp_18;\nhess[6] = tmp_16;\nhess[7] = tmp_18;\nhess[8] = tmp_17*tmp_6 - tmp_11*tmp_7/pow(q3, 4);\n\n\n\n-- LongMuraliBarPotential --\ndouble tmp_0 = cos(alpha);\ndouble tmp_1 = tmp_0*x;\ndouble tmp_2 = sin(alpha);\ndouble tmp_3 = tmp_2*y;\ndouble tmp_4 = tmp_1 + tmp_3;\ndouble tmp_5 = a + tmp_4;\ndouble tmp_6 = tmp_0*tmp_5;\ndouble tmp_7 = tmp_0*y - tmp_2*x;\ndouble tmp_8 = tmp_2*tmp_7;\ndouble tmp_9 = -tmp_8;\ndouble tmp_10 = tmp_6 + tmp_9;\ndouble tmp_11 = pow(z, 2);\ndouble tmp_12 = pow(c, 2) + tmp_11;\ndouble tmp_13 = sqrt(tmp_12);\ndouble tmp_14 = b + tmp_13;\ndouble tmp_15 = pow(tmp_14, 2);\ndouble tmp_16 = tmp_15 + pow(tmp_7, 2);\ndouble tmp_17 = tmp_16 + pow(tmp_5, 2);\ndouble tmp_18 = sqrt(tmp_17);\ndouble tmp_19 = 1.0/tmp_18;\ndouble tmp_20 = tmp_10*tmp_19;\ndouble tmp_21 = tmp_18 + tmp_5;\ndouble tmp_22 = 1.0/tmp_21;\ndouble tmp_23 = a - tmp_1 - tmp_3;\ndouble tmp_24 = tmp_0*tmp_23;\ndouble tmp_25 = -tmp_24 + tmp_9;\ndouble tmp_26 = tmp_16 + pow(tmp_23, 2);\ndouble tmp_27 = sqrt(tmp_26);\ndouble tmp_28 = 1.0/tmp_27;\ndouble tmp_29 = tmp_25*tmp_28;\ndouble tmp_30 = tmp_0 + tmp_29;\ndouble tmp_31 = -tmp_0;\ndouble tmp_32 = -tmp_20 + tmp_31;\ndouble tmp_33 = pow(tmp_21, -2);\ndouble tmp_34 = -a + tmp_27 + tmp_4;\ndouble tmp_35 = tmp_33*tmp_34;\ndouble tmp_36 = tmp_22*tmp_30 + tmp_32*tmp_35;\ndouble tmp_37 = (1.0/2.0)*G*m/a;\ndouble tmp_38 = tmp_37/tmp_34;\ndouble tmp_39 = tmp_36*tmp_38;\ndouble tmp_40 = tmp_21*tmp_37/pow(tmp_34, 2);\ndouble tmp_41 = tmp_36*tmp_40;\ndouble tmp_42 = tmp_32*tmp_33;\ndouble tmp_43 = pow(tmp_0, 2) + pow(tmp_2, 2);\ndouble tmp_44 = tmp_28*tmp_43;\ndouble tmp_45 = pow(tmp_26, -3.0/2.0);\ndouble tmp_46 = tmp_25*tmp_45;\ndouble tmp_47 = -tmp_19*tmp_43;\ndouble tmp_48 = pow(tmp_17, -3.0/2.0);\ndouble tmp_49 = tmp_10*tmp_48;\ndouble tmp_50 = tmp_34/pow(tmp_21, 3);\ndouble tmp_51 = tmp_32*tmp_50;\ndouble tmp_52 = tmp_21*tmp_38;\ndouble tmp_53 = tmp_0*tmp_7;\ndouble tmp_54 = tmp_2*tmp_5;\ndouble tmp_55 = tmp_53 + tmp_54;\ndouble tmp_56 = tmp_19*tmp_55;\ndouble tmp_57 = tmp_2 + tmp_56;\ndouble tmp_58 = -tmp_2;\ndouble tmp_59 = tmp_2*tmp_23;\ndouble tmp_60 = tmp_53 - tmp_59;\ndouble tmp_61 = tmp_28*tmp_60;\ndouble tmp_62 = tmp_58 - tmp_61;\ndouble tmp_63 = -tmp_53;\ndouble tmp_64 = tmp_59 + tmp_63;\ndouble tmp_65 = tmp_22*tmp_46;\ndouble tmp_66 = -tmp_56 + tmp_58;\ndouble tmp_67 = tmp_33*tmp_66;\ndouble tmp_68 = tmp_2 + tmp_61;\ndouble tmp_69 = -tmp_54 + tmp_63;\ndouble tmp_70 = tmp_35*tmp_49;\ndouble tmp_71 = -2*tmp_2 - 2*tmp_56;\ndouble tmp_72 = tmp_39*tmp_57 + tmp_41*tmp_62 + tmp_52*(tmp_30*tmp_67 + tmp_42*tmp_68 + tmp_51*tmp_71 + tmp_64*tmp_65 - tmp_69*tmp_70);\ndouble tmp_73 = 1.0/tmp_13;\ndouble tmp_74 = tmp_14*tmp_73;\ndouble tmp_75 = tmp_74*z;\ndouble tmp_76 = tmp_19*tmp_75;\ndouble tmp_77 = tmp_28*tmp_75;\ndouble tmp_78 = tmp_19*tmp_33;\ndouble tmp_79 = tmp_75*tmp_78;\ndouble tmp_80 = 2*tmp_76;\ndouble tmp_81 = tmp_39*tmp_76 - tmp_41*tmp_77 + tmp_52*(-tmp_30*tmp_79 + tmp_42*tmp_77 - tmp_51*tmp_80 - tmp_65*tmp_75 + tmp_70*tmp_75);\ndouble tmp_82 = tmp_22*tmp_68 + tmp_35*tmp_66;\ndouble tmp_83 = tmp_38*tmp_82;\ndouble tmp_84 = tmp_40*tmp_82;\ndouble tmp_85 = tmp_45*tmp_60;\ndouble tmp_86 = tmp_50*tmp_66;\ndouble tmp_87 = tmp_48*tmp_55;\ndouble tmp_88 = tmp_52*(-tmp_22*tmp_75*tmp_85 + tmp_35*tmp_75*tmp_87 + tmp_67*tmp_77 - tmp_68*tmp_79 - tmp_80*tmp_86) + tmp_76*tmp_83 - tmp_77*tmp_84;\ndouble tmp_89 = tmp_22*tmp_28;\ndouble tmp_90 = tmp_14*tmp_89;\ndouble tmp_91 = tmp_73*tmp_90;\ndouble tmp_92 = tmp_19*tmp_35;\ndouble tmp_93 = tmp_74*tmp_92;\ndouble tmp_94 = tmp_91*z - tmp_93*z;\ndouble tmp_95 = tmp_11/tmp_12;\ndouble tmp_96 = tmp_11/pow(tmp_12, 3.0/2.0);\ndouble tmp_97 = tmp_15*tmp_95;\ndouble tmp_98 = 2*tmp_97;\nhess[0] = tmp_39*(tmp_0 + tmp_20) + tmp_41*(-tmp_29 + tmp_31) + tmp_52*(tmp_22*(tmp_44 + tmp_46*(tmp_24 + tmp_8)) + 2*tmp_30*tmp_42 + tmp_35*(tmp_47 - tmp_49*(-tmp_6 + tmp_8)) + tmp_51*(-2*tmp_0 - 2*tmp_20));\nhess[1] = tmp_72;\nhess[2] = tmp_81;\nhess[3] = tmp_72;\nhess[4] = tmp_52*(tmp_22*(tmp_44 + tmp_64*tmp_85) + tmp_35*(tmp_47 - tmp_69*tmp_87) + 2*tmp_67*tmp_68 + tmp_71*tmp_86) + tmp_57*tmp_83 + tmp_62*tmp_84;\nhess[5] = tmp_88;\nhess[6] = tmp_81;\nhess[7] = tmp_88;\nhess[8] = tmp_38*tmp_76*tmp_94 - tmp_40*tmp_77*tmp_94 + tmp_52*(tmp_14*tmp_92*tmp_96 - tmp_22*tmp_45*tmp_97 - tmp_28*tmp_78*tmp_98 + tmp_35*tmp_48*tmp_97 + tmp_89*tmp_95 - tmp_90*tmp_96 + tmp_91 - tmp_92*tmp_95 - tmp_93 + tmp_50*tmp_98/tmp_17);\n\n\n\n-- MiyamotoNagaiPotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = pow(b, 2) + tmp_2;\ndouble tmp_4 = sqrt(tmp_3);\ndouble tmp_5 = a + tmp_4;\ndouble tmp_6 = pow(tmp_5, 2);\ndouble tmp_7 = tmp_0 + tmp_1 + tmp_6;\ndouble tmp_8 = G*m;\ndouble tmp_9 = tmp_8/pow(tmp_7, 3.0/2.0);\ndouble tmp_10 = 3*tmp_8/pow(tmp_7, 5.0/2.0);\ndouble tmp_11 = tmp_10*x;\ndouble tmp_12 = -tmp_11*y;\ndouble tmp_13 = tmp_5/tmp_4;\ndouble tmp_14 = tmp_13*z;\ndouble tmp_15 = -tmp_11*tmp_14;\ndouble tmp_16 = -tmp_10*tmp_14*y;\ndouble tmp_17 = 1.0/tmp_3;\ndouble tmp_18 = tmp_2*tmp_9;\nhess[0] = -tmp_0*tmp_10 + tmp_9;\nhess[1] = tmp_12;\nhess[2] = tmp_15;\nhess[3] = tmp_12;\nhess[4] = -tmp_1*tmp_10 + tmp_9;\nhess[5] = tmp_16;\nhess[6] = tmp_15;\nhess[7] = tmp_16;\nhess[8] = -tmp_10*tmp_17*tmp_2*tmp_6 + tmp_13*tmp_9 + tmp_17*tmp_18 - tmp_18*tmp_5/pow(tmp_3, 3.0/2.0);\n\n\n\n"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "-- NFWPotential --\ndouble tmp_0 = pow(a, -2);\ndouble tmp_1 = G*m;\ndouble tmp_2 = tmp_0*tmp_1;\ndouble tmp_3 = pow(x, 2);\ndouble tmp_4 = pow(b, -2);\ndouble tmp_5 = pow(y, 2);\ndouble tmp_6 = pow(c, -2);\ndouble tmp_7 = pow(z, 2);\ndouble tmp_8 = tmp_0*tmp_3 + tmp_4*tmp_5 + tmp_6*tmp_7;\ndouble tmp_9 = pow(tmp_8, -3.0/2.0);\ndouble tmp_10 = 1.0/r_s;\ndouble tmp_11 = tmp_10*sqrt(tmp_8) + 1;\ndouble tmp_12 = log(tmp_11);\ndouble tmp_13 = tmp_12*tmp_9;\ndouble tmp_14 = tmp_3/pow(a, 4);\ndouble tmp_15 = 3*tmp_1;\ndouble tmp_16 = tmp_12/pow(tmp_8, 5.0/2.0);\ndouble tmp_17 = tmp_15*tmp_16;\ndouble tmp_18 = tmp_10/tmp_11;\ndouble tmp_19 = tmp_18/tmp_8;\ndouble tmp_20 = tmp_9/(pow(r_s, 2)*pow(tmp_11, 2));\ndouble tmp_21 = tmp_1*tmp_20;\ndouble tmp_22 = tmp_18/pow(tmp_8, 2);\ndouble tmp_23 = tmp_15*tmp_22;\ndouble tmp_24 = tmp_4*y;\ndouble tmp_25 = tmp_2*x;\ndouble tmp_26 = 3*tmp_25;\ndouble tmp_27 = tmp_16*tmp_26;\ndouble tmp_28 = tmp_20*tmp_25;\ndouble tmp_29 = tmp_22*tmp_26;\ndouble tmp_30 = -tmp_24*tmp_27 + tmp_24*tmp_28 + tmp_24*tmp_29;\ndouble tmp_31 = tmp_6*z;\ndouble tmp_32 = -tmp_27*tmp_31 + tmp_28*tmp_31 + tmp_29*tmp_31;\ndouble tmp_33 = tmp_1*tmp_13;\ndouble tmp_34 = tmp_5/pow(b, 4);\ndouble tmp_35 = tmp_1*tmp_19;\ndouble tmp_36 = tmp_24*tmp_31;\ndouble tmp_37 = -tmp_17*tmp_36 + tmp_21*tmp_36 + tmp_23*tmp_36;\ndouble tmp_38 = tmp_7/pow(c, 4);\nhess[0] = tmp_13*tmp_2 - tmp_14*tmp_17 + tmp_14*tmp_21 + tmp_14*tmp_23 - tmp_19*tmp_2;\nhess[1] = tmp_30;\nhess[2] = tmp_32;\nhess[3] = tmp_30;\nhess[4] = -tmp_17*tmp_34 + tmp_21*tmp_34 + tmp_23*tmp_34 + tmp_33*tmp_4 - tmp_35*tmp_4;\nhess[5] = tmp_37;\nhess[6] = tmp_32;\nhess[7] = tmp_37;\nhess[8] = -tmp_17*tmp_38 + tmp_21*tmp_38 + tmp_23*tmp_38 + tmp_33*tmp_6 - tmp_35*tmp_6;\n\n\n\n-- PlummerPotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = pow(b, 2) + tmp_0 + tmp_1 + tmp_2;\ndouble tmp_4 = G*m;\ndouble tmp_5 = tmp_4/pow(tmp_3, 3.0/2.0);\ndouble tmp_6 = 3*tmp_4/pow(tmp_3, 5.0/2.0);\ndouble tmp_7 = tmp_6*x;\ndouble tmp_8 = -tmp_7*y;\ndouble tmp_9 = -tmp_7*z;\ndouble tmp_10 = -tmp_6*y*z;\nhess[0] = -tmp_0*tmp_6 + tmp_5;\nhess[1] = tmp_8;\nhess[2] = tmp_9;\nhess[3] = tmp_8;\nhess[4] = -tmp_1*tmp_6 + tmp_5;\nhess[5] = tmp_10;\nhess[6] = tmp_9;\nhess[7] = tmp_10;\nhess[8] = -tmp_2*tmp_6 + tmp_5;\n\n\n\n-- PowerLawCutoffPotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = tmp_0 + tmp_1 + tmp_2;\ndouble tmp_4 = (1.0/2.0)*alpha;\ndouble tmp_5 = -tmp_4;\ndouble tmp_6 = tmp_5 + 1.5;\ndouble tmp_7 = pow(r_c, -2);\ndouble tmp_8 = tmp_3*tmp_7;\ndouble tmp_9 = G*m;\ndouble tmp_10 = tmp_9/tgamma(tmp_5 + 2.5);\ndouble // Not supported in C:\n// lowergamma\ntmp_11 = tmp_10*lowergamma(tmp_6, tmp_8);\ndouble tmp_12 = tmp_11/pow(tmp_3, 5.0/2.0);\ndouble tmp_13 = (9.0/2.0)*tmp_12;\ndouble tmp_14 = exp(-tmp_8);\ndouble tmp_15 = tmp_0*tmp_14;\ndouble tmp_16 = pow(tmp_8, -tmp_4)*tmp_9/tgamma(tmp_6);\ndouble tmp_17 = 4*tmp_16/pow(r_c, 5);\ndouble tmp_18 = alpha*tmp_0;\ndouble tmp_19 = (3.0/2.0)*tmp_12;\ndouble tmp_20 = 6*tmp_15;\ndouble tmp_21 = pow(r_c, -4);\ndouble tmp_22 = tmp_5 + 0.5;\ndouble tmp_23 = pow(tmp_8, tmp_22);\ndouble tmp_24 = tmp_10*tmp_23/sqrt(tmp_3);\ndouble tmp_25 = tmp_21*tmp_24;\ndouble tmp_26 = pow(tmp_3, -3.0/2.0);\ndouble tmp_27 = tmp_10*tmp_23*tmp_26*tmp_7;\ndouble tmp_28 = tmp_20*tmp_27;\ndouble tmp_29 = 2*tmp_14;\ndouble tmp_30 = tmp_18*tmp_29;\ndouble tmp_31 = tmp_16*tmp_29/pow(r_c, 3);\ndouble tmp_32 = tmp_31/tmp_3;\ndouble tmp_33 = tmp_27*tmp_30;\ndouble tmp_34 = tmp_11*tmp_26;\ndouble tmp_35 = tmp_14*tmp_24;\ndouble tmp_36 = tmp_35*tmp_7;\ndouble tmp_37 = alpha*tmp_36 + tmp_31 - tmp_34*tmp_4 + (3.0/2.0)*tmp_34 - 3*tmp_36;\ndouble tmp_38 = tmp_13*x;\ndouble tmp_39 = alpha*tmp_19;\ndouble tmp_40 = x*y;\ndouble tmp_41 = tmp_14*tmp_17;\ndouble tmp_42 = alpha*tmp_40;\ndouble tmp_43 = tmp_21*tmp_35;\ndouble tmp_44 = 6*tmp_40;\ndouble tmp_45 = tmp_14*tmp_27;\ndouble tmp_46 = tmp_44*tmp_45;\ndouble tmp_47 = tmp_29*tmp_42;\ndouble tmp_48 = tmp_27*tmp_47;\ndouble tmp_49 = -tmp_22*tmp_46 + tmp_22*tmp_48 - tmp_25*tmp_47 - tmp_32*tmp_42 - tmp_38*y + tmp_39*tmp_40 - tmp_40*tmp_41 + tmp_43*tmp_44 + tmp_46 - tmp_48;\ndouble tmp_50 = x*z;\ndouble tmp_51 = alpha*tmp_50;\ndouble tmp_52 = 6*tmp_50;\ndouble tmp_53 = tmp_45*tmp_52;\ndouble tmp_54 = tmp_29*tmp_51;\ndouble tmp_55 = tmp_27*tmp_54;\ndouble tmp_56 = -tmp_22*tmp_53 + tmp_22*tmp_55 - tmp_25*tmp_54 - tmp_32*tmp_51 - tmp_38*z + tmp_39*tmp_50 - tmp_41*tmp_50 + tmp_43*tmp_52 + tmp_53 - tmp_55;\ndouble tmp_57 = 6*tmp_1;\ndouble tmp_58 = tmp_45*tmp_57;\ndouble tmp_59 = alpha*tmp_1;\ndouble tmp_60 = tmp_29*tmp_59;\ndouble tmp_61 = tmp_27*tmp_60;\ndouble tmp_62 = y*z;\ndouble tmp_63 = alpha*tmp_62;\ndouble tmp_64 = 6*tmp_62;\ndouble tmp_65 = tmp_45*tmp_64;\ndouble tmp_66 = tmp_29*tmp_63;\ndouble tmp_67 = tmp_27*tmp_66;\ndouble tmp_68 = -tmp_13*tmp_62 - tmp_22*tmp_65 + tmp_22*tmp_67 - tmp_25*tmp_66 - tmp_32*tmp_63 + tmp_39*tmp_62 - tmp_41*tmp_62 + tmp_43*tmp_64 + tmp_65 - tmp_67;\ndouble tmp_69 = 6*tmp_2;\ndouble tmp_70 = tmp_45*tmp_69;\ndouble tmp_71 = alpha*tmp_2;\ndouble tmp_72 = tmp_29*tmp_71;\ndouble tmp_73 = tmp_27*tmp_72;\nhess[0] = -tmp_0*tmp_13 - tmp_15*tmp_17 + tmp_18*tmp_19 - tmp_18*tmp_32 + tmp_20*tmp_25 - tmp_22*tmp_28 + tmp_22*tmp_33 - tmp_25*tmp_30 + tmp_28 - tmp_33 + tmp_37;\nhess[1] = tmp_49;\nhess[2] = tmp_56;\nhess[3] = tmp_49;\nhess[4] = -tmp_1*tmp_13 + tmp_1*tmp_39 - tmp_1*tmp_41 - tmp_22*tmp_58 + tmp_22*tmp_61 - tmp_25*tmp_60 - tmp_32*tmp_59 + tmp_37 + tmp_43*tmp_57 + tmp_58 - tmp_61;\nhess[5] = tmp_68;\nhess[6] = tmp_56;\nhess[7] = tmp_68;\nhess[8] = -tmp_13*tmp_2 + tmp_2*tmp_39 - tmp_2*tmp_41 - tmp_22*tmp_70 + tmp_22*tmp_73 - tmp_25*tmp_72 - tmp_32*tmp_71 + tmp_37 + tmp_43*tmp_69 + tmp_70 - tmp_73;\n\n\n\n-- SatohPotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = pow(b, 2) + tmp_2;\ndouble tmp_4 = sqrt(tmp_3);\ndouble tmp_5 = a*(a + 2*tmp_4) + tmp_0 + tmp_1 + tmp_2;\ndouble tmp_6 = G*m;\ndouble tmp_7 = tmp_6/pow(tmp_5, 3.0/2.0);\ndouble tmp_8 = tmp_6/pow(tmp_5, 5.0/2.0);\ndouble tmp_9 = 3*tmp_8;\ndouble tmp_10 = -tmp_9*x*y;\ndouble tmp_11 = 3*z;\ndouble tmp_12 = a/tmp_4;\ndouble tmp_13 = tmp_8*(-tmp_11*tmp_12 - tmp_11);\ndouble tmp_14 = tmp_13*x;\ndouble tmp_15 = tmp_13*y;\nhess[0] = -tmp_0*tmp_9 + tmp_7;\nhess[1] = tmp_10;\nhess[2] = tmp_14;\nhess[3] = tmp_10;\nhess[4] = -tmp_1*tmp_9 + tmp_7;\nhess[5] = tmp_15;\nhess[6] = tmp_14;\nhess[7] = tmp_15;\nhess[8] = -tmp_13*(-tmp_12*z - z) - tmp_7*(a*tmp_2/pow(tmp_3, 3.0/2.0) - tmp_12 - 1);\n\n\n\n-- StonePotential --\ndouble tmp_0 = pow(r_h, 2);\ndouble tmp_1 = 1.0/tmp_0;\ndouble tmp_2 = pow(x, 2);\ndouble tmp_3 = pow(y, 2);\ndouble tmp_4 = pow(z, 2);\ndouble tmp_5 = tmp_2 + tmp_3 + tmp_4;\ndouble tmp_6 = tmp_1*tmp_5 + 1;\ndouble tmp_7 = 1.0/tmp_6;\ndouble tmp_8 = 3/pow(tmp_5, 2);\ndouble tmp_9 = tmp_2*tmp_8;\ndouble tmp_10 = pow(r_c, 2);\ndouble tmp_11 = 1.0/tmp_10;\ndouble tmp_12 = tmp_11*tmp_5 + 1;\ndouble tmp_13 = 1.0/tmp_12;\ndouble tmp_14 = tmp_10 + tmp_5;\ndouble tmp_15 = pow(tmp_14, -2);\ndouble tmp_16 = 8*tmp_15;\ndouble tmp_17 = tmp_0 + tmp_5;\ndouble tmp_18 = 8*tmp_17/pow(tmp_14, 3);\ndouble tmp_19 = 2/tmp_14;\ndouble tmp_20 = 2*tmp_15*tmp_17;\ndouble tmp_21 = tmp_19 - tmp_20;\ndouble tmp_22 = 1.0/tmp_17;\ndouble tmp_23 = 0.5*tmp_14*tmp_22;\ndouble tmp_24 = tmp_19*x - tmp_20*x;\ndouble tmp_25 = 1.0*tmp_22;\ndouble tmp_26 = tmp_24*tmp_25;\ndouble tmp_27 = sqrt(tmp_5);\ndouble tmp_28 = r_c*atan(tmp_27/r_c);\ndouble tmp_29 = 3/pow(tmp_5, 5.0/2.0);\ndouble tmp_30 = tmp_2*tmp_29;\ndouble tmp_31 = 1.0/tmp_5;\ndouble tmp_32 = 2*tmp_31;\ndouble tmp_33 = tmp_2*tmp_32;\ndouble tmp_34 = tmp_1/pow(tmp_6, 2);\ndouble tmp_35 = tmp_11/pow(tmp_12, 2);\ndouble tmp_36 = r_h*atan(tmp_27/r_h);\ndouble tmp_37 = 1.0*tmp_14/pow(tmp_17, 2);\ndouble tmp_38 = tmp_24*tmp_37;\ndouble tmp_39 = pow(tmp_5, -3.0/2.0);\ndouble tmp_40 = -tmp_13*tmp_31 + tmp_28*tmp_39 + tmp_31*tmp_7 - tmp_36*tmp_39;\ndouble tmp_41 = 2*G*m/(-3.1415926535897931*r_c + 3.1415926535897931*r_h);\ndouble tmp_42 = x*y;\ndouble tmp_43 = tmp_42*tmp_8;\ndouble tmp_44 = tmp_29*tmp_42;\ndouble tmp_45 = tmp_32*tmp_42;\ndouble tmp_46 = tmp_16*x;\ndouble tmp_47 = -tmp_41*(tmp_13*tmp_43 + tmp_23*(tmp_18*tmp_42 - tmp_46*y) + tmp_26*y - tmp_28*tmp_44 - tmp_34*tmp_45 + tmp_35*tmp_45 + tmp_36*tmp_44 - tmp_38*y - tmp_43*tmp_7);\ndouble tmp_48 = x*z;\ndouble tmp_49 = tmp_48*tmp_8;\ndouble tmp_50 = tmp_29*tmp_48;\ndouble tmp_51 = tmp_32*tmp_48;\ndouble tmp_52 = -tmp_41*(tmp_13*tmp_49 + tmp_23*(tmp_18*tmp_48 - tmp_46*z) + tmp_26*z - tmp_28*tmp_50 - tmp_34*tmp_51 + tmp_35*tmp_51 + tmp_36*tmp_50 - tmp_38*z - tmp_49*tmp_7);\ndouble tmp_53 = tmp_3*tmp_8;\ndouble tmp_54 = tmp_19*y - tmp_20*y;\ndouble tmp_55 = tmp_25*tmp_54;\ndouble tmp_56 = tmp_29*tmp_3;\ndouble tmp_57 = tmp_3*tmp_32;\ndouble tmp_58 = tmp_37*tmp_54;\ndouble tmp_59 = y*z;\ndouble tmp_60 = tmp_59*tmp_8;\ndouble tmp_61 = tmp_29*tmp_59;\ndouble tmp_62 = tmp_32*tmp_59;\ndouble tmp_63 = -tmp_41*(tmp_13*tmp_60 + tmp_23*(-tmp_16*tmp_59 + tmp_18*tmp_59) - tmp_28*tmp_61 - tmp_34*tmp_62 + tmp_35*tmp_62 + tmp_36*tmp_61 + tmp_55*z - tmp_58*z - tmp_60*tmp_7);\ndouble tmp_64 = tmp_4*tmp_8;\ndouble tmp_65 = z*(tmp_19*z - tmp_20*z);\ndouble tmp_66 = tmp_29*tmp_4;\ndouble tmp_67 = tmp_32*tmp_4;\nhess[0] = -tmp_41*(tmp_13*tmp_9 + tmp_23*(-tmp_16*tmp_2 + tmp_18*tmp_2 + tmp_21) + tmp_26*x - tmp_28*tmp_30 + tmp_30*tmp_36 - tmp_33*tmp_34 + tmp_33*tmp_35 - tmp_38*x + tmp_40 - tmp_7*tmp_9);\nhess[1] = tmp_47;\nhess[2] = tmp_52;\nhess[3] = tmp_47;\nhess[4] = -tmp_41*(tmp_13*tmp_53 + tmp_23*(-tmp_16*tmp_3 + tmp_18*tmp_3 + tmp_21) - tmp_28*tmp_56 - tmp_34*tmp_57 + tmp_35*tmp_57 + tmp_36*tmp_56 + tmp_40 - tmp_53*tmp_7 + tmp_55*y - tmp_58*y);\nhess[5] = tmp_63;\nhess[6] = tmp_52;\nhess[7] = tmp_63;\nhess[8] = -tmp_41*(tmp_13*tmp_64 + tmp_23*(-tmp_16*tmp_4 + tmp_18*tmp_4 + tmp_21) + tmp_25*tmp_65 - tmp_28*tmp_66 - tmp_34*tmp_67 + tmp_35*tmp_67 + tmp_36*tmp_66 - tmp_37*tmp_65 + tmp_40 - tmp_64*tmp_7);\n\n\n\n"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T13:05:57.608639Z",
          "start_time": "2020-10-14T13:05:49.202161Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "print(\"Flattened NFW\")\npot = gp.NFWPotential\n\nPhi, v, p = pot.to_sympy()\n\nHess = sy.hessian(Phi, (v['x'], v['y'], v['z']))\nHess = Hess.subs({p['a']: 1, p['b']: 1})\nHess.simplify()\n\nCSE_results = sy.cse(Hess, sy.numbered_symbols(\"tmp_\"))\n\nfor helper in CSE_results[0]:\n    print(\"double \" + sy.ccode(helper[1], helper[0]))\n\ntmp = sy.MatrixSymbol('hess', 9, 1)\nprint(sy.ccode(CSE_results[1][0].reshape(9, 1), assign_to=tmp))",
      "execution_count": 14,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "Flattened NFW\n-- StonePotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(z, 2);\ndouble tmp_2 = pow(c, 2);\ndouble tmp_3 = pow(y, 2);\ndouble tmp_4 = tmp_1 + tmp_2*(tmp_0 + tmp_3);\ndouble tmp_5 = pow(tmp_4, 4);\ndouble tmp_6 = 3*tmp_0;\ndouble tmp_7 = tmp_4/tmp_2;\ndouble tmp_8 = sqrt(tmp_7);\ndouble tmp_9 = r_s + tmp_8;\ndouble tmp_10 = pow(c, 8);\ndouble tmp_11 = tmp_10*tmp_9;\ndouble tmp_12 = tmp_11*pow(tmp_7, 7.0/2.0);\ndouble tmp_13 = pow(tmp_4, 3);\ndouble tmp_14 = pow(tmp_9, 2);\ndouble tmp_15 = tmp_14*log(tmp_9/r_s);\ndouble tmp_16 = tmp_15*tmp_2;\ndouble tmp_17 = -tmp_11*pow(tmp_7, 9.0/2.0) + tmp_15*tmp_5;\ndouble tmp_18 = G*m/tmp_14;\ndouble tmp_19 = tmp_18/pow(tmp_7, 11.0/2.0);\ndouble tmp_20 = tmp_19/tmp_10;\ndouble tmp_21 = pow(c, 4);\ndouble tmp_22 = pow(tmp_4, 2);\ndouble tmp_23 = 3*tmp_9;\ndouble tmp_24 = pow(tmp_7, 3.0/2.0);\ndouble tmp_25 = tmp_18*x;\ndouble tmp_26 = tmp_21*tmp_25*y*(-3*tmp_15*tmp_21*tmp_24 + tmp_21*pow(tmp_7, 5.0/2.0) + tmp_22*tmp_23)/tmp_5;\ndouble tmp_27 = 3*tmp_16;\ndouble tmp_28 = tmp_2*z*(tmp_2*tmp_24 + tmp_23*tmp_4 - tmp_27*tmp_8)/tmp_13;\ndouble tmp_29 = tmp_25*tmp_28;\ndouble tmp_30 = tmp_18*tmp_28*y;\nhess[0] = tmp_20*(tmp_0*tmp_5 + tmp_12*tmp_6 - tmp_13*tmp_16*tmp_6 + tmp_17);\nhess[1] = tmp_26;\nhess[2] = tmp_29;\nhess[3] = tmp_26;\nhess[4] = tmp_20*(3*tmp_12*tmp_3 - tmp_13*tmp_27*tmp_3 + tmp_17 + tmp_3*tmp_5);\nhess[5] = tmp_30;\nhess[6] = tmp_29;\nhess[7] = tmp_30;\nhess[8] = tmp_13*tmp_19*(tmp_1*tmp_2*tmp_23*tmp_8 - tmp_1*tmp_27 + tmp_1*tmp_4 + tmp_16*tmp_4 - tmp_22*tmp_9/tmp_8)/pow(c, 12);\n"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T13:08:13.857828Z",
          "start_time": "2020-10-14T13:08:05.520976Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "print(\"Spherical NFW\")\npot = gp.NFWPotential\n\nPhi, v, p = pot.to_sympy()\n\nHess = sy.hessian(Phi, (v['x'], v['y'], v['z']))\nHess = Hess.subs({p['a']: 1, p['b']: 1, p['c']: 1})\nHess.simplify()\n\nCSE_results = sy.cse(Hess, sy.numbered_symbols(\"tmp_\"))\n\nfor helper in CSE_results[0]:\n    print(\"double \" + sy.ccode(helper[1], helper[0]))\n\ntmp = sy.MatrixSymbol('hess', 9, 1)\nprint(sy.ccode(CSE_results[1][0].reshape(9, 1), assign_to=tmp))",
      "execution_count": 15,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "Spherical NFW\n-- StonePotential --\ndouble tmp_0 = pow(x, 2);\ndouble tmp_1 = pow(y, 2);\ndouble tmp_2 = pow(z, 2);\ndouble tmp_3 = tmp_0 + tmp_1 + tmp_2;\ndouble tmp_4 = pow(tmp_3, 7);\ndouble tmp_5 = 3*tmp_0;\ndouble tmp_6 = sqrt(tmp_3);\ndouble tmp_7 = r_s + tmp_6;\ndouble tmp_8 = pow(tmp_3, 13.0/2.0)*tmp_7;\ndouble tmp_9 = pow(tmp_7, 2);\ndouble tmp_10 = 1.0/r_s;\ndouble tmp_11 = tmp_9*log(tmp_10*tmp_7);\ndouble tmp_12 = tmp_11*pow(tmp_3, 6);\ndouble tmp_13 = tmp_11*tmp_4 - pow(tmp_3, 15.0/2.0)*tmp_7;\ndouble tmp_14 = G*m;\ndouble tmp_15 = tmp_14/tmp_9;\ndouble tmp_16 = tmp_15/pow(tmp_3, 17.0/2.0);\ndouble tmp_17 = x*y;\ndouble tmp_18 = 4*tmp_15/pow(tmp_3, 3.0/2.0);\ndouble tmp_19 = 3*tmp_17;\ndouble tmp_20 = r_s*tmp_15/pow(tmp_3, 2);\ndouble tmp_21 = tmp_14*log(tmp_10*tmp_6 + 1)/pow(tmp_3, 5.0/2.0);\ndouble tmp_22 = tmp_17*tmp_18 + tmp_19*tmp_20 - tmp_19*tmp_21;\ndouble tmp_23 = x*z;\ndouble tmp_24 = 3*tmp_20;\ndouble tmp_25 = 3*tmp_21;\ndouble tmp_26 = tmp_18*tmp_23 + tmp_23*tmp_24 - tmp_23*tmp_25;\ndouble tmp_27 = 3*tmp_8;\ndouble tmp_28 = 3*tmp_12;\ndouble tmp_29 = y*z;\ndouble tmp_30 = tmp_18*tmp_29 + tmp_24*tmp_29 - tmp_25*tmp_29;\nhess[0] = tmp_16*(tmp_0*tmp_4 - tmp_12*tmp_5 + tmp_13 + tmp_5*tmp_8);\nhess[1] = tmp_22;\nhess[2] = tmp_26;\nhess[3] = tmp_22;\nhess[4] = tmp_16*(tmp_1*tmp_27 - tmp_1*tmp_28 + tmp_1*tmp_4 + tmp_13);\nhess[5] = tmp_30;\nhess[6] = tmp_26;\nhess[7] = tmp_30;\nhess[8] = tmp_16*(tmp_13 + tmp_2*tmp_27 - tmp_2*tmp_28 + tmp_2*tmp_4);\n"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:52:45.760737Z",
          "start_time": "2020-10-14T18:52:45.757460Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "import mpmath\nfrom astropy.utils.misc import isiterable\ndef lowergamma(a, x):\n    if isiterable(x):\n        return np.array([float(mpmath.gammainc(a, b=x[i]))\n                         for i in range(len(x))])\n    else:\n        return float(mpmath.gammainc(a, b=x))",
      "execution_count": 3,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true,
        "ExecuteTime": {
          "start_time": "2020-10-15T01:45:47.652197Z",
          "end_time": "2020-10-15T01:45:51.087926Z"
        }
      },
      "cell_type": "code",
      "source": "failures = []\nfor name in potential_names[1:]:\n    Potential = getattr(gp, name)\n    \n    try:\n        Phi, v, p = Potential.to_sympy()\n    except NotImplementedError:\n        continue\n        \n    print(name)\n    \n    modules = ['numpy', \n               {'atan': np.arctan, 'lowergamma': lowergamma, 'gamma': gamma}]\n    \n    e_func = sy.lambdify(list(p.values()) + list(v.values()), Phi,\n                         modules=modules)\n\n    grad = sy.derive_by_array(Phi, list(v.values()))\n    grad_func = sy.lambdify(list(p.values()) + list(v.values()), grad,\n                            modules=modules)\n    \n    Hess = sy.hessian(Phi, list(v.values()))\n    Hess_func = sy.lambdify(list(p.values()) + list(v.values()), Hess,\n                            modules=modules)\n                            \n    par_vals = {}\n    for k in Potential._parameters:\n        if k == 'm':\n            par_vals[k] = 1e10    \n        elif k == 'alpha':\n            par_vals[k] = 2.1\n        elif k == 'r_c':\n            par_vals[k] = 3.1\n        elif k == 'r_h':\n            par_vals[k] = 10.1\n        elif k == 'b':\n            par_vals[k] = 0.9\n        elif k == 'c':\n            par_vals[k] = 0.846\n        else:\n            par_vals[k] = 1\n            \n    pot = Potential(**par_vals, units=galactic)\n    \n    N = 64\n    trial_x = np.random.uniform(-3., 3., size=(pot.ndim, N))\n    x_dict = {k: v for k, v in zip(['x', 'y', 'z'], trial_x)}\n    \n    f_sympy = e_func(G=pot.G, **par_vals, **x_dict)\n    G_sympy = grad_func(G=pot.G, **par_vals, **x_dict)\n    H_sympy = Hess_func(G=pot.G, **par_vals, **x_dict)\n    \n    f_gala = pot.energy(trial_x).value\n    e_close = np.allclose(f_gala, f_sympy)\n                            \n    G_gala = pot.gradient(trial_x).value\n    g_close = np.allclose(G_gala, G_sympy)\n                            \n    H_gala = pot.hessian(trial_x).value\n    h_close = np.allclose(H_gala, H_sympy)\n                            \n    print(f'energy {e_close}, gradient {g_close}, hessian {h_close}\\n')\n    if not all([e_close, g_close, h_close]):\n        failures.append(name)",
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "text": "HernquistPotential\nenergy True, gradient True, hessian True\n\nIsochronePotential\nenergy True, gradient True, hessian True\n\nJaffePotential\nenergy True, gradient True, hessian True\n\nKeplerPotential\nenergy True, gradient True, hessian True\n\nLogarithmicPotential\nenergy True, gradient True, hessian True\n\nLongMuraliBarPotential\nenergy True, gradient True, hessian True\n\nMiyamotoNagaiPotential\nenergy True, gradient True, hessian True\n\nNFWPotential\nenergy True, gradient True, hessian True\n\nPlummerPotential\nenergy True, gradient True, hessian True\n\nPowerLawCutoffPotential\nenergy True, gradient True, hessian True\n\nSatohPotential\nenergy True, gradient True, hessian True\n\nStonePotential\nenergy True, gradient True, hessian True\n\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "---\n\n### Gamma function debugging crap below - don't look at this!"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:39:20.851911Z",
          "start_time": "2020-10-14T18:39:20.841147Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "CSE_results = sy.cse(Phi, sy.numbered_symbols(\"tmp_\"))\n\nfor helper in CSE_results[0]:\n    print(\"double \" + sy.ccode(helper[1], helper[0]))\n\n# tmp = sy.MatrixSymbol('hess', 9, 1)\nprint(sy.ccode(CSE_results[1][0]))",
      "execution_count": 21,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "double tmp_0 = (1.0/2.0)*alpha;\ndouble tmp_1 = -tmp_0;\ndouble tmp_2 = tmp_1 + 1.5;\ndouble tmp_3 = pow(x, 2) + pow(y, 2) + pow(z, 2);\ndouble tmp_4 = tmp_3/pow(r_c, 2);\ndouble tmp_5 = G*m;\ndouble // Not supported in C:\n// lowergamma\ntmp_6 = tmp_5*lowergamma(tmp_2, tmp_4)/(sqrt(tmp_3)*tgamma(tmp_1 + 2.5));\n// Not supported in C:\n// lowergamma\ntmp_0*tmp_6 - 3.0/2.0*tmp_6 + tmp_5*lowergamma(tmp_1 + 1, tmp_4)/(r_c*tgamma(tmp_2))\n"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:36:25.053077Z",
          "start_time": "2020-10-14T18:36:24.862132Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "Phi",
      "execution_count": 19,
      "outputs": [
        {
          "data": {
            "text/latex": "$\\displaystyle \\frac{G \\alpha m \\gamma\\left(1.5 - \\frac{\\alpha}{2}, \\frac{x^{2} + y^{2} + z^{2}}{r_{c}^{2}}\\right)}{2 \\sqrt{x^{2} + y^{2} + z^{2}} \\Gamma\\left(2.5 - \\frac{\\alpha}{2}\\right)} - \\frac{3 G m \\gamma\\left(1.5 - \\frac{\\alpha}{2}, \\frac{x^{2} + y^{2} + z^{2}}{r_{c}^{2}}\\right)}{2 \\sqrt{x^{2} + y^{2} + z^{2}} \\Gamma\\left(2.5 - \\frac{\\alpha}{2}\\right)} + \\frac{G m \\gamma\\left(1 - \\frac{\\alpha}{2}, \\frac{x^{2} + y^{2} + z^{2}}{r_{c}^{2}}\\right)}{r_{c} \\Gamma\\left(1.5 - \\frac{\\alpha}{2}\\right)}$",
            "text/plain": "G*alpha*m*lowergamma(1.5 - alpha/2, (x**2 + y**2 + z**2)/r_c**2)/(2*sqrt(x**2 + y**2 + z**2)*gamma(2.5 - alpha/2)) - 3*G*m*lowergamma(1.5 - alpha/2, (x**2 + y**2 + z**2)/r_c**2)/(2*sqrt(x**2 + y**2 + z**2)*gamma(2.5 - alpha/2)) + G*m*lowergamma(1 - alpha/2, (x**2 + y**2 + z**2)/r_c**2)/(r_c*gamma(1.5 - alpha/2))"
          },
          "execution_count": 19,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:30:47.297888Z",
          "start_time": "2020-10-14T18:30:47.291653Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "f_gala",
      "execution_count": 6,
      "outputs": [
        {
          "data": {
            "text/plain": "array([3.15433872, 3.16574385, 3.14378289, 3.32737306, 3.13413157,\n       3.21269569, 3.13500574, 3.18922066, 3.14985041, 3.13272827,\n       3.15920551, 3.27160331, 3.14126244, 3.15781781, 3.1787848 ,\n       3.12818331, 3.16971518, 3.17561444, 3.15850503, 3.14649336,\n       3.16323922, 3.16924558, 3.13322471, 3.13960862, 3.35809196,\n       3.13599589, 3.13262894, 3.16814016, 3.3349993 , 3.13051177,\n       3.13738946, 3.14424797, 3.13611046, 3.13318257, 3.14389834,\n       3.14025543, 3.13956221, 3.14850932, 3.17397065, 3.15284405,\n       3.40466389, 3.13092589, 3.20075893, 3.15184996, 3.29718763,\n       3.27699005, 3.16648554, 3.20870378, 3.37658457, 3.14526664,\n       3.52288265, 3.17528366, 3.12901874, 3.21689608, 3.20503863,\n       3.21581988, 3.14277545, 3.76994585, 3.20710442, 3.14300249,\n       3.19520901, 3.1485141 , 3.15059167, 3.14754267])"
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:30:47.820775Z",
          "start_time": "2020-10-14T18:30:47.816448Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "f_sympy",
      "execution_count": 7,
      "outputs": [
        {
          "data": {
            "text/plain": "array([-0.16721369, -0.1683995 , -0.1659389 , -0.18003474, -0.1644733 ,\n       -0.17236617, -0.16462638, -0.17050795, -0.16669887, -0.16421459,\n       -0.16773807, -0.17649074, -0.16559425, -0.16759168, -0.16961164,\n       -0.16320925, -0.16878184, -0.1693277 , -0.16766446, -0.16628907,\n       -0.16815123, -0.1687373 , -0.16430811, -0.16535608, -0.18190395,\n       -0.16479335, -0.16419559, -0.16863176, -0.18050315, -0.16376429,\n       -0.16501828, -0.16600038, -0.16481227, -0.16430026, -0.16595422,\n       -0.16545048, -0.16534924, -0.16653801, -0.16917801, -0.16704595,\n       -0.18466509, -0.16385309, -0.1714441 , -0.16693242, -0.17814608,\n       -0.17684428, -0.16847192, -0.17206236, -0.18300934, -0.16613296,\n       -0.19143796, -0.16929772, -0.16342107, -0.17268141, -0.17177951,\n       -0.17260105, -0.16580356, -0.20513821, -0.17193941, -0.16583433,\n       -0.17100006, -0.16653859, -0.1667863 , -0.16641977])"
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:26:41.323833Z",
          "start_time": "2020-10-14T18:26:41.318429Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "def safe_lowergamma(a, x):\n    if a < 0:\n        N = int(np.ceil(-a))\n        A = np.prod([a + n for n in range(N)])\n        \n        B = np.sum([x**(a+n) * np.exp(-x) * np.prod([a+m for m in range(N-1, n, -1)])\n                    for n in range(N)], axis=0)\n        \n        return (B + safe_lowergamma(a + N, x)) / A\n    \n    else:\n        return lowergamma(a, x)",
      "execution_count": 81,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:26:41.466511Z",
          "start_time": "2020-10-14T18:26:41.464590Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "a = -.35\nx = 0.56",
      "execution_count": 82,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:26:41.630408Z",
          "start_time": "2020-10-14T18:26:41.627032Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "lowergamma(a, x)",
      "execution_count": 83,
      "outputs": [
        {
          "data": {
            "text/plain": "-4.451220157231834"
          },
          "execution_count": 83,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:26:42.109712Z",
          "start_time": "2020-10-14T18:26:42.105838Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "safe_lowergamma(a, x)",
      "execution_count": 84,
      "outputs": [
        {
          "data": {
            "text/plain": "-4.451220157231835"
          },
          "execution_count": 84,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:26:42.673857Z",
          "start_time": "2020-10-14T18:26:42.669980Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "%load_ext cython",
      "execution_count": 85,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "The cython extension is already loaded. To reload it, use:\n  %reload_ext cython\n"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:26:44.185314Z",
          "start_time": "2020-10-14T18:26:43.160950Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "%%cython -lgsl -lgslcblas\n\nfrom libc.math cimport exp\nfrom cython_gsl.gsl_gamma cimport *\n\ncpdef double safe_gamma_inc(double a, double x):\n    cdef:\n        int N, m, n;\n        double A = 1.;\n        double B = 0.;\n        double tmp;\n\n    if (a > 0):\n        return gsl_sf_gamma_inc_P(a, x) * gsl_sf_gamma(a);\n    \n    else:\n        N = <int> ceil(-a);\n\n        for n in range(0, N):\n            A = A * (a + n);\n\n            tmp = 1.;\n            for m in range(N-1, n, -1):\n                tmp = tmp * (a + m)\n            B = B + pow(x, a+n) * exp(-x) * tmp;        \n        return (B + gsl_sf_gamma_inc_P(a + N, x) * gsl_sf_gamma(a + N)) / A ;",
      "execution_count": 86,
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": "warning: /Users/apricewhelan/anaconda/lib/python3.8/site-packages/cython_gsl/gsl_integration.pxd:65:9: 'GSL_EMAXITER' redeclared \nwarning: /Users/apricewhelan/anaconda/lib/python3.8/site-packages/cython_gsl/gsl_integration.pxd:67:9: 'GSL_EROUND' redeclared \nwarning: /Users/apricewhelan/anaconda/lib/python3.8/site-packages/cython_gsl/gsl_integration.pxd:69:9: 'GSL_ESING' redeclared \nwarning: /Users/apricewhelan/anaconda/lib/python3.8/site-packages/cython_gsl/gsl_integration.pxd:71:9: 'GSL_EDIVERGE' redeclared \n"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T18:26:44.195306Z",
          "start_time": "2020-10-14T18:26:44.189107Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "safe_gamma_inc(a, x)",
      "execution_count": 87,
      "outputs": [
        {
          "data": {
            "text/plain": "-4.451220157231832"
          },
          "execution_count": 87,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T15:51:51.540744Z",
          "start_time": "2020-10-14T15:51:51.535944Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "f_sympy",
      "execution_count": 8,
      "outputs": [
        {
          "data": {
            "text/plain": "array([-0.16437906, -0.1716586 , -0.16952406, -0.16372658, -0.16403998,\n       -0.17038593, -0.19075267, -0.16406653, -0.16408614, -0.16164308,\n       -0.16786076, -0.16922065, -0.16393867, -0.16687255, -0.16985756,\n       -0.16626539, -0.165896  , -0.16487296, -0.16888773, -0.16505674,\n       -0.16571324, -0.17453405, -0.16343362, -0.17049447, -0.16424199,\n       -0.16772346, -0.17234982, -0.17198436, -0.16654504, -0.16487241,\n       -0.16856596, -0.16758286, -0.16949195, -0.16369496, -0.17472793,\n       -0.17179556, -0.17275124, -0.16253627, -0.16387925, -0.16828466,\n       -0.17555244, -0.16619143, -0.16663711, -0.16563687, -0.16679925,\n       -0.17209439, -0.17492289, -0.16786929, -0.16512356, -0.16395431,\n       -0.16861984, -0.17312716, -0.16479444, -0.16773551, -0.17265879,\n       -0.17171899, -0.16492944, -0.16481624, -0.16833961, -0.16381814,\n       -0.16304819, -0.16936516, -0.1656083 , -0.1655758 ])"
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T15:40:21.041846Z",
          "start_time": "2020-10-14T15:40:21.038834Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "pot",
      "execution_count": 18,
      "outputs": [
        {
          "data": {
            "text/plain": "<PowerLawCutoffPotential: m=1.00e+10, alpha=2.10, r_c=3.10 (kpc,Myr,solMass,rad)>"
          },
          "execution_count": 18,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T15:41:10.647863Z",
          "start_time": "2020-10-14T15:41:10.633171Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "Phi.subs({'G': pot.G, 'm': 1e10, 'y': 0, 'z': 0, 'alpha': 2.1, 'r_c': 3.1})",
      "execution_count": 20,
      "outputs": [
        {
          "data": {
            "text/latex": "$\\displaystyle 0.00737311563282076 \\gamma\\left(-0.05, 0.104058272632674 x^{2}\\right) - \\frac{0.0228566584617443 \\gamma\\left(0.45, 0.104058272632674 x^{2}\\right)}{\\sqrt{x^{2}}}$",
            "text/plain": "0.00737311563282076*lowergamma(-0.05, 0.104058272632674*x**2) - 0.0228566584617443*lowergamma(0.45, 0.104058272632674*x**2)/sqrt(x**2)"
          },
          "execution_count": 20,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T15:42:28.413372Z",
          "start_time": "2020-10-14T15:42:28.409989Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "lowergamma(-0.05, 0.10405 * 3**2), lowergamma(0.45, 0.10405)",
      "execution_count": 25,
      "outputs": [
        {
          "data": {
            "text/plain": "(nan, 0.7775462921976992)"
          },
          "execution_count": 25,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T15:39:35.920693Z",
          "start_time": "2020-10-14T15:39:35.917404Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "f_sympy",
      "execution_count": 16,
      "outputs": [
        {
          "data": {
            "text/plain": "array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])"
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T15:37:59.582477Z",
          "start_time": "2020-10-14T15:37:59.578701Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "pot",
      "execution_count": 5,
      "outputs": [
        {
          "data": {
            "text/plain": "<PowerLawCutoffPotential: m=1.00e+10, alpha=2.10, r_c=0.50 (kpc,Myr,solMass,rad)>"
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "pot.value([15., 0, 1])",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "end_time": "2020-10-14T15:30:59.970496Z",
          "start_time": "2020-10-14T15:30:59.966841Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "failures",
      "execution_count": 5,
      "outputs": [
        {
          "data": {
            "text/plain": "[]"
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "conda-root-py",
      "display_name": "Python [conda env:root] *",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.8.3",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "toc": {
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "base_numbering": 1,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": false,
      "toc_position": {},
      "toc_section_display": true,
      "toc_window_display": false
    },
    "gist": {
      "id": "",
      "data": {
        "description": "projects/gala-notebooks/Make-all-Hessians.ipynb",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}