APPLawD
Accurate Potentials in Power Law Disks
J.M. Huré (Université Bordeaux 1 & LAB/OASU) — 2006, November

Email : jean-marc.hure@obs.u-bordeaux1.fr
Web page : www.obs.u-bordeaux1.fr/radio/JMHure/intro2applawd.html

Current version: APPLawD V1: 2006, november

  What is it ?

APPLawD (Accurate Disk Potentials for Power Law Surface densities) determines the gravitational potential in the equatorial plane of a flat axially symmetric disk (inside and outside) with finite size and power law surface density profile. Potential values are computed on the basis of the density splitting method (Pierens & Huré, 2005, ApJ, 624, 289), were the residual Poisson kernel is expanded over the modulus of the complete elliptic integral of the first kind. In contrast with classical multipole expansions of potential theory, the residual series converges linearly inside sources, leading to very accurate potential values for low order truncations of the series.

The code is easy-to-use; it is written in Fortran 90 (with no external dependancies).
The code works under variable precision.

Current version: APPLawD V1, 006, november

  Reference

Huré J.-M., Pelat D., Pierens A., 2007,
"Generation of potential/surface density pairs in finite size disks. The case of power law distributions",
Astronomy & Astrophysics, accepted for publication

Abstract. We report a simple method to generate potential/surface density pairs in flat axially symmetric finite size disks. Potential/surface density pairs consist of a ``homogeneous'' pair (a closed form expression) corresponding to a uniform disk, and a ``residual'' pair. This residual component is converted into an infinite series of integrals over the radial extent of the disk. For a certain class of surface density distributions (like power laws of the radius), this series is fully analytical. The extraction of the homogeneous pair is equivalent to a convergence acceleration technique, in a matematical sense. In the case of power law distributions (i.e. surface densities of the form Σ(R) ~ Rs), the convergence rate of the residual series is shown to be cubic inside the source. As a consequence, very accurate potential values are obtained by low order truncation of the series. At zero order, relative errors on potential values do not exceed a few percent typically, and scale with the order N of trunctation as 1/N3. This method is superior to the classical multipole expansion whose very slow convergence is often critical for most practical applications.

→ for more details, see Huré, Pelat & Pierens, 2007, A&A, accepted for publication (A&A paper or arXiv:0706.3616)

  Properties and typical calling sequence

Precision : single, double or quadruple.
Computing time : linear with the expansion order.
Reliability : potential are finite everywhere; no singularities.
Accuracy : computer precision reachable with a few orders inside and outside the disk (more terms are need at the disk edges; see the graphical examples).

Typical F90 calling sequence for APPLawD V1: 2006, november (current version):

RIN=0.1D+0 ! disk inner edge ROUT=1.D+0 ! disk outer edge S=-1.2301D+0 ! power index SIGMAROUT=1.01D+0 ! surface density at the outer edge ORDER=4 ! expansion order R=0.13D+0 ! radius where the potential is requested CALL APPLAWD(RIN,ROUT,SIGMARIN,S,ORDER,R,PSI) PRINT*,"The potential due to this flat disk at",R,"is about",PSI

Typical F90 calling sequence for APPLawD V2: 2007, june (coming soon):

RIN=0.1D+0 ! disk inner edge ROUT=1.D+0 ! disk outer edge S=-1.2301D+0 ! power index SIGMAROUT=1.01D+0 ! surface density at the outer edge R=0.13D+0 ! radius where the potential is requested SWITCH=.FALSE. ORDER=10 ! expansion order ! computes the potential, fixed order ORDER CALL APPLAWD(RIN,ROUT,SIGMAROUT,S,R,ORDER,PSIRERROR,SWITCH,PSI) PRINT*,"The potential due to this flat disk at",R,"is about",PSI PRINT*,"Relative accuracy (estimate) at order",ORDER," is ",PSIRERROR SWITCH=.TRUE. PSIRERROR=1.E-06 ! threshold in the relative error ! computes the potential, fixed relative accuracy less than PSIRERROR CALL APPLAWD(RIN,ROUT,SIGMAROUT,S,R,ORDER,PSIRERROR,SWITCH,PSI) PRINT*,"The potential due to this flat disk at",R,"is about",PSI PRINT*,"Relative accuracy (estimate) at order",ORDER," is ",PSIRERROR

  Copyright


APPLawD is protected according to the terms and conditions of the GNU General Public License.

  Credits

If you use APPLawD, please cite the following reference : Huré J.-M., Pelat D., Pierens A., 2007, A&A, accepted

  Portability and optimization

