!
! DRIVERAPPLAWD is a Fortran-90 driver program for APPLAWD.
! Dependancy : APPLAWDPACK
! 
! 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
!

PROGRAM DRIVERAPPLAWD

  USE APPLAWDPACK

  IMPLICIT NONE

  ! computational boxes
  INTEGER(KIND=4),PARAMETER::NBOX=16
  REAL(KIND=AP),DIMENSION(NBOX)::RBOX
  REAL(KIND=AP)::RBOXIN,RBOXOUT

  ! field point
  REAL(KIND=AP)::R,DELTAR

  ! disk potential
  INTEGER(KIND=4)::ORDER,NEWORDER
  REAL(KIND=AP)::S,PSI,BETTERPSI

  ! disk inner and outer edge, and inner surface density
  REAL(KIND=AP)::SIGMARIN,RIN,ROUT

  ! misc.
  INTEGER(KIND=4)::I

  ! statements 
  PRINT '(/1X,"ADPPLawS Fortran 90 driver program",&
       &/,1X,"GNU GPL protected, V1.0, 2006, Nov.",&
       &/,1X,"Copyright (C) 2006 Jean-Marc Huré")'

  ! prints the actual precision
  CALL WHICHPRECISION()

  ! defines mathematical constants
  CALL MATHS()

  ! opens output file
  OPEN(UNIT=1,FILE="driverapplawd.dat")
  PRINT '(/,1X,"OUTPUT FILE : driverapplawd.dat")'

  ! inner edge
  RIN=0.1_AP

  ! outer edge
  ROUT=1._AP

  ! power law index
  S=-1._AP

  ! surface density at the outer edge
  SIGMARIN=1._AP

  PRINT '(/,1X,"DISK EDGES", &
       & /,1X,"Inner edge",1X,SP,1PE21.14,SS,&
       &/,1X,"Outer edge",1X,SP,1PE21.14,SS)',RIN,ROUT
  PRINT '(/,1X,"SURFACE DENSITY AT THE INNER EDGE:",1X,SP,1PE21.14,SS)',SIGMARIN
  PRINT '(1X,"POWER LAW INDEX S :",1X,SP,1PE21.14,SS)',S

  ORDER=0
  PRINT '(/,1X,"EXPANSION ORDER N :",1X,I)',ORDER

  ! inner grid point
  RBOXIN=0._AP

  ! outer grid point
  RBOXOUT=1.5_AP

  PRINT '(/,1X,"COMPUTATIONAL BOX : ",1X,I6.6,1X,"points")',NBOX
  PRINT '(1X,"From",1X,SP,1PE21.14," to ",1PE21.14,SS)',RBOXIN,RBOXOUT
  IF (NBOX==1) THEN
     RBOX(1)=RBOXIN
  ELSE
     DELTAR=(RBOXOUT-RBOXIN)/REAL(NBOX-1,AP)
     RBOX(1)=RBOXIN
     RBOX(2:NBOX-1)=(/(RBOXIN+REAL(I-1,AP)*DELTAR,I=2,NBOX-1)/)
     RBOX(NBOX)=RBOXOUT
  ENDIF

  PRINT '(/,2X,"R, PSI(R), relative error DPSI/PSI",/,1X,"---")'

  DO I=1,NBOX

     R=RBOX(I)
     ! computes the potential
     CALL APPLAWD(RIN,ROUT,SIGMARIN,S,R,ORDER,PSI)

     ! computes a better approximation (about 10%) for the potential
     NEWORDER=MAX(ORDER+1,INT(REAL(ORDER/0.75_AP,AP)))
     CALL APPLAWD(RIN,ROUT,SIGMARIN,S,R,NEWORDER,BETTERPSI)

     ! prints and save the results
     SELECT CASE (AP)
     CASE (4)
        PRINT 100,R,PSI,1._AP-PSI/BETTERPSI
        WRITE(1,100) R,PSI,1._AP-PSI/BETTERPSI
     CASE(8)
        PRINT 110,R,PSI,1._AP-PSI/BETTERPSI
        WRITE(1,110) R,PSI,1._AP-PSI/BETTERPSI
     CASE(16)
        PRINT 120,R,PSI,1._AP-PSI/BETTERPSI
        WRITE(1,120) R,PSI,1._AP-PSI/BETTERPSI
     END SELECT

  ENDDO

  CLOSE(1)

  PRINT '(1X,"---",//,1X,"Done.",/)'

100 FORMAT(1X,SP,2(1PE14.7,1X),1PE9.2,SS)
110 FORMAT(1X,SP,2(1PE21.14,1X),1PE9.2,SS)
120 FORMAT(1X,SP,2(1PE28.21,1X),1PE9.2,SS)

END PROGRAM DRIVERAPPLAWD