PROGRAM PLSAPT
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
C
C* Begin Prologue PLSAPT
C
C ---------------------------------------------------------------------
C
C* Purpose           Post-processor using the global solution generated
C                    by a previous run of MEXX with
C                    example pendulum (see driver PLSADR)
C                     - dense output mode will be read from stdin
C
C* File              plsapt.f
C* Version           1.1
C* Latest Change     96/07/19 (1.4)
C* Modification history
C    1.1             New free format starts text output with ' '.
C* Library           CodeLib
C* Code              Fortran 77
C                    Double Precision
C* Environment       Standard version for FORTRAN77 environments on
C                    PCs, workstations, and hosts
C
C* Subroutines supplied
C
C     PLSIN          Provide initial values
C     PLDOUT         Output interpolated solution values
C      
C  Problem dimensions for a maximum of 30 bodies.
C  The actual number of bodies will be requested from standard input.
C
      PARAMETER
     $     ( NBODY  = 30
     $     , NPDIM  = 3*NBODY
     $     , NVDIM  = NPDIM
     $     , NUDIM  = 1
     $     , NGDIM  = NVDIM
     $     , NLDIM  = NGDIM
     $     , NZGDIM = 8*NBODY - 4
     $     , NPTSMX = 2000
     $     , NDMAX  = 1000
     $     , KDMAX  = 15
     $     , LDENS  = NDMAX*KDMAX
     $     )
C
C  Workspace dimension
C
      EXTERNAL PLSIN, PLDOUT
C
      LOGICAL QPRTIM
      DIMENSION YIP(NDMAX), DOUTS(NPTSMX), DENS(LDENS)
      DIMENSION INFO(50), INFOH(50)
      DIMENSION INDP(NPDIM), INDV(NVDIM), INDA(NVDIM), INDU(NUDIM),
     $     INDL(NLDIM)
      DIMENSION P(NPDIM), V(NVDIM), U(NUDIM), A(NVDIM), RLAM(NLDIM)
C
      COMMON / PEND / PMASS(NBODY), PLEN(NBODY), NPEND
C
      CHARACTER PROTXT*30, FILNAM*10
C
      SAVE
C
      CHARACTER*67 CMODUL
      DATA CMODUL /'@(#) plsapt.f 1.4 from 96/07/19 (mexx 1.2)
     $     '/
C
      DATA FILNAM /'PLSADR.GLO'/
C
C  Machine constants (logical units stdin and stdout)
      LUSIN = 5
      LUSOUT = 6
      LUSTAT = 0
C
C -- Start message
      WRITE (LUSOUT,*) ' '
      WRITE (LUSOUT,*)
     $     '===================================================='
      WRITE (LUSOUT,*) CMODUL
      WRITE (LUSOUT,*)
     $     '===================================================='
      WRITE (LUSOUT,*) ' '
C
C -- Initialize problem
C
C
C
      NP = NPDIM
      NV = NVDIM
      NL = NLDIM
      NG = NGDIM
      NU = NUDIM
      NBLK = 0
      NMRC = 0
      MMODE = 0
      NZGMAX = NZGDIM
      NZFMAX = 0
C
      ICALL=0
      DO 1000 I=1,50
 1000 INFO(I)=0
      QPRTIM=.FALSE.
C
      CALL PLSIN (NP, NV, NL, NG, NU, NBLK, NMRC, NZGMAX,
     $     NZFMAX, NSWMAX, T0, P, V, U, A, RLAM, PROTXT)
C
C  Loop over integration points
 1010 CONTINUE
C
      CALL MXDRE1(FILNAM, INFOH, ND, NDP, NDV, NDU, NDA, NDR,
     $            KDENS, ILAM, IRHO, T, TOLD, H, DENS,
     $            YIP, T0, TFIN, DOUTS, ICALL,
     $            INDP, INDV, INDU, INDA, INDL, IRTRN)
      ICALL=ICALL+1
      INFO(41)=INFOH(41)
      INFO(42)=INFOH(42)
      IF(ICALL.EQ.1) THEN
         INFO(31)=INFOH(31)
         INFO(32)=INFOH(32)
         INFO(33)=INFOH(33)
         INFO(3)=INFOH(3)
         MDOUT = INFOH(31)
         NDOUT = INFOH(32)
      ENDIF
      IF(IRTRN .NE. 0) THEN
         WRITE (LUSOUT,*)
     $        '- PLSAPT - eof reached or read error in ', FILNAM
         GOTO 1120
      ENDIF
C
      IF(ICALL.EQ.1) THEN
         WRITE (LUSOUT,*) 'The dense output mode for MEXX was:',MDOUT
         WRITE (LUSOUT,*) 'Enter zero to use this mode'
