!
! APPLAWDPACK is a Fortran-90 module for APPLAWD.
!
! APPLAWD determines the gravitational potential inside and outside
! a flat axially symmetric disk with finite size and power law surface
! density profile from the splitting method. Potential values are computed
! from a truncated series so that accuracy is easily controlled through the
! expansion order.
! 
! see Huré, Pelat & Pierens (2006), A&A, submitted
!
! Copyright (C) 2006 Jean-Marc Huré.
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License
! as published by the Free Software Foundation; either version 2
! of the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! Web site currently at http://www.gnu.org/licenses/gpl.html
!
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
! Author address:   Jean-Marc Huré,
!                   L3AB/OASU,
!                   2, rue de l'Observatoire,
!                   B.P. 89,
!                   33270 Floirac,
!                   France
!
!             and   Université Bordeaux 1
!                   351, cours de la Libération,
!                   33405 Talence cedex
!                   France
!
!                   email: jean-marc.hure@obs.u-bordeaux1.fr
!

MODULE APPLAWDPACK

  ! global variables

  ! actual precision
  ! set AP=SP for simple precision (1.1920929E-07)
  ! set AP=DP for double precision (2.220446049250313E-016)
  ! set AP=QP for quadruple precision (1.925929944387235853055977942584927E-0034)
  !     (may not work on all computers)
  INTEGER,PARAMETER::SP=KIND(1.00E+00)
  INTEGER,PARAMETER::DP=KIND(1.00D+00)
  ! comment next line if quadruple precision is not supported by your computer,
  ! uncomment otherwise
  !INTEGER,PARAMETER::QP=KIND(1.00Q+00)
  INTEGER,PARAMETER::AP=DP

  ! mathematical constants at the actual precision
  !  PI=3.1415...
  !  TWOPI=2*PI
  !  FOURPI=4*PI
  !  PIOVER2=PI/2
  !  PI2=PI*PI
  !  SQRPI=sqare root of PI
  !  ONETHIRD=1/3
  REAL(KIND=AP)::PI,TWOPI,FOURPI,PIOVER2,PI2,SQRPI,ONETHIRD
  REAL(KIND=AP)::INFINITY,EPSMACH

CONTAINS

  SUBROUTINE WHICHPRECISION()
    !
    ! WHICHPRECISION is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! prints the informations connected with the actual precision
    IMPLICIT NONE
    ! statements
    PRINT '(/,1X,"ACTUAL PRECISION F90, real kind",1X,I2)',AP
    PRINT '(1X,"Maximal precision",1X,SP,1PE21.14,SS)',EPSILON(0._AP)
  END SUBROUTINE WHICHPRECISION

  SUBROUTINE MATHS()
    !
    ! MATHS is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! defines a few mathematical constants at the actual precision,
    IMPLICIT NONE
    ! statements
    PIOVER2=ACOS(0._AP)
    PI=ACOS(-1._AP)
    TWOPI=PI+PI
    FOURPI=2*TWOPI
    PI2=PI*PI
    SQRPI=SQRT(PI)
    ONETHIRD=1._AP/3._AP
    INFINITY=HUGE(0._AP)
    EPSMACH=EPSILON(0._AP)
  END SUBROUTINE MATHS

  FUNCTION SURFDENSITY1(RINTERN,REXTERN,SIGMARIN,A,S)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! surface density as a power law
    IMPLICIT NONE
    REAL(KIND=AP),INTENT(IN)::RINTERN,REXTERN,SIGMARIN,A,S
    REAL(KIND=AP)::SURFDENSITY1
    ! statements
    IF ((A>=RINTERN).AND.(A<=REXTERN)) THEN
       SURFDENSITY1=SIGMARIN*(A/RINTERN)**S
    ELSE
       SURFDENSITY1=0._AP
    ENDIF
  END FUNCTION SURFDENSITY1

  FUNCTION EIK(M)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! complete elliptical integral of the first kind K(M) with modulus M
    ! adapted from the Slatec numerical library
    IMPLICIT NONE
    REAL(KIND=AP),INTENT(IN)::M
    REAL(KIND=AP)::EIK
    ! local variables
    REAL(KIND=AP)::M12,DRFVP
    REAL(KIND=AP)::LOLIM,UPLIM,EPS,ERRTOL
    REAL(KIND=AP)::C1,C2,C3,E2,E3,LAMDA
    REAL(KIND=AP)::MU,S,X,XN,XNDEV
    REAL(KIND=AP)::XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
    ! statements
    IF (M==1._AP) THEN
       PRINT*,'***EIK::ATTEMPT TO COMPUTE K(1)',HUGE(0._AP)
       STOP
    ELSE
       M12=1._AP-M**2
       DRFVP=0._AP
       ERRTOL=(2._AP*EPSILON(0._AP))**(1._AP/6._AP)
       LOLIM=5._AP*TINY(0._AP)
       UPLIM=HUGE(0._AP)/5._AP
       X=0._AP
       Y=M12
       Z=1._AP
       IF (MIN(X,Y,Z)<0._AP) STOP
       IF (MAX(X,Y,Z)>UPLIM) STOP
       IF (MIN(X+Y,X+Z,Y+Z)<LOLIM) STOP
       C1=1._AP/24._AP
       C2=3._AP/44._AP
       C3=1._AP/14._AP
       XN=X
       YN=Y
       ZN=Z
