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) ~ R
s), 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/N
3. 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
- APPLawD V1: 2006, november (current version)
Input: fixed expansion order N.
Output: potential Ψ.
Accuracy ΔΨ/Ψ on potential values not directly available (accuracy can then be estimated from a second call to ADDPLAWS with order ≥ N+1).
- APPLawD V2: 2007, june (coming soon)
Two options :
- Fixed order (setting the flag SWITCH to .FALSE.)
Input: fixed expansion order N.
Output: potential Ψ and corresponding accuracy ΔΨ/Ψ (estimated from N+1-order expansion).
- Fixed accuracy (setting the flag SWITCH to .TRUE.)
Input: minimum accuracy |ΔΨ/Ψ| on potential values.
Output: Potential Ψ and corresponding order N.
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
- Download the following two files :
- Build an object file with the command (f90 refer to your favorite Fortran
90 compiler) :
f90 -c applawdpack.f90
- Compile and link the driver program with the command :
f90 applawdpack.o driverapplawd.f90 -o driverapplawd.exe
- 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.