C
         WRITE (LUSOUT,*) 'Enter MDOUT (0):'
C         READ (LUSIN,*, ERR=1020, END=1020) MDOUTH
C         GOTO 1030
 1020    MDOUTH = 2
 1030    WRITE (LUSOUT,*) 'Given:                   ',
     $        MDOUTH
         WRITE (LUSOUT,*) ' '
         IF(MDOUTH.EQ.0) GOTO 1110
         MDOUT=MDOUTH
         IF(MDOUT.EQ.3) THEN
            DOUTS(1) = 10.D0
            WRITE (LUSOUT,*) 'Enter delta-t-max (', DOUTS(1), '):'
            READ (LUSIN,*, ERR=1040, END=1040) DOUTS(1)
            GOTO 1050
 1040       DOUTS(1) = 10.D0
 1050       WRITE (LUSOUT,*) 'Given:                   ',
     $           DOUTS(1)
            WRITE (LUSOUT,*) ' '
         ENDIF
         IF(MDOUT.NE.3) THEN
            WRITE (LUSOUT,*) 'Enter NDOUT(9):'
C            READ (LUSIN,*, ERR=1060, END=1060) NDOUT
C            GOTO 1070
 1060       NDOUT = 2
 1070       WRITE (LUSOUT,*) 'Given:                   ',
     $           NDOUT
            WRITE (LUSOUT,*) ' '
         ENDIF
         IF(MDOUT.EQ.4) THEN
            WRITE (LUSOUT,*) 'Enter halting points:'
            DO 1100 I=1,NDOUT
 1080          WRITE (LUSOUT,*) 'No.:',I
               READ (LUSIN,*, ERR=1090, END=1090) DOUTS(I)
               WRITE (LUSOUT,*) 'Given:                   ',
     $              DOUTS(I)
               WRITE (LUSOUT,*) ' '
               GOTO 1100
 1090          CONTINUE
               WRITE (LUSOUT,*)
     $              '+++ PLSAPT: Illegal input for halting points.'
               INP = INP + 1
               IF (INP .GE. MAXINP) THEN
                  WRITE (LUSOUT,*)
     $                 '            You have had', INP,
     $                 'tries. Program stopped.'
                  STOP
               ENDIF
               GOTO 1080
 1100       CONTINUE
         ENDIF
         INFO(31)=MDOUT
         INFO(32)=NDOUT
      ENDIF
C
 1110 CONTINUE
C
      CALL MXDOUT (INFO, ND, NDP, NDV, NDU, NDA, NDR,
     $            KDENS, ILAM, IRHO, T, TOLD, H, DENS,
     $            YIP, T0, TFIN, DOUTS, PLDOUT,
     $            ND, NDP, NDV, NDU, NDA, NDR,
     $            INDP, INDV, INDU, INDA, INDL, QPRTIM, IRTRN)
C
      IF (INFOH(41) .NE. 2 .AND. INFOH(41) .NE. 4) GOTO 1010
C
 1120 CONTINUE
      WRITE (LUSOUT,*) ' '
      WRITE (LUSOUT,*) 'End of post processing for example ', PROTXT
      WRITE (LUSOUT,*) 'Bye'
C
      STOP
C
 9000 FORMAT (2A)
 9100 FORMAT (A, 3I8)
 9200 FORMAT (A, 1PD8.1)
 9210 FORMAT ((4(1PD19.10)))   
C
C  end of program PLSAPT
C
      END
C
      SUBROUTINE PLSIN (NP, NV, NL, NG, NU, NBLK, NMRC, NZGMAX,
     $     NZFMAX, NSWMAX, T0, P, V, U, A, RLAM, PROTXT)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      INTEGER NP, NV, NL, NG, NU, NBLK, NMRC, NZGMAX, NZFMAX
      DIMENSION P(*), V(*), U(*), A(*), RLAM(*)
      CHARACTER*(*) PROTXT
C
C **********************************************************************
C
C   Pendulum (variable number of bodies < 31)
C
C
C   Parameter:  masses                  pmass(i) 
C               lengths                 plen(i)  
C
C
C **********************************************************************
C
C* File              plsin.f
C* Version           1.1
C* Latest Change     96/07/19 (1.4)
C* Modification history
C    1.1             Correct initial values of RLAM.
C                    New free format starts text output with ' '.
C
C -------------------------------------------------------
C
      PARAMETER (NBMAX = 30)
      COMMON / PEND / PMASS(NBMAX), PLEN(NBMAX), NPEND
C
      SAVE
C
      CHARACTER*67 CMODUL
      DATA CMODUL /'@(#) plsin.f 1.4 from 96/07/19 (mexx 1.2)
     $     '/
