Current version: APPLawD V1: 2006, november
The code is easytouse; it is written in Fortran 90 (with no external dependancies).
The code works under variable precision.
Current version: APPLawD V1, 006, november
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 
→ for more details, see Huré, Pelat & Pierens, 2007, A&A, accepted for publication (A&A paper or arXiv:0706.3616)
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.E06 ! 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
The code is not optimized. Depending on applications, computing time can be saved by specific recoding.
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).
Two options :

Here is the corresponding plot Ψ(R;s).
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 
f90 c applawdpack.f90
f90 applawdpack.o driverapplawd.f90 o driverapplawd.exe
With default parameters, the data file driverapplawd.dat should be created after runtime../driverapplawd.exe
APPLawD Fortran 90 driver program GNU GPL protected, V1.0, 2006, Nov. Copyright (C) 2006 JeanMarc Huré ACTUAL PRECISION F90, real kind 8 Maximal precision +2.22044604925031E16 OUTPUT FILE : driverapplawd.dat DISK EDGES Inner edge +1.00000000000000E01 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.00000000000000E01 1.65933787410814E+00 3.99E02 +2.00000000000000E01 1.60518267308335E+00 2.64E02 +3.00000000000000E01 1.44105995746879E+00 1.90E02 +4.00000000000000E01 1.29610878026641E+00 1.24E02 +5.00000000000000E01 1.16960021956944E+00 5.94E03 +6.00000000000000E01 1.05688000363628E+00 +2.95E04 +7.00000000000000E01 9.53872605544343E01 +6.13E03 +8.00000000000000E01 8.56910372208067E01 +1.14E02 +9.00000000000000E01 7.61481122030367E01 +1.58E02 +1.00000000000000E+00 6.51323470492979E01 +1.96E02 +1.10000000000000E+00 5.61694399578068E01 +1.71E02 +1.20000000000000E+00 5.03825085254129E01 +1.48E02 +1.30000000000000E+00 4.58809971775975E01 +1.28E02 +1.40000000000000E+00 4.22047609255025E01 +1.11E02 +1.50000000000000E+00 3.91190306939325E01 +9.78E03  Done.
Please, contact me if problems or questions.