Improved documentation

Examples.ipynb

@@ -59,8 +59,8 @@
       "Number of points in the n=1 grid  6250000\n",
       "Number of points in the n=2 grid  6250000\n",
       "File  ./Results/LensingBands_a_0.94_i_17.h5  created.\n",
-      "CPU times: user 814 ms, sys: 274 ms, total: 1.09 s\n",
-      "Wall time: 52 s\n"
+      "CPU times: user 929 ms, sys: 324 ms, total: 1.25 s\n",
+      "Wall time: 56.1 s\n"
      ]
     }
    ],

README.md

@@ -6,15 +6,15 @@ The code, described in detail in Ref. [1], implements all the relevant equations
 
 The use of AART in scientific publications must be properly acknowledged. Please cite:
 
-	Cardenas-Avendano, A., Zhu, H. & Lupsasca, A. 
-	Adaptive Analytical Ray-Tracing of Black Hole Photon Rings. 
-	[arXiv:Please come back soon for the reference!]
+_______
+Cardenas-Avendano, A., Lupsasca, A. & Zhu, H. "Adaptive Analytical Ray Tracing of Black Hole Photon Rings." [arXiv:2211.07469](https://arxiv.org/abs/2211.07469) 
+_______
 
 We also request that AART modifications or extensions leading to a scientific publication be made public as free software. 
 
 <center> <em>Feel free to use images and movies produced with this code (with attribution) for your next presentation! </em> </center>
 
-Last updated: 07.03.2022
+Last updated: 11.15.2022
 
 ## AART's Components ##
 
@@ -37,11 +37,15 @@ All the dependencies are located in the <em>init.py</em> file. Most of the libra
 
 To install any missing packages, run
 
-  <code> >> pip install "package_name" </code>
+<code> pip install "package_name" </code>
 
-  or, if using anaconda, search for the missing packages and run, e.g. for h5py (Read and write HDF5 files from Python,) 
+or, if using anaconda, search for the missing packages and run, e.g. for h5py (Read and write HDF5 files from Python,) 
 
-  <code> conda install -c anaconda h5py</code>
+<code> conda install -c anaconda h5py</code>
+
+Sometimes scipy does not update automatically to the latest version. If that is the case, you may want to type 
+
+<code> pip install -U scipy</code>
 
 
 ## How to run AART ##
@@ -56,7 +60,7 @@ We present some examples in the notebook:
 
 The lensing bands are computed by simply running
 
-  <code> >> python lensingbands.py </code>
+  <code> python lensingbands.py </code>
 
 The result will be stored in a HDF5 file that contains the values of the Bardeen's coordinates withing each lensing band. The datasets inside the resulting file are:
 
@@ -77,7 +81,7 @@ This image is produced in the example code:
 
 To compute the equitorial radius, angle, and emission time of a photon, we perform a backward ray-tracing from the observer plane. By running the following, we evaluate the source radius, angle, and time within the grid from each lensing bands:
 
-  <code> >> python raytracing.py </code>
+  <code> python raytracing.py </code>
 
 The result will be stored in a HDF5 file that contains source radius, angle, time, as well as the radial component of the four momentum at the equitorial plane, for lensing bands n=0,1,2. The datasets inside the resulting file are:
 
@@ -96,7 +100,7 @@ This image is produced in the example code:
 
 Once the lensing bands and the rays have been computed, an image can be produced using a defined analytical profile by simply running
 
-  <code> >> python radialintensity.py </code>
+  <code> python radialintensity.py </code>
 
 The datasets inside the resulting file are:
 
@@ -114,11 +118,11 @@ As the dataset produced after ray tracing contains all the information of the BL
 
 One can also use a precomputed equatorial profile. AART currently implements profiles computed with inoisy. The example includes a test case (<em>inoisy.h5</em>), for which one can simply run by 
 
-  <code> >> python iImages.py </code>
+  <code> python iImages.py </code>
 
   or 
 
-  <code> >> python iMovies.py </code>
+  <code> python iMovies.py </code>
 
 to produce images or a set of images, respectively. Images can be produced by using a single equatorial profile, i.e., in the mode "stationary," or using the entire history of the equatorial structure, i.e, in the mode "dynamical." When movies are made, the dynamical version is assumed. In both cases, the resulting datasets inside the resulting file are:
 
@@ -132,7 +136,7 @@ This gif is produced in the example code:
 
 With the images created using radial intensity prifiles, one may then calculate the visibility of the image projected onto a baseline. This function first performs radon transforms of the image at a set of specified angles (radonangles in <em>params.py</em>), and then compute the visibility amplitude by 
 
-  <code> >> python visamp.py </code>
+  <code> python visamp.py </code>
 
 This function creates a set of h5 files, one for each basline angle. These files contains the visibility amplitude as well as the frequency (baseline length in G$\lambda$). The resulting datasets inside the resulting file are:
 
@@ -152,7 +156,7 @@ This image is produced in the example code:
 
 The linear polarization of a given configuration of the magnetic field can be computed by
 
-  <code> >> python polarization.py </code>
+  <code> python polarization.py </code>
 
   The resulting datasets inside the resulting file are:
 
@@ -179,11 +183,11 @@ The linear polarization of a given configuration of the magnetic field can be co
 
 ## References ##
 
-[1] Cardenas-Avendano, A., Zhu, H. & Lupsasca, A. Adaptive Analytical Ray-Tracing of Black Hole Photon Rings. 
+[1] Cardenas-Avendano, A., Lupsasca, A. & Zhu, H. Adaptive Analytical Ray Tracing of Black Hole Photon Rings. [arXiv:2211.07469](https://arxiv.org/abs/2211.07469)
 
-[2] Gralla, S. E., & Lupsasca, A. (2020). Lensing by Kerr black holes. Physical Review D, 101(4), 044031.(arXiv:1910.12873)
+[2] Gralla, S. E., & Lupsasca, A. (2020). Lensing by Kerr black holes. Physical Review D, 101(4), 044031.
 
-[3] Gralla, S. E., & Lupsasca, A. (2020). Null geodesics of the Kerr exterior. Physical Review D, 101(4), 044032.(arXiv:1910.12881)
+[3] Gralla, S. E., & Lupsasca, A. (2020). Null geodesics of the Kerr exterior. Physical Review D, 101(4), 044032.
 
 [4] Gralla, S. E., Lupsasca, A., & Marrone, D. P. (2020). The shape of the black hole photon ring: A precise test of strong-field general relativity. Physical Review D, 102(12), 124004.
 

aart_func/intensity_f.py

@@ -1,17 +1,89 @@
 from aart_func import *
 from params import * 
 
-def gDisk(r,a,lamb):
+def Delta(r,a):
     """
-    Calculates the redshift factor for a photon outside the inner-most stable circular orbit(isco) (assume circular orbit)
+    Calculates the Kerr metric function \Delta(t)
     :param r: radius of the source
     :param a: spin of the black hole
-    :param lamb: angular momentum
+    """
+    return r**2-2*r+a**2
 
-    :return: the redshift factor associated with the ray
+def PIF(r,a):
+    """
+    Calculates PI(r) (Eq. B6 P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    """
+    return (r**2+a**2)**2-a**2*Delta(r,a)
+
+def urbar(r,a):
+    """
+    Calculates the r (contravariant) component of the four velocity for radial infall
+    (Eq. B34b P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    """
+    return -np.sqrt(2*r*(r**2+a**2))/(r**2)
+
+def Omegabar(r,a):
+    """
+    Calculates the angular velocity of the radial infall
+    (Eq. B32a P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    """
+    return (2*a*r)/PIF(r,a)
+
+def Omegahat(r,a,laux):
+    """
+    Calculates the angular velocity of the sub-Keplerian orbit
+    (Eq. B39 P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    """
+    return (a+(1-2/r)*(laux-a))/(PIF(r,a)/(r**2)-(2*a*laux)/r)
+
+def uttilde(r, a,urT,OT):
+    """
+    Calculates the t (contravariant) component of the general four velocity
+    (Eq. B52 P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    :param urT: r (contravariant) component of the general four velocity
+    :param OT: Angular velocity of the general four velocity
+    """
+    return np.sqrt((1 + urT**2*r**2/Delta(r,a))/(1-(r**2+a**2)*OT**2-(2/r)*(1-a*OT)**2))
+
+def Ehat(r,a,laux):
+    """
+    Calculates the orbital energy of the sub-Keplerian flow
+    (Eq. B44a P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    :param laux: sub-Keplerian specific angular momentum
+    """
+    return np.sqrt(Delta(r,a)/(PIF(r,a)/(r**2)-(4*a*laux)/r-(1-2/r)*laux**2))
+
+def nuhat(r,a,laux,Ehataux):
+    """
+    Calculates the radial velocity of the sub-Keplerian flow
+    (Eq. B45 P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    :param laux: sub-Keplerian specific angular momentum
+    :param Ehataux: sub-Keplerian orbital energy
+    """
+    return r/Delta(r,a)*np.sqrt(np.abs(PIF(r,a)/(r**2)-(4*a*laux)/r-(1-2/r)*laux**2-Delta(r,a)/(Ehataux**2)))
+
+def lhat(r,a):
     """
-    #Eqns (P4 B7)
-    return sqrt(r**3 - 3*r**2 + 2*a*r**(3/2))/(r**(3/2) - (lamb- a))
+    Calculates the rspecific angular momentum of the sub-Keplerian flow
+    (Eq. B44b P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    """
+    return sub_kep*(r**2+a**2-2*a*np.sqrt(r))/(np.sqrt(r)*(r-2)+a)
 
 def Rint(r,a,lamb,eta):
     """
@@ -26,75 +98,96 @@ def Rint(r,a,lamb,eta):
     #Eqns (P2 5)
     return (r**2 + a**2 - a*lamb)**2 - (r**2 - 2*r + a**2)*(eta + (lamb - a)**2)
 
-def gGas(r,b,a,lamb,eta):
+def gDisk(r,a,b,lamb,eta):
+    """
+    Calculates the redshift factor for a photon outside the inner-most stable circular orbit(isco) (assume circular orbit)
+    (Eq. B13 P1)
+    :param r: radius of the source
+    :param a: spin of the black hole
+    :param lamb: angular momentum
+    :param eta: Carter constant
+
+    :return: the redshift factor associated with the ray
+    """
+
+    OH=Omegahat(r,a,lhat(r,a))
+    OT=OH+(1-betaphi)*(Omegabar(r,a)-OH)
+    ur=(1-betar)*urbar(r,a)
+    ut=uttilde(r,a,ur,OT)
+    uphi=ut*OT
+
+    return 1/(ut*(1-b*np.sign(ur)*sqrt(np.abs(Rint(r,a,lamb,eta)*ur**2))/Delta(r,a)/ut-lamb*uphi/ut))
+
+def gGas(r,a,b,lamb,eta):
     """
     Calculates the redshift factor for a photon inside the isco (assume infalling orbit)
+    (Eq. B13 P1)
     :param r: radius of the source
-    :param b: sign for the redshift
     :param a: spin of the black hole
+    :param b: sign for the redshift
     :param lamb: angular momentum
     :param eta: carter constant
 
     :return: the redshift factor associated with the ray
     """
-    #calculate radius of the inner-most stable circular orbit
+    #Calculate radius of the inner-most stable circular orbit
     isco=rms(a)
 
-    #Eqns (P2 2)
-    Delta=r**2 - 2*r + a**2
+    lms=lhat(isco,a)
+    OH=Omegahat(r,a,lms)
+    OT=OH+(1-betaphi)*(Omegabar(r,a)-OH)
 
-    #Eqns (P3 B13)
-    lambe=((isco**2 - 2*a*sqrt(isco) + a**2))/(isco**(3/2) - 2*sqrt(isco) + a)
-    #Eqns (P3 B12)
-    H=(2*r - a*lambe)/Delta
+    Ems=Ehat(isco,a,lms)
+    urhat=-Delta(r,a)/(r**2)*nuhat(r, a, lms ,Ems)*Ems
+    ur=urhat+(1-betar)*(urbar(r,a)-urhat)
+    ut=uttilde(r,a,ur,OT)
+    uphi=OT*ut
 
-    #Eqns (P3 B14)
-    gamma=sqrt(1 - 2/3 *1/isco)
-
-    #Eqns (P3 B9-B11)
-    ut=gamma*(1 + (2)/r *(1 + H))
-    uphi=gamma/r**2*(lambe + a*H)
-    ur=-np.sqrt(2/3/isco)*(isco/r - 1)**(3/2)
-    #Eqns (P3 B15)
-    return 1/(ut - uphi*lamb - ur*Delta**(-1)*b*sqrt(Rint(r,a,lamb,eta)))
+    return 1/(ut*(1-b*np.sign(ur)*sqrt(np.abs(Rint(r,a,lamb,eta)*ur**2))/Delta(r,a)/ut-lamb*uphi/ut))
 
 #calculate the observed brightness for a purely radial profile
 def bright_radial(grid,mask,redshift_sign,a,rs,isco,thetao):
     """
     Calculate the brightness of a rotationally symmetric disk
+    (Eq. 50 P1)
     :param grid: alpha and beta grid on the observer plane on which we evaluate the observables
     :param mask: mask out the lensing band, see lb_f.py for detail
     :param redshift_sign: sign of the redshift
     :param a: black hole spin
     :param rs: source radius
     :param isco: radius of the inner-most stable circular orbit
     :param thetao: observer inclination
-    :hidden param profile: radial profile of the intensity
 
     :return: image of a lensed equitorial source with only radial dependence. 
     """
     alpha = grid[:,0][mask]
     beta = grid[:,1][mask]
+
     rs = rs[mask]
+
     lamb,eta = rt.conserved_quantities(alpha,beta,thetao,a)
+
     brightness = np.zeros(rs.shape[0])
     redshift_sign = redshift_sign[mask]
-    brightness[rs>=isco]= gDisk(rs[rs>=isco],a,lamb[rs>=isco])**gfactor*ilp.profile(rs[rs>=isco],a,gammap,mup,sigmap)
 
-    brightness[rs<isco]= gGas(rs[rs<isco],redshift_sign[rs<isco],a,lamb[rs<isco],eta[rs<isco])**gfactor*ilp.profile(rs[rs<isco],a,gammap,mup,sigmap)
+    brightness[rs>=isco]= gDisk(rs[rs>=isco],a,redshift_sign[rs>=isco],lamb[rs>=isco],eta[rs>=isco])**gfactor*ilp.profile(rs[rs>=isco],a,gammap,mup,sigmap)
+    brightness[rs<isco]= gGas(rs[rs<isco],a,redshift_sign[rs<isco],lamb[rs<isco],eta[rs<isco])**gfactor*ilp.profile(rs[rs<isco],a,gammap,mup,sigmap)
 
     r_p = 1+np.sqrt(1-a**2)
     brightness[rs<=r_p] = 0
 
     I = np.zeros(mask.shape)
     I[mask] = brightness
+
     return(I)
 
+
 #calculate the observed brightness for an arbitrary profile, passed in as the interpolation object
 #but ignoring the time delay due to lensing
 def fast_light(grid,mask,redshift_sign,a,isco,rs,th,interpolation,thetao):
     """
     Calculate the black hole image ignoring the time delay due to lensing or geometric effect
+    (Eq. 116 P1)
     :param grid: alpha and beta grid on the observer plane on which we evaluate the observables
     :param mask: mask out the lensing band, see lb_f.py for detail
     :param redshift_sign: sign of the redshift
@@ -115,12 +208,11 @@ def fast_light(grid,mask,redshift_sign,a,isco,rs,th,interpolation,thetao):
     brightness = np.zeros(rs.shape[0])
     redshift_sign = redshift_sign[mask]
 
-    #Eqs. 60 and 61 in 0706.0622
-    x_aux=np.sqrt(rs**2+a**2)*np.cos(th)
-    y_aux=np.sqrt(rs**2+a**2)*np.sin(th)
+    x_aux=rs*np.cos(th)
+    y_aux=rs*np.sin(th)
 
-    brightness[rs>=isco]= gDisk(rs[rs>=isco],a,lamb[rs>=isco])**gfactor*interpolation(np.vstack([x_aux[rs>=isco],y_aux[rs>=isco]]).T)
-    brightness[rs<isco]= gGas(rs[rs<isco],redshift_sign[rs<isco],a,lamb[rs<isco],eta[rs<isco])**gfactor*interpolation(np.vstack([x_aux[rs<isco],y_aux[rs<isco]]).T)
+    brightness[rs>=isco]= gDisk(rs[rs>=isco],a,redshift_sign[rs>=isco],lamb[rs>=isco],eta[rs>=isco])**gfactor*interpolation(np.vstack([x_aux[rs>=isco],y_aux[rs>=isco]]).T)
+    brightness[rs<isco]= gGas(rs[rs<isco],a,redshift_sign[rs<isco],lamb[rs<isco],eta[rs<isco])**gfactor*interpolation(np.vstack([x_aux[rs<isco],y_aux[rs<isco]]).T)
 
     r_p = 1+np.sqrt(1-a**2)
     brightness[rs<=r_p] = 0
@@ -133,6 +225,8 @@ def fast_light(grid,mask,redshift_sign,a,isco,rs,th,interpolation,thetao):
 def slow_light(grid,mask,redshift_sign,a,isco,rs,th,ts,interpolation,thetao):
     """
     Calculate the black hole image including the time delay due to lensing and geometric effect
+    (Eq. 50 P1)
+
     :param grid: alpha and beta grid on the observer plane on which we evaluate the observables
     :param mask: mask out the lensing band, see lb_f.py for detail
     :param redshift_sign: sign of the redshift
@@ -155,12 +249,11 @@ def slow_light(grid,mask,redshift_sign,a,isco,rs,th,ts,interpolation,thetao):
     brightness = np.zeros(rs.shape[0])
     redshift_sign = redshift_sign[mask]
 
-    #Eqs. 60 and 61 in 0706.0622
-    x_aux=np.sqrt(rs**2+a**2)*np.cos(th)
-    y_aux=np.sqrt(rs**2+a**2)*np.sin(th)
+    x_aux=rs*np.cos(th)
+    y_aux=rs*np.sin(th)
 
-    brightness[rs>=isco]= gDisk(rs[rs>=isco],a,lamb[rs>=isco])**gfactor*interpolation(np.vstack([ts[rs>=isco],x_aux[rs>=isco],y_aux[rs>=isco]]).T)
-    brightness[rs<isco]= gGas(rs[rs<isco],redshift_sign[rs<isco],a,lamb[rs<isco],eta[rs<isco])**gfactor*interpolation(np.vstack([ts[rs<isco],x_aux[rs<isco],y_aux[rs<isco]]).T)
+    brightness[rs>=isco]= gDisk(rs[rs>=isco],a,redshift_sign[rs>=isco],lamb[rs>=isco],eta[rs>=isco])**gfactor*interpolation(np.vstack([ts[rs>=isco],x_aux[rs>=isco],y_aux[rs>=isco]]).T)
+    brightness[rs<isco]= gGas(rs[rs<isco],a,redshift_sign[rs<isco],lamb[rs<isco],eta[rs<isco])**gfactor*interpolation(np.vstack([ts[rs<isco],x_aux[rs<isco],y_aux[rs<isco]]).T)
 
     r_p = 1+np.sqrt(1-a**2)
     brightness[rs<=r_p] = 0

aart_func/misc.py

@@ -13,6 +13,7 @@ def cbrt(x):
 def rms(a):
     '''
     ISCO value
+    (Eq. B16 P1)
     :param a: BH spin
     '''
     Z1=1 + (1 - a**2)**(1/3) *((1 + a)**(1/3) + (1 - a)**(1/3))