30     MU=(XN+YN+ZN)/3._AP
       XNDEV=2._AP-(MU+XN)/MU
       YNDEV=2._AP-(MU+YN)/MU
       ZNDEV=2._AP-(MU+ZN)/MU
       EPS=MAX(ABS(XNDEV),ABS(YNDEV),ABS(ZNDEV))
       IF (EPS<ERRTOL) GOTO 40
       XNROOT=SQRT(XN)
       YNROOT=SQRT(YN)
       ZNROOT=SQRT(ZN)
       LAMDA=XNROOT*(YNROOT+ZNROOT)+YNROOT*ZNROOT
       XN=(XN+LAMDA)*0.25_AP
       YN=(YN+LAMDA)*0.25_AP
       ZN=(ZN+LAMDA)*0.25_AP
       GOTO 30
40     E2=XNDEV*YNDEV-ZNDEV*ZNDEV
       E3=XNDEV*YNDEV*ZNDEV
       S=1._AP+(C1*E2-0.1_AP-C2*E3)*E2+C3*E3
       DRFVP=S/SQRT(MU)
       EIK=DRFVP
    ENDIF
  END FUNCTION EIK

  FUNCTION EIE(M)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! complete elliptical integral of the second kind E(M) with modulus M
    ! adapted from the Slatec numerical library
    IMPLICIT NONE
    REAL(KIND=AP),INTENT(IN)::M
    REAL(KIND=AP)::EIE
    ! local variables
    REAL(KIND=AP)::M12,DRDVP
    REAL(KIND=AP)::LOLIM,TUPLIM,UPLIM,EPS,ERRTOL
    REAL(KIND=AP)::C1,C2,C3,C4,EA,EB,EC,ED,EF,LAMDA
    REAL(KIND=AP)::MU,POWER4,SIGMA,S1,S2,X,XN,XNDEV
    REAL(KIND=AP)::XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
    ! statements
    IF (M==1._AP) THEN
       EIE=1._AP
    ELSE
       M12=1._AP-M**2
       DRDVP=0._AP
       ERRTOL=(EPSILON(0._AP)/6._AP)**(1._AP/6._AP)
       LOLIM=2._AP/(HUGE(0._AP))**(2._AP/3._AP)
       TUPLIM=TINY(0._AP)**(1._AP/3._AP)
       TUPLIM=(0.1_AP*ERRTOL)**(1._AP/3._AP)/TUPLIM
       UPLIM=TUPLIM*TUPLIM
       C1=3._AP/14._AP
       C2=1._AP/6._AP
       C3=9._AP/22._AP
       C4=3._AP/26._AP
       X=0._AP
       Y=M12
       Z=1._AP
       IF (MIN(X,Y,Z)<0._AP) STOP
       IF (MAX(X,Y,Z)>UPLIM) STOP
       IF (MIN(X+Y,X+Z,Y+Z)<LOLIM) STOP
       XN=X
       YN=Y
       ZN=Z
       SIGMA=0._AP
       POWER4=1._AP
