PROGRAM GRAPH
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (LD=50000,LDPAR=300)
      DIMENSION G(LD),W(LD),WT(LD),NV1(LD)
      DIMENSION PAR(LDPAR),WPP(5),A(2)
      CHARACTER*10 HKERN(3),HRANGE(2)
      LOGICAL LORDER
      DATA HKERN/'   Uniform','Triangular','  Biweight'/
      DATA HRANGE/'          ','  Quantile'/
      DATA CON3/.9375D0/
C  Read in and check user information
      WRITE(6,10)
   10 FORMAT(' Preparation of posterior density for plotting')
      READ(5,*)KERN,HFRAC,KRANGE,A(1),A(2),NPTS,NBURN,IFUNC
      CALL PINTI('Kernel type',KERN,1,3)
      CALL PINTR('Smoothing faction parameter',HFRAC,1.0D-12,1.0D0)
      CALL PINTI('Range type',KRANGE,1,2)
      IF(.NOT.LORDER(2,A))CALL PEND('Inverted range limits')
      IF(KRANGE.NE.2)GO TO 20
      CALL PINTR('Range limit',A(1),0.0D0,A(2))
      CALL PINTR('Range limit',A(2),A(1),1.0D0)
   20 CALL PINTI('Points produced',NPTS,1,10000000)
      WRITE(6,30),HKERN(KERN),HFRAC,HRANGE(KRANGE),A(1),A(2),NPTS,
     1  NBURN,IFUNC
   30 FORMAT('                                  Kernel type:',A10,/,
     1       ' Window width fraction of interquartile range:',F10.6,/,       
     2                         20X,A10,' ordinate range:',2(1PE15.5),/,
     3       '  Number of density ordinates written to file:',I10,/,
     4       '                 Number of burn-in iterations:',I10,/,
     5       '                         Function of interest:',I10)
      CALL RDSIM0(1,LDPAR,NITER,NPARS)
      CALL PINTI('Burn-in iterations',NBURN,0,NITER)
      CALL PINTI('Function of interest',IFUNC,1,NPARS)
      CALL FOPEN('Smoothed ordinates',2,.TRUE.)
      NUSE=NITER-NBURN
      NSKIP=(NUSE-1)/LD+1
      IF(NBURN.EQ.0)GO TO 45
      DO 40 IBURN=1,NBURN
   40 CALL RDSIM(1,JTER,WPP,PAR)
C  Pass through file to get weights and values
   45 M=0
      DO 50 ITER=1,NUSE,NSKIP
      M=M+1
      CALL RDSIM(1,JTER,WPP,PAR)
      WTLOG=WPP(1)
      IF(ITER.EQ.1)WTMAX=WTLOG
      IF(WTLOG.GT.WTMAX)WTMAX=WTLOG
      WT(M)=WTLOG
   50 G(M)=PAR(IFUNC)
C  Convert from log weights to weights
      MTOT=M
      WSUM=0.0D0
      DO 60 M=1,MTOT
      WT(M)=DEXP(WT(M)-WTMAX)
      NV1(M)=M
   60 WSUM=WSUM+WT(M)
C  Sort weights and values in ascending order of values; get range points
      CALL DSVRGP(MTOT,G,G,NV1)
      B1=0.25D0*WSUM
      B2=0.75D0*WSUM
      R1=A(1)
      R2=A(2)
      IF(KRANGE.EQ.1)GO TO 70
      R1=A(1)*WSUM
      R2=A(2)*WSUM
   70 WRUN=0.0D0
      DO 80 M=1,MTOT
      W(M)=WT(NV1(M))
      WRUN1=WRUN+W(M)
      IF((WRUN.LT.B1).AND.(WRUN1.GE.B1))Q1=G(M)
      IF((WRUN.LT.B2).AND.(WRUN1.GE.B2))Q2=G(M)
      IF(KRANGE.EQ.1)GO TO 75
      IF((WRUN.LT.R1).AND.(WRUN1.GE.R1))R1=G(M)
      IF((WRUN.LT.R2).AND.(WRUN1.GE.R2))R2=G(M)
   75 WRUN=WRUN1
   80 W(M)=W(M)/WSUM
C  Compute density at ordinates and write to file
      H=HFRAC*(Q2-Q1)
      IF(KRANGE.EQ.2)WRITE(6,85)R1,R2
   85 FORMAT(/,31X,'Ordinate range:',2(1PE15.5))
      WRITE(6,86)H
   86 FORMAT(33X,'Window width:',1PE15.5)
      H1=1.0D0/H
      MBASE=1
      Z=R1
      ZINC=(R2-R1)/(NPTS-1)
      AK=0.5D0
      DO 130 IPTS=1,NPTS
      ZLOWER=Z-H
      ZUPPER=Z+H
      DEN=0.0D0
   90 M=MBASE
      IF(G(M).GT.ZLOWER)GO TO 100
      MBASE=MBASE+1
      GO TO 90
  100 IF(G(M).GE.ZUPPER)GO TO 110
      IF(KERN.EQ.2)AK=1.0D0-DABS(H1*(G(M)-Z))
      IF(KERN.EQ.3)AK=CON3*(1.0D0-(H1*(G(M)-Z))**2)**2
      DEN=DEN+W(M)*AK
      M=M+1
      IF(M.GT.MTOT)GO TO 110
      GO TO 100
  110 WRITE(2,120)Z,H1*DEN
  120 FORMAT(2(1PE20.8))
  130 Z=Z+ZINC
      END