check_ELQ.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RRprime_circ(E,Lz,r,th,a):\n",
    "    Delta = r**2 - 2 *r + a**2\n",
    "    Q = np.cos(th)**2 *(a**2 * (1 - E**2) + Lz**2 /np.sin(th)**2 ) #令Vth=0\n",
    "    R = (E*(r**2 + a**2) - a*Lz)**2 - Delta*(r**2 + (Lz - a*E)**2 + Q) \n",
    "    Rprime= (-(-2 + 2*r))*(((-a)*E + Lz)**2 + Q + r**2) - 2*r*(a**2 - 2*r + r**2) + 4*E*r*((-a)*Lz + E*(a**2 + r**2))\n",
    "    return R,Rprime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r=10; th=np.pi/4; spin=0.9\n",
    "def myfunc(x):\n",
    "    E=x[0]\n",
    "    Lz=x[1]\n",
    "    return RRprime_circ(E,Lz,r,th,spin)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "output=fsolve(myfunc,[0.9,3.0],full_output=True,xtol=1e-20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "E=output[0][0];Lz=output[0][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Q = np.cos(th)**2 *(spin**2 * (1 - E**2) + Lz**2 /np.sin(th)**2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.95312958061646835, 2.497639186615801, 6.2752768275531601)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E,Lz,Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getELQ(e,p,iota,spin):\n",
    "    if(e<1e-5):\n",
    "        r=p\n",
    "        th=np.pi/2-iota\n",
    "        def myfunc(x):\n",
    "            E=x[0]\n",
    "            Lz=x[1]\n",
    "            return RRprime_circ(E,Lz,r,th,spin)\n",
    "        output=fsolve(myfunc,[0.9,3.0],full_output=True,xtol=1e-20)\n",
    "        E=output[0][0];Lz=output[0][1]\n",
    "        Q = np.cos(th)**2 *(spin**2 * (1 - E**2) + Lz**2 /np.sin(th)**2 )\n",
    "        return np.array([E,Lz,Q])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from KRZmetric import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.95312958061646835"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "getELQ(0.0000001,10,np.pi/4,0.9)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "10/0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0006185"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu=0.01\n",
    "-6.1850e-4 *mu*mu*1e4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from KRZmetric import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.016179590241289451, 0.020931698727445294, 0.022037893624304687)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "freq3_dt(0.5,10,np.pi/4,0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7853981633974483"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.pi/4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from KRZmetric import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "spin = 0.500000;current_p = 6.500000;current_e = 0.500000;current_iota = 0.000000;\n",
      "E=0.9485149461;Lz=3.1862999457,Q=0.0000000000;\n"
     ]
    }
   ],
   "source": [
    "e=0.5\n",
    "p=6.5\n",
    "iota= 0  *np.pi/180\n",
    "spin=0.5\n",
    "\n",
    "[E,L,Q]=getELQ(e,p,iota,spin)\n",
    "print('spin = %f;current_p = %f;current_e = %f;current_iota = %f;\\nE=%.10f;Lz=%.10f,Q=%.10f;'%(spin,p,e,iota,E,L,Q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(0,6.8,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.945196560826 3.13046788241\n",
      "krz_d[1]= 0.000000; spin = 0.500000;current_p = 6.000000;current_e = 0.500000;current_iota = 0.000000;\n",
      "E=0.9451965608;Lz=3.1304678824; Q=0.0000000000;\n"
     ]
    }
   ],
   "source": [
    "#赤道面 非Kerr 轨道\n",
    "p=6;\n",
    "ecc=0.5\n",
    "d1=0.0\n",
    "spin=0.5\n",
    "\n",
    "rmax = p / (1 - ecc);\n",
    "rmin = p / (1 + ecc);\n",
    "r0 = rmax;\n",
    "th0=np.pi/2.0\n",
    "invgmax=metric_KRZ_inverse(spin, d1, rmax, th0);\n",
    "invgmin=metric_KRZ_inverse(spin, d1, rmin, th0);\n",
    "\n",
    "\n",
    "EoverL = ((invgmax[3][0] - invgmin[3][0]) + sqrt((invgmax[3][0] - invgmin[3][0]) *(invgmax[3][0] - invgmin[3][0]) - (invgmax[0][0] - invgmin[0][0])*(invgmax[3][3] - invgmin[3][3]))) / (invgmax[0][0] - invgmin[0][0]);\n",
    "Lz = sqrt((invgmax[3][0] - invgmin[3][0]) / (EoverL*EoverL*(invgmin[3][0] * invgmax[0][0] - invgmax[3][0] * invgmin[0][0]) + (invgmin[3][0] * invgmax[3][3] - invgmax[3][0] * invgmin[3][3])));\n",
    "\n",
    "E = EoverL*Lz;\n",
    "iota=0;Q=0;\n",
    "print(E,Lz)\n",
    "print('krz_d[1]= %f; spin = %f;current_p = %f;current_e = %f;current_iota = %f;\\nE=%.10f;Lz=%.10f; Q=%.10f;'%(d1,spin,p,ecc,iota,E,Lz,Q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}