In principle, APPLawD is fully portable and requires no external libraries.

The code is not optimized. Depending on applications, computing time can be saved by specific recoding.

  Versions

  Graphical examples

  • Here is the plot Ψ(R) corresponding to the execution of the driver program with default parameters, namely: disk inner edge at 0.1, outer edge at 1, power law index -1, and zero-order expansion (i.e. one term only).

  • Graph showing the rise of the averaged relative precision ΔΨ/Ψ on potential values as the truncation order N increases, in a typical case (double precision computations). Precision averaged over the radius.

  • Two movies showing the rise of the relative precision ΔΨ/Ψ (local precision):
    • case with negative power law index s=-1.5 (dense inner disk).
    • case with positive power law index s=+1.5 (dense outer disk).
  • Here (also on the left) is a polar plot showing the potential Ψ versus the power law exponent s (as the polar angle) ranging from -3 to 3. This graphs was used as the code logotype

  Grids of potential

The following data files contains potential values Ψ (2nd column) versus the radius R (1st column) in the range [0,1.25] computed for different power law exponents s. The radius is given in units of the outer edge. Potential values are normalized with respect to the central value Ψ(0) (given in the file header). The disk parameters are : inner and outer edges: 0.1 and 1, surface density at the inner edge : 1. The expansion order N (4th column) has been determined so that the relative precision ΔΨ/Ψ (3th column; decimal log.) is better or equal to 10-6.

Here is the corresponding plot Ψ(R;s).

  Download, install and run APPLawD

  1. Download the following two files :

    V1: 2006, november (current version) V2: 2007, june (coming soon)
    Driver program driverapplawd.f90 driverapplawdV2.f90
    APPLawD package applawdpack.f90 applawdpackV2.f90
    Ouput file (with default parameters) driverapplawd.dat driverapplawdV2.dat

  2. Build an object file with the command (f90 refer to your favorite Fortran 90 compiler) :
    f90 -c applawdpack.f90
  3. Compile and link the driver program with the command :
    f90 applawdpack.o driverapplawd.f90 -o driverapplawd.exe
  4. Execution is launched with the command :
    ./driverapplawd.exe
    With default parameters, the data file driverapplawd.dat should be created after runtime.
    On screen (or default output), you should read this message :
    APPLawD Fortran 90 driver program GNU GPL protected, V1.0, 2006, Nov. Copyright (C) 2006 Jean-Marc Huré ACTUAL PRECISION F90, real kind 8 Maximal precision +2.22044604925031E-16 OUTPUT FILE : driverapplawd.dat DISK EDGES Inner edge +1.00000000000000E-01 Outer edge +1.00000000000000E+00 SURFACE DENSITY AT THE INNER EDGE : +1.00000000000000E+00 POWER LAW INDEX S : -1.00000000000000E+00 EXPANSION ORDER N : 0 COMPUTATIONAL BOX : 000016 points From +0.00000000000000E+00 to +1.50000000000000E+00 R, PSI(R), relative error DPSI/PSI --- +0.00000000000000E+00 -1.44675688248309E+00 +0.00E+00 +1.00000000000000E-01 -1.65933787410814E+00 -3.99E-02 +2.00000000000000E-01 -1.60518267308335E+00 -2.64E-02 +3.00000000000000E-01 -1.44105995746879E+00 -1.90E-02 +4.00000000000000E-01 -1.29610878026641E+00 -1.24E-02 +5.00000000000000E-01 -1.16960021956944E+00 -5.94E-03 +6.00000000000000E-01 -1.05688000363628E+00 +2.95E-04 +7.00000000000000E-01 -9.53872605544343E-01 +6.13E-03 +8.00000000000000E-01 -8.56910372208067E-01 +1.14E-02 +9.00000000000000E-01 -7.61481122030367E-01 +1.58E-02 +1.00000000000000E+00 -6.51323470492979E-01 +1.96E-02 +1.10000000000000E+00 -5.61694399578068E-01 +1.71E-02 +1.20000000000000E+00 -5.03825085254129E-01 +1.48E-02 +1.30000000000000E+00 -4.58809971775975E-01 +1.28E-02 +1.40000000000000E+00 -4.22047609255025E-01 +1.11E-02 +1.50000000000000E+00 -3.91190306939325E-01 +9.78E-03 --- Done.
That's all, in principle.
Feel free to change the input parameters (at you own risk!).

Please, contact me if problems or questions.