C
      PROTXT = 'Pendulum A (sparse)'
C
C  Machine constants (logical units stdin and stdout)
      LUSIN = 5
      LUSOUT = 6
C
      WRITE (LUSOUT,*) 'Start example ', PROTXT
      WRITE (LUSOUT,*) ' '
C
C  Dimensions
C  Get problem parameters
C
      WRITE (LUSOUT,*) 'Enter number of bodies (default 8):'
C      READ (LUSIN,*, ERR=1000, END=1000) NPEND
C      GOTO 1010
 1000 NPEND = 8
 1010 WRITE (LUSOUT,*) 'Given:', NPEND
      IF (NPEND .LT. 1 .OR. NPEND .GT. NBMAX) THEN
         WRITE (LUSOUT,*)
     $        'Number of bodies is out of range (default value chosen)'
         NPEND = 8
      ENDIF
      WRITE (LUSOUT,*) ' '
C
      NP = 3*NPEND
      NV = NP
      NG = 2*NPEND
      NL = NG
      NU = 0
      NMRC = 1
      NBLK = NP
      NZGMAX = 8*NPEND - 4
      NZFMAX = 0
      NSWMAX = NPEND - 1
C
C  Model parameters
      T0 = 0.D0
      PMASSH = 1.0D0/DBLE(NPEND)
      PLENH = 0.5D0/DBLE(NPEND)
      DO 1020 I=1,NPEND
         PMASS(I) = PMASSH
         PLEN(I) = PLENH  
 1020 CONTINUE
C     
C  Initial values
      PIHALF = ASIN(1.D0)
      P(1) = PLEN(1)
      P(2) = 0.D0 
      P(3) = PIHALF  
C     
      V(1) = 0.D0
      V(2) = 0.D0
      V(3) = 0.D0
C     
      RLAM(1) = 0.D0
      RLAM(2) = 0.D0
C     
      DO 1030 I=2,NPEND
         L = (I-1)*3 + 1
         P(L) = P(L-3) + PLEN(I-1) + PLEN(I)
         P(L+1) = P(2) 
         P(L+2) = P(3)
         V(L) = 0.D0
         V(L+1) = 0.D0
         V(L+2) = 0.D0
         L = (I-1)*2 + 1
         RLAM(L) = 0.D0
         RLAM(L+1) = 0.D0
 1030 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE PLDOUT (NDP, NDV, NDU, NDA, NDL, T,
     $                   DP, DV, DU, DA, DL, INDP,
     $                   INDV, INDU, INDA, INDL, INFO, IFAIL)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      INTEGER NDP, NDV, NDU, NDA, NDL
      DOUBLE PRECISION T, DP(*), DV(*), DU(*), DA(*), DL(*)
      INTEGER INDP(*), INDV(*), INDU(*), INDA(*), INDL(*), INFO(50),
     $     IFAIL 