30     MU=(XN+YN+3._AP*ZN)*0.2_AP
       XNDEV=(MU-XN)/MU
       YNDEV=(MU-YN)/MU
       ZNDEV=(MU-ZN)/MU
       EPS=MAX(ABS(XNDEV),ABS(YNDEV),ABS(ZNDEV))
       IF (EPS<ERRTOL) GOTO 40
       XNROOT=SQRT(XN)
       YNROOT=SQRT(YN)
       ZNROOT=SQRT(ZN)
       LAMDA=XNROOT*(YNROOT+ZNROOT)+YNROOT*ZNROOT
       SIGMA=SIGMA+POWER4/(ZNROOT*(ZN+LAMDA))
       POWER4=POWER4*0.25_AP
       XN=(XN+LAMDA)*0.25_AP
       YN=(YN+LAMDA)*0.25_AP
       ZN=(ZN+LAMDA)*0.25_AP
       GOTO 30
40     EA=XNDEV*YNDEV
       EB=ZNDEV*ZNDEV
       EC=EA-EB
       ED=EA-6._AP*EB
       EF=ED+EC+EC
       S1=ED*(-C1+0.25_AP*C3*ED-1.5_AP*C4*ZNDEV*EF)
       S2= ZNDEV*(C2*EF+ZNDEV*(-C3*EC+ZNDEV*C4*EA))
       DRDVP=3._AP*SIGMA+POWER4*(1._AP+S1+S2)/(MU*SQRT(MU))
       EIE=EIK(M)-M**2*DRDVP/3._AP
    ENDIF
  END FUNCTION EIE

  FUNCTION PSIHOMO(R,RINTERN,REXTERN,SIGMA0)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! potential in the plane of a flat axially symmetric homogeneous disk
    IMPLICIT NONE
    REAL(KIND=AP),INTENT(IN)::R,SIGMA0,RINTERN,REXTERN
    REAL(KIND=AP)::PSIHOMO
    ! local variables
    REAL(KIND=AP)::ROVERROUT,RINOVERR,ROUTOVERR,ROVERRIN
    ! statements
    IF (SIGMA0==0._AP) THEN
       PSIHOMO=0._AP
    ELSE
       IF (R==0._AP) THEN
          PSIHOMO=-2._AP*PI*(REXTERN-RINTERN)*SIGMA0
       ELSE
          IF (R==RINTERN) THEN
             ROVERROUT=R/REXTERN
             PSIHOMO=-4._AP*R*SIGMA0*(EIE(ROVERROUT)/ROVERROUT-1._AP)
          ELSE
             IF (R==REXTERN) THEN
                RINOVERR=RINTERN/R
                PSIHOMO=-4._AP*R*SIGMA0*(1._AP-EIE(RINOVERR) &
                     +(1._AP-RINOVERR**2)*EIK(RINOVERR))
             ELSE
                IF (R<RINTERN) THEN
                   ROVERROUT=R/REXTERN
                   ROVERRIN=R/RINTERN
                   PSIHOMO=-4._AP*R*SIGMA0 &
                        *(EIE(ROVERROUT)/ROVERROUT-EIE(ROVERRIN)/ROVERRIN)
                ELSE
                   IF (R>REXTERN) THEN
                      ROUTOVERR=REXTERN/R
                      RINOVERR=RINTERN/R
                      PSIHOMO=-4._AP*R*SIGMA0*( &
                           EIE(ROUTOVERR)-(1._AP-ROUTOVERR**2)*EIK(ROUTOVERR) &
                           -(EIE(RINOVERR)-(1._AP-RINOVERR**2)*EIK(RINOVERR)))
                   ELSE
                      ROVERROUT=R/REXTERN
                      RINOVERR=RINTERN/R
                      PSIHOMO=4._AP*R*SIGMA0*(EIE(RINOVERR) &
                           -(1._AP-RINOVERR**2)*EIK(RINOVERR) &
                           -EIE(ROVERROUT)/ROVERROUT)
                   ENDIF
                ENDIF
             ENDIF
          ENDIF
       ENDIF
    ENDIF
  END FUNCTION PSIHOMO

  SUBROUTINE APPLAWD(RINTERN,REXTERN,SIGMARIN,S,R,ORDER,PSI)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    IMPLICIT NONE
    REAL(KIND=AP),INTENT(IN)::R,RINTERN,REXTERN,SIGMARIN,S
    INTEGER(KIND=4),INTENT(IN)::ORDER
    REAL(KIND=AP),INTENT(OUT)::PSI
    ! local variables
    REAL(KIND=AP)::ANPSIVALUE,SIGMA0
    ! statement
    IF (R==0._AP) THEN
       IF (S==-1._AP) THEN
          ANPSIVALUE=-TWOPI*SURFDENSITY1(RINTERN,REXTERN,SIGMARIN,REXTERN,S)*REXTERN*LOG(REXTERN/RINTERN)
       ELSE
          ANPSIVALUE=-TWOPI/(1._AP+S)*(SURFDENSITY1(RINTERN,REXTERN,SIGMARIN,REXTERN,S)*REXTERN-SURFDENSITY1(RINTERN,REXTERN,SIGMARIN,RINTERN,S)*RINTERN)
       ENDIF
    ELSE
       ! local surface density
       IF (R>REXTERN) THEN
          SIGMA0=SURFDENSITY1(RINTERN,REXTERN,SIGMARIN,REXTERN,S)
       ELSE
          IF (R>=RINTERN) THEN
             ! local value
             SIGMA0=SURFDENSITY1(RINTERN,REXTERN,SIGMARIN,R,S)
          ELSE
             SIGMA0=SURFDENSITY1(RINTERN,REXTERN,SIGMARIN,RINTERN,S)
          ENDIF
       ENDIF
       ! homogeneous component
       ANPSIVALUE=PSIHOMO(R,RINTERN,REXTERN,SIGMA0)
       ! residual component
       ANPSIVALUE=ANPSIVALUE+DELTAPSI(R,RINTERN,REXTERN,SIGMA0,S,ORDER)
    ENDIF
    PSI=ANPSIVALUE
  END SUBROUTINE APPLAWD

  FUNCTION DELTAPSI(R,RINTERN,REXTERN,SIGMA0,S,ORDER)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! residual potential
    IMPLICIT NONE
    INTEGER(KIND=4),INTENT(IN)::ORDER
    REAL(KIND=AP),INTENT(IN)::R,RINTERN,REXTERN,SIGMA0,S
    ! local variables
    REAL(KIND=AP)::VALUE,DELTAPSI,GAMMAI
    INTEGER(4)::I
    ! statement
    VALUE=0._AP
    DO I=0,ORDER
       IF (I==0) THEN
          GAMMAI=1._AP
       ELSE
          GAMMAI=GAMMAI*(REAL(2*I-1,AP)/REAL(2*I,AP))**2
       ENDIF
       IF (R>REXTERN) THEN
          VALUE=VALUE+GAMMAI*EDELTAPSIN(R,RINTERN,REXTERN,S,I)
       ELSE
          IF (R>=RINTERN) THEN
             VALUE=VALUE+GAMMAI*DELTAPSIN(R,RINTERN,REXTERN,S,I)
          ELSE
             VALUE=VALUE+GAMMAI*IDELTAPSIN(R,RINTERN,REXTERN,S,I)
          ENDIF
       ENDIF
    END DO
    DELTAPSI=TWOPI*SIGMA0*R*VALUE
  END FUNCTION DELTAPSI

  FUNCTION DELTAPSIN(R,RINTERN,REXTERN,S,ORDER)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! terms of the series; inside the disk
    IMPLICIT NONE
    INTEGER(KIND=4),INTENT(IN)::ORDER
    REAL(KIND=AP),INTENT(IN)::R,RINTERN,REXTERN,S
    ! local variables
    REAL(KIND=AP)::INDEX1,INDEX2,INDEX3,INDEX4,UIN,VOUT,DELTAPSIN
    ! statement
    ! normalized edges
    UIN=RINTERN/R
    VOUT=R/REXTERN
    ! indexes leading to special integrals
    INDEX1=REAL(2*(ORDER+1),AP)
    INDEX2=REAL(2*ORDER-1,AP)-S
    INDEX3=REAL(2*ORDER-1,AP)
    INDEX4=S+REAL(2*(ORDER+1),AP)
    IF (INDEX4==0) THEN
       DELTAPSIN=(1._AP-UIN**INDEX1)/INDEX1+LOG(UIN)&
            +(VOUT**INDEX2-1._AP)/INDEX2-(VOUT**INDEX3-1._AP)/INDEX3
    ELSE
       IF (INDEX2==0) THEN
          DELTAPSIN=(1._AP-UIN**INDEX1)/INDEX1-(1._AP-UIN**INDEX4)/INDEX4&
               +LOG(VOUT)-(VOUT**INDEX3-1._AP)/INDEX3
       ELSE
          DELTAPSIN=(1._AP-UIN**INDEX1)/INDEX1-(1._AP-UIN**INDEX4)/INDEX4&
               +(VOUT**INDEX2-1._AP)/INDEX2-(VOUT**INDEX3-1._AP)/INDEX3
       ENDIF
    ENDIF
  END FUNCTION DELTAPSIN

  FUNCTION EDELTAPSIN(R,RINTERN,REXTERN,S,ORDER)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! terms of the series; rightwards to the outer edge
    IMPLICIT NONE
    INTEGER(KIND=4),INTENT(IN)::ORDER
    REAL(KIND=AP),INTENT(IN)::R,RINTERN,REXTERN,S
    ! local variables
    REAL(KIND=AP)::INDEX1,INDEX4,UIN,UOUT,VOUT,EDELTAPSIN
    ! statement
    ! normalized edges
    UIN=RINTERN/R
    UOUT=REXTERN/R
    VOUT=R/REXTERN
    ! indexes leading to special integrals
    INDEX1=REAL(2*(ORDER+1),AP)
    INDEX4=S+REAL(2*(ORDER+1),AP)
    IF (INDEX4==0) THEN
       EDELTAPSIN=(UOUT**INDEX1-UIN**INDEX1)/INDEX1+LOG(UIN/UOUT)*VOUT**S
    ELSE
       EDELTAPSIN=(UOUT**INDEX1-UIN**INDEX1)/INDEX1-(UOUT**INDEX4-UIN**INDEX4)/INDEX4*VOUT**S
    ENDIF
  END FUNCTION EDELTAPSIN

  FUNCTION IDELTAPSIN(R,RINTERN,REXTERN,S,ORDER)
    !
    ! APPLAWD is part of APPLAWDPACK
    !
    ! APPLAWDPACK is a Fortran-90 module for APPLAWD.
    !
    ! APPLAWD determines the gravitational potential inside and outside
    ! a flat axially symmetric disk with finite size and power law surface
    ! density profile from the splitting method. Potential values are computed
    ! from a truncated series so that accuracy is easily controlled through the
    ! expansion order.
    ! 
    ! see Huré, Pelat & Pierens (2006), A&A, submitted
    ! 
    ! Copyright (C) 2006 Jean-Marc Huré.
    !
    ! This program is free software; you can redistribute it and/or
    ! modify it under the terms of the GNU General Public License
    ! as published by the Free Software Foundation; either version 2
    ! of the License, or (at your option) any later version.
    !
    ! This program is distributed in the hope that it will be useful,
    ! but WITHOUT ANY WARRANTY; without even the implied warranty of
    ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    ! GNU General Public License for more details.
    !
    ! Web site currently at http://www.gnu.org/licenses/gpl.html
    !
    ! You should have received a copy of the GNU General Public License
    ! along with this program; if not, write to the Free Software
    ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
    !
    ! Author address:   Jean-Marc Huré,
    !                   L3AB/OASU,
    !                   2, rue de l'Observatoire,
    !                   B.P. 89,
    !                   33270 Floirac,
    !                   France
    !
    !             and   Université Bordeaux 1
    !                   351, cours de la Libération,
    !                   33405 Talence cedex
    !                   France
    !
    !                   email: jean-marc.hure@obs.u-bordeaux1.fr
    !
    ! terms of the series; leftwards to the inner edge
    IMPLICIT NONE
    INTEGER(KIND=4),INTENT(IN)::ORDER
    REAL(KIND=AP),INTENT(IN)::R,RINTERN,REXTERN,S
    ! local variables
    REAL(KIND=AP)::INDEX2,INDEX3,VIN,VOUT,IDELTAPSIN
    ! statement
    ! normalized edges
    VIN=R/RINTERN
    VOUT=R/REXTERN
    ! indexes leading to special integrals
    INDEX2=REAL(2*ORDER-1,AP)-S
    INDEX3=REAL(2*ORDER-1,AP)
    IF (INDEX2==0) THEN
       IDELTAPSIN=LOG(VOUT/VIN)*VIN**S-(VOUT**INDEX3-VIN**INDEX3)/INDEX3
    ELSE
       IDELTAPSIN=(VOUT**INDEX2-VIN**INDEX2)/INDEX2*VIN**S-(VOUT**INDEX3-VIN**INDEX3)/INDEX3
    ENDIF
  END FUNCTION IDELTAPSIN

END MODULE APPLAWDPACK