PROGRAM PGEX4
C-----------------------------------------------------------------------
C Demonstration program for PGPLOT: draw histograms.
C-----------------------------------------------------------------------
      REAL PI
      PARAMETER (PI=3.14159265359)
      INTEGER  I, ISEED, PGOPEN
      REAL     DATA(1000), X(620), Y(620)
      REAL     PGRNRM
C
C Call PGOPEN to initiate PGPLOT and open the output device; PGOPEN
C will prompt the user to supply the device name and type. Always
C check the return code from PGOPEN.
C
      IF (PGOPEN('?') .LE. 0) STOP
C
C Call PGRNRM to obtain 1000 samples from a normal distribution.
C
      ISEED = -5678921
      DO I=1,1000
          DATA(I) = PGRNRM(ISEED)
      ENDDO
C
C Draw a histogram of these values.
C
      CALL PGSAVE
      CALL PGHIST(1000,DATA,-3.1,3.1,31,0)
C
C Samples from another normal distribution.
C
      DO I=1,200
          DATA(I) = 1.0+0.5*PGRNRM(ISEED)
      ENDDO
C
C Draw another histogram (filled) on same axes.
C
      CALL PGSCI(15)
      CALL PGHIST(200,DATA,-3.1,3.1,31,3)
      CALL PGSCI(0)
      CALL PGHIST(200,DATA,-3.1,3.1,31,1)
      CALL PGSCI(1)
C
C Redraw the box which may have been clobbered by the histogram.
C
      CALL PGBOX('BST', 0.0, 0, ' ', 0.0, 0)
C
C Label the plot.
C
      CALL PGLAB('Variate', ' ',
     $             'PGPLOT Example 4:  Histograms (Gaussian)')
C
C Superimpose the theoretical distribution.
C
      DO I=1,620
          X(I) = -3.1 + 0.01*(I-1)
          Y(I) = 0.2*1000./SQRT(2.0*PI)*EXP(-0.5*X(I)*X(I))
      ENDDO
      CALL PGLINE(620,X,Y)
      CALL PGUNSA
C
C Finally, call PGCLOS to terminate things properly.
C
      CALL PGCLOS
C
      END


C=======================================================================


      REAL FUNCTION PGRNRM (ISEED)
      INTEGER ISEED
C-----------------------------------------------------------------------
C Returns a normally distributed deviate with zero mean and unit 
C variance. The routine uses the Box-Muller transformation of uniform
C deviates. For a more efficient implementation of this algorithm,
C see Press et al., Numerical Recipes, Sec. 7.2.
C
C Arguments:
C  ISEED  (in/out) : seed used for PGRAND random-number generator.
C
C Subroutines required:
C  PGRAND -- return a uniform random deviate between 0 and 1.
C
C History:
C  1995 Dec 12 - TJP.
C-----------------------------------------------------------------------
      REAL R, X, Y, PGRAND
C
 10   X = 2.0*PGRAND(ISEED) - 1.0
      Y = 2.0*PGRAND(ISEED) - 1.0
      R = X**2 + Y**2
      IF (R.GE.1.0) GOTO 10
      PGRNRM = X*SQRT(-2.0*LOG(R)/R)
C-----------------------------------------------------------------------
      END


C=======================================================================


      REAL FUNCTION PGRAND(ISEED)
      INTEGER ISEED
C-----------------------------------------------------------------------
C Returns a uniform random deviate between 0.0 and 1.0.
C
C NOTE: this is not a good random-number generator; it is only
C intended for exercising the PGPLOT routines.
C
C Based on: Park and Miller's "Minimal Standard" random number
C   generator (Comm. ACM, 31, 1192, 1988)
C
C Arguments:
C  ISEED  (in/out) : seed.
C-----------------------------------------------------------------------
      INTEGER   IM, IA, IQ, IR
      PARAMETER (IM=2147483647)
      PARAMETER (IA=16807, IQ=127773, IR= 2836)
      REAL      AM
      PARAMETER (AM=128.0/IM)
      INTEGER   K
C-
      K = ISEED/IQ
      ISEED = IA*(ISEED-K*IQ) - IR*K
      IF (ISEED.LT.0) ISEED = ISEED+IM
      PGRAND = AM*(ISEED/128)
      RETURN
      END