C
C
C* File              pldout.f
C* Version           1.1
C* Latest Change     96/04/30 (1.5)
C* Modification history
C    1.1             Avoid free format because of compiler dependency.
C      
C* Input arguments:
C
C   NDP       Int   Number of dense output components of position 
C                   vector P. Dimension of arrays DP and INDP. The
C                   current value is either the user prescibed value
C                   IWK(51) (if MDIND=INFO(33)=1 was set) or NP, the
C                   dimension of the position vector P (set by MEXX if
C                   INFO(33)=0 was set).
C
C   NDV       Int   as NDP but now for the velocity vector V.
C
C   NDU       Int   as NDP but now for the vector of dynamic 
C                   variables U.
C
C   NDA       Int   as NDP but now for the acceleration vector A.
C
C   NDL       Int   as NDP but now for the vector of Lagrange 
C                   multipliers.
C
C   T         Dble  Current time
C
C   DP        Dble  the dense output position vector. 
C
C   DV        Dble  the dense output velocity vector. 
C
C   DU        Dble  the dense output dynamic variable vector. 
C
C   DA        Dble  the dense output acceleration vector. 
C
C   DL        Dble  the dense output Lagrange multiplier vector. 
C
C   INDP(NDP) Int   indicator array of length NP. List of components of
C                   the position vecor p which was selected for dense
C                   output. INDP(I)=J means that component J of the
C                   vector P is available as component I of vector DP.
C                   (I.e. DP(I)=P(INDP(I))
C                   Note, that MEXX automatically provides INDP(I)=I
C                   if INFO(33)=0 was set.
C
C   INDV(NDV) Int   as INDP but now for the velocity vector V.
C
C   INDU(NDU) Int   as INDP but now the vector of dynamic variables U.
C
C   INDU(NDA) Int   as INDP but now for the acceleration vector A.
C
C   INDL(NDL) Int   as INDP but now the vector of Lagrange multipliers.
C  
C   INFO(50)  Int   Info array
C                   INFO(1) .. INFO(40) see MXJOB(1) .. MXJOB(40) 
C                                           in MEXX
C                   INFO(41): Type of call: 
C                             = 0: the first call
C                             = 1: an intermediate call
C                             = 2: the last call at TFIN
C                             = 3: at a root of switching function
C                                  and integration continues
C                             = 4: at a root of switching function
C                                  and integration will terminate
C                                  (this is the last call)
C
C                   INFO(42): number of current integration step
C
C                   INFO(43): interpolation indicator
C                             = 0: current t is an integration point
C                             = 1: the solution values are interpolated
C                                  values
C
C                   INFO(44): call counter
C
C   IFAIL     Int   Error indicator. IFAIL=0 on input
C   
C   Do NOT alter any of the input arguments.
C
C
C* Output arguments:
C
C   IFAIL     Int  Error indicator
C                  If IFAIL.LT.0  on return, MEXX will terminate
C 
C----------------------------------------------------------------------
C
C
      INTEGER LUNIT
C
      CHARACTER*11 FNAME
C
      PARAMETER ( LUNIT = 11 , FNAME = 'PLSADR.DNO' )
C
      SAVE
C
      IFAIL = 0
C
      IF(INFO(41).EQ.0) THEN
         OPEN (UNIT=LUNIT,FILE=FNAME)
         WRITE(LUNIT,'(A)') 
     $    'C Dense Output Data from MEXX by routine pldout'
         WRITE(LUNIT,'(A)') 
     $    'C  the user selected dense output options are:'
         WRITE(LUNIT,'(A)') 'C    MDOUT     NDOUT     MDIND     MMODE'
         WRITE(LUNIT,'(3I10)') INFO(31),INFO(32),INFO(33),INFO(3)
         WRITE(LUNIT,'(A)')  'C  the dense output dimensions are:'
         WRITE(LUNIT,'(A)') 'C NDP  NDV  NDU  NDA  NDL'
         WRITE(LUNIT,'(5I5)') NDP,NDV,NDU,NDA,NDL 
         WRITE(LUNIT,'(A)') 'C  the selected components of P are:'
         WRITE(LUNIT,'(10I7)') (INDP(I),I=1,NDP)   
         WRITE(LUNIT,'(A)') 'C  the selected components of V are:'
         WRITE(LUNIT,'(10I7)') (INDV(I),I=1,NDV)   
         WRITE(LUNIT,'(A)') 'C  the selected components of U are:'
         WRITE(LUNIT,'(10I7)') (INDU(I),I=1,NDU)   
         WRITE(LUNIT,'(A)') 'C  the selected components of A are:'
         WRITE(LUNIT,'(10I7)') (INDA(I),I=1,NDV)   
         WRITE(LUNIT,'(A)') 'C  the selected components of LAMBDA are:'
         WRITE(LUNIT,'(10I7)') (INDL(I),I=1,NDL)   
      ENDIF
C
      WRITE(LUNIT,'(A)') 'C  No                   T Type Ipol'
      WRITE(LUNIT,'(I5,1PD20.10,I5,I5)') INFO(44),T,INFO(41),INFO(43)
      WRITE(LUNIT,'(A)') 'C  DP:'
      WRITE(LUNIT,'(4(1PD20.10))') (DP(I),I=1,NDP)   
      WRITE(LUNIT,'(A)') 'C  DV:'
      WRITE(LUNIT,'(4(1PD20.10))') (DV(I),I=1,NDV)   
      WRITE(LUNIT,'(A)') 'C  DU:'
      WRITE(LUNIT,'(4(1PD20.10))') (DU(I),I=1,NDU)   
      WRITE(LUNIT,'(A)') 'C  DA:'
      WRITE(LUNIT,'(4(1PD20.10))') (DA(I),I=1,NDA)   
      WRITE(LUNIT,'(A)') 'C  DL:'
      WRITE(LUNIT,'(4(1PD20.10))') (DL(I),I=1,NDL)   
C
      IF(INFO(41).EQ.2) THEN
         WRITE(LUNIT,'(A)') 
     $    'C MEXX terminates at final point of integration'
      ENDIF
      IF(INFO(41).EQ.4) THEN
         WRITE(LUNIT,'(A)') 
     $    'C MEXX terminates at a root of the switching functions'
      ENDIF
      IF(INFO(41).EQ.2 .OR. INFO(41).EQ.4) THEN
         CLOSE (UNIT=LUNIT)
      ENDIF
C
C
C
C
      RETURN 
C
C end of subroutine PLDOUT
C
      END