C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*******************************
C**     PROGRAM 'AFSTATS'     **
C*******************************
C
C
C PURPOSE: Statistics on .afout file as function of wavelength
C
C AUTHOR:  Marjorie Harding
C
C
C====== Parameter statements
C
 
C
C====== GLOBALS:
C
      REAL WVMIN
      REAL WVMAX
      INTEGER NUMBINS
      INTEGER IFW
      COMMON /VALUES/ WVMIN,WVMAX,NUMBINS,IFW
C
C====== FUNCTION TYPES:
C
      LOGICAL CCPEXS 
C
C====== LOCALS:
C
      REAL BIN(200)
      INTEGER NBIN(200)
      REAL WAVE(200)
      REAL FAV(200)
      REAL RBIN(200)
      CHARACTER*120 DATIN
      CHARACTER*120 DATOUT
 
C
C-------------------------------------------------------------------------------
C
C
C====== Get input and output file names and open files
C
      WRITE (6,6010)
5     READ (5,6020) DATIN
      IF (DATIN.EQ.' ') STOP
      IF (.NOT.CCPEXS(DATIN)) THEN
         WRITE (6,6080)
         GO TO 5
      END IF
      IFAIL = 0
      CALL CCPDPN(-8,DATIN,'READONLY','F',0,IFAIL)
10    WRITE (6,6030)
      READ (5,6020) DATOUT
      IF (DATOUT.EQ.' ') STOP
      IFAIL = 0
      CALL CCPDPN(-9,DATOUT,'NEW','F',0,IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (6,6090)
         GO TO 10
      END IF
C
C====== Get required option and parameters
C
      WRITE (6,6040)
      READ (5,*) IFW
      IF (IFW.EQ.1) WRITE (6,6050)
      IF (IFW.EQ.2) WRITE (6,6060)
      WRITE (9,6070) DATIN(1:LENSTR(DATIN))
      READ (5,*) WVMIN , WVMAX , NUMBINS
      WVMIN = WVMIN - (WVMAX-WVMIN)/NUMBINS
      NUMBINS = NUMBINS + 1
C
C====== Accumulate and display/print statistics
C
      CALL STATS(BIN,NBIN,WAVE,FAV,RBIN)
      CLOSE (8)
      CLOSE (9)
      STOP
c
C====== Format statements
C
 6010 FORMAT (' Enter name of input file: ',$)
 6020 FORMAT (A)
 6030 FORMAT (' Name of output file: ',$)
 6040 FORMAT (' Statistics on wavelength (1) or d-spacing (2): ',$)
 6050 FORMAT (' Enter wvmin, wvmax and numbins: ',$)
 6060 FORMAT (' Enter dmin, dmax, and numbins: ',$)
 6070 FORMAT ('Intensity statistics,  input data from file: ',A,//)
 6080 FORMAT (' **Error** file does not exist')
 6090 FORMAT (' **Error** Unable to open output file')
      END
C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************
C**     STATS       **
C*********************
C
C
      SUBROUTINE STATS(BIN, NBIN, WAVE, FAV, RBIN)  
C
C
C PURPOSE:
C
C AUTHOR:  Marjorie Harding
C
C ARGUMENTS:
C
      REAL BIN(*)
      INTEGER NBIN(*)
      REAL WAVE(*)
      REAL FAV(*)
      REAL RBIN(*)
C
C====== Parameter statements
C
 
C
C====== GLOBALS:
C
      REAL WVMIN
      REAL WVMAX
      INTEGER NUMBINS
      INTEGER IFW
      COMMON /VALUES/ WVMIN,WVMAX,NUMBINS,IFW
C
C====== FUNCTION TYPES:
C
 
C
C====== LOCALS:
C
      CHARACTER*4 ANS
      CHARACTER*80 STR
C
C-------------------------------------------------------------------------------
C
      WRITE (6,6010)
      READ (5,*) SIGTEST
C
C====== Clear arrays
C
      DO 20 I = 1 , NUMBINS
         BIN(I) = 0
         NBIN(I) = 0
 20   CONTINUE
      NOBS = 0
      NGTSIG = 0
C
C====== Output to log file
C
      WRITE (9,'(1X)')
      STR = 'Using reflections with I > '
      CALL STR_ADDF (STR,SIGTEST,1,1)
      CALL STR_ADDS (STR,'sigma(I)',1)
      WRITE (9,'(A)') STR
      WRITE (9,'(1X)')
      WVBIN = (WVMAX-WVMIN)/NUMBINS
C
C====== Read data from .afout file
C
 30   READ (8,*,END=40) IH , K , L , WVL , THETA , FF , SIGFF , XF , YF
      IF (FF.GT.SIGFF*SIGTEST) THEN
C
C====== If d-spacing stats required set wvl = lambda / 2sintheta
C
         IF (IFW.EQ.2) THEN
            D = WVL/(2*SIN(THETA))
            WVL = D
         END IF

         IF (WVL.LT.WVMIN) GO TO 30
         IF (WVL.GT.WVMAX) GO TO 30
         I = IFIX((WVL-WVMIN)/WVBIN)
         BIN(I) = BIN(I) + FF
         NBIN(I) = NBIN(I) + 1
         NGTSIG = NGTSIG + 1
      END IF
      NOBS = NOBS + 1
      GO TO 30

 40   STR= ' '
      CALL STR_ADDI (STR,NOBS,0)
      CALL STR_ADDS (STR,'reflections observed',1)
      WRITE (6,'(A)') STR
      STR =' '
      CALL STR_ADDI (STR,NGTSIG,0)
      CALL STR_ADDS (STR,'reflections with I > sigtest',1)
      WRITE (6,'(A)') STR
      STR =' '
      CALL STR_ADDI (STR,NOBS,0)
      CALL STR_ADDS (STR,'reflections read in from file',1)
      WRITE (9,'(A)') STR
      STR =' '
      CALL STR_ADDI (STR,NGTSIG,0)
      CALL STR_ADDS (STR,'reflections with I >',1)
      CALL STR_ADDF (STR,SIGTEST,1,1)
      CALL STR_ADDS (STR,'* sig(I)',1)
      WRITE (9,'(A)') STR
      WRITE (9,'(1X)')
      IF (IFW.EQ.1) WRITE (9,'(A)') ' wavelength     number   Imean'
      IF (IFW.EQ.2) WRITE (9,'(A)') ' d-spacing      number   Imean'
      DO 50 I = 1 , NUMBINS
         J = NBIN(I)
         IF (J.EQ.0) J = 1
         FMEAN = BIN(I)/J
         IF (NBIN(I).EQ.0) FMEAN = 0
         WVL = WVMIN + WVBIN*I
         WAVE(I) = WVL
         RBIN(I) = FLOAT(NBIN(I))
         FAV(I) = FMEAN
         WRITE (9,6020) WVL , (-WVL-WVBIN) , RBIN(I) , FAV(I)
 50   CONTINUE

      IBIN = 0
      IFMEAN = 0
      IHISTO = 0
      IGRAPH = 0

 60   WRITE (6,6030)
      READ (5,6040) ANS
      IF (ANS.EQ.'n' .OR. ANS.EQ.'N') IBIN = 1
      IF (ANS.EQ.'f' .OR. ANS.EQ.'F') IFMEAN = 1
      IF (IBIN.EQ.0.AND.IFMEAN.EQ.0) GO TO 60

      WRITE (6,6050)
      READ (5,6040) ANS
      IF (ANS.EQ.'g' .OR. ANS.EQ.'G') IGRAPH = 1
      IF (ANS.EQ.'h' .OR. ANS.EQ.'H') IHISTO = 1


      IF (IGRAPH.EQ.1) THEN
         IF (IBIN.EQ.1) CALL GRAPH(WAVE,RBIN,1,NUMBINS)
         IF (IFMEAN.EQ.1) CALL GRAPH(WAVE,FAV,1,NUMBINS)
      ELSE
         IF (IBIN.EQ.1) CALL HISTO(WAVE,RBIN,1,NUMBINS,NUMBINS)
         IF (IFMEAN.EQ.1) CALL HISTO(WAVE,FAV,1,NUMBINS,NUMBINS)
      END IF

      RETURN
C
C====== Format statements
C
 6010 FORMAT (' Use I > const * sigma(I)  ;  give const: ',$)
 6020 FORMAT (2F6.3,F10.0,F10.2)
 6030 FORMAT (' Output nbin (n) or FFmean (f): ',$)
 6040 FORMAT (A)
 6050 FORMAT (' Output graph (g) or histogram (h): ',$)
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C********************
C**     HISTO      **
C********************
C
C
      SUBROUTINE HISTO (X, Y, NSTART, NEND, NUMBINS)
C
C
C PURPOSE: Graphical routine using qscreen graphics library
C
C AUTHOR:  Paul Carr (qsceen_library  John Campbell)
C
C
C ARGUMENTS:
C
      REAL X(NSTART:NEND)
      REAL Y(NSTART:NEND)
      INTEGER NSTART
      INTEGER NEND
      INTEGER NUMBINS
C
C====== Parameter statements
C 

C
C====== GLOBALS:
C

C
C====== FUNCTION TYPES:
C

C
C====== LOCALS:
C
      CHARACTER*1 IYN
      CHARACTER*20 TITLE
      CHARACTER*20 X_LABEL
      CHARACTER*20 Y_LABEL
      CHARACTER*121 FNAME
      CHARACTER*121 STRING
      INTEGER TEMP_UNIT
      INTEGER IN_UNIT
      INTEGER OUT_UNIT
      REAL MINX
      REAL MAXX
      REAL MINY
      REALMAXY
      REAL MM_TO_PIXEL
C
C-------------------------------------------------------------------------------
C
C
C====== Initialise graphics
C
      TEMP_UNIT = 9
      IN_UNIT = 5
      OUT_UNIT = 6
      IXTERM = -1
      OPEN (TEMP_UNIT,FILE='QGTEMP',STATUS='UNKNOWN')
20    WRITE (OUT_UNIT,6000)
      READ (IN_UNIT,6010) IYN
      IF (IYN.EQ.'Y' .OR. IYN.EQ.'y') IXTERM = 1
      IF (IYN.EQ.'N' .OR. IYN.EQ.'n') IXTERM = 0
      IF (IXTERM.NE.1 .AND. IXTERM.NE.0) GO TO 20

      CALL QGINIT(TEMP_UNIT,OUT_UNIT,IXTERM,IN_UNIT,OUT_UNIT)
C
C====== Get max and min values of x and y
C
      CALL GETMAX(X,Y,NSTART,NEND,MINX,MAXX,MINY,MAXY,NVALS)
C
C====== Dialogue with user for titles,labels etc
C
 
      IJOIN = 0
      IPOST = 0
      WRITE (OUT_UNIT,6040)
      READ (IN_UNIT,6010) TITLE
      WRITE (OUT_UNIT,6050)
      READ (IN_UNIT,6010) X_LABEL
      WRITE (OUT_UNIT,6060)
      READ (IN_UNIT,6010) Y_LABEL
C
C====== Draw axes and labels
C
      CHSIZE = 20
      NTICK = 5
      NDECX = 3
      NDECY = 2
      MM_TO_PIXEL = 3.5
      CALL QGENTR
      CALL QGPSET(2)
      CALL QGTEXT(75.0*MM_TO_PIXEL,210.0*MM_TO_PIXEL,1.5*CHSIZE,TITLE)
      CALL QGTEXT(230*MM_TO_PIXEL,15*MM_TO_PIXEL,CHSIZE,X_LABEL)
      CALL QGTEXT(10*MM_TO_PIXEL,190*MM_TO_PIXEL,CHSIZE,Y_LABEL)
      CALL QGXAX1(20*MM_TO_PIXEL,35*MM_TO_PIXEL,250*MM_TO_PIXEL,MINX,
     +            MAXX,NTICK,NDECX,CHSIZE/2,CHSIZE)
      CALL QGYAX1(20*MM_TO_PIXEL,35*MM_TO_PIXEL,150*MM_TO_PIXEL,MINY,
     +            MAXY,NTICK,NDECY,CHSIZE/2,CHSIZE)
C
C====== Map graph area to X and Y space
C
      CALL QGMAPC(20*MM_TO_PIXEL,270*MM_TO_PIXEL,35*MM_TO_PIXEL,
     +            185*MM_TO_PIXEL,MINX,MAXX,MINY,MAXY)
C
C====== Draw historgram
C
      XWIDTH = (X(NEND)-X(NSTART))/NUMBINS

      CALL QGMOVE(MINX,MINY)
      DO 30 I = 1 , NUMBINS
         CALL QGMOVE(X(I),MINY)
         CALL QGDRAW(X(I),Y(I))
         CALL QGDRAW(X(I)+XWIDTH,Y(I))
         CALL QGDRAW(X(I)+XWIDTH,MINY)
 30   CONTINUE
C 
C====== Check if hard copy required
C
      CALL QGPSET(0)
      CALL QGPIXC
      CALL QGMOVE(20*MM_TO_PIXEL,15*MM_TO_PIXEL)
      WRITE (STRING,6020)
      CALL QGTEXT(20*MM_TO_PIXEL,15*MM_TO_PIXEL,CHSIZE,STRING)
      READ (IN_UNIT,6010) IYN
      IF (IYN.EQ.'Y' .OR. IYN.EQ.'y') IPOST = 1
      IF (IPOST.EQ.1) THEN
         WRITE (STRING,6030)
         CALL QGTEXT(20*MM_TO_PIXEL,10*MM_TO_PIXEL,CHSIZE,STRING)
         READ (IN_UNIT,6010) FNAME
         CALL QGTEXT(20*MM_TO_PIXEL,5*MM_TO_PIXEL,CHSIZE,FNAME)
         OPEN (10,FILE=FNAME,STATUS='new')
         CALL QGWRPS(10,1)
      END IF
      CALL QGALPH
      CALL QGCLOS
      RETURN
C
C====== Format statements
C
 6000 FORMAT (' Is this a unix xterm window (y/n): ',$)
 6010 FORMAT (A)
 6020 FORMAT (' Do you want a hard copy plot (y/n): ',$)
 6030 FORMAT (' Name of output postscript file: ',$)
 6040 FORMAT (' Input title for plot(max 20 chars): ',$)
 6050 FORMAT (' Input label for x-axis (max 20 chars): ',$)
 6060 FORMAT (' Input label for y-axis(max 20 chars): ',$)
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************
C**     GETMAX      **
C*********************
C
C
      SUBROUTINE GETMAX (X, Y, NSTART, NEND, MINX, MAXX, MINY, MAXY,
     +                   NVALS)
C
C
C PURPOSE: Find max and min of two arrays
C
C AUTHOR:  Marjorie Harding
C
C ARGUMENTS:
C
      REAL X(NSTART:NEND)
      REAL Y(NSTART:NEND)
      INTEGER NSTART
      INTEGER NEND
      REAL MINX
      REAL MAXX
      REAL MINY
      REAL MAXY
      INTEGER NVALS
C
C====== Parameter statements
C 

C
C====== GLOBALS:
C

C
C====== FUNCTION TYPES:
C

C
C====== LOCALS:
C

C
C-------------------------------------------------------------------------------
C
C
C====== Initialisations
C
      MINX = 99999.0
      MAXX = -99999.0
      MINY = 99999.9
      MAXY = -99999.9
C
C====== Loop through values
C
      DO 20 I = NSTART , NEND
         IF (X(I).LT.MINX) MINX = X(I)
         IF (X(I).GT.MAXX) MAXX = X(I)
         IF (Y(I).LT.MINY) MINY = Y(I)
         IF (Y(I).GT.MAXY) MAXY = Y(I)
 20   CONTINUE
      NVALS = NEND - NSTART + 1
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************
C**     GRAPH       **
C*********************
C
C
      SUBROUTINE GRAPH (X, Y, NSTART, NEND)
C
C
C PURPOSE: Graphical routine using qscreen graphics library
C
C AUTHOR:  Paul Carr (using qscreen_graphics John Campbell)
C
C ARGUMENTS:
C
      REAL X(NSTART:NEND)
      REAL Y(NSTART:NEND)
      INTEGER NSTART
      INTEGER NEND
C
C====== Parameter statements
C 

C
C====== GLOBALS:
C

C
C====== FUNCTION TYPES:
C

C
C====== LOCALS:
C
      CHARACTER*1 IYN
      CHARACTER*20 TITLE
      CHARACTER*20 X_LABEL
      CHARACTER*20 Y_LABEL
      CHARACTER*12 FNAME
      CHARACTER*121 STRING
      INTEGER TEMP_UNIT, IN_UNIT, OUT_UNIT
      REAL MINX,MAXX,MINY,MAXY,MM_TO_PIXEL

C
C-------------------------------------------------------------------------------
C
C
C====== Initialise graphics
C
      TEMP_UNIT = 9
      IN_UNIT = 5
      OUT_UNIT = 6
      IXTERM = -1
      OPEN (TEMP_UNIT,FILE='QGTEMP',STATUS='UNKNOWN')
20    WRITE (OUT_UNIT,6000)
      READ (IN_UNIT,6010) IYN
      IF (IYN.EQ.'Y' .OR. IYN.EQ.'y') IXTERM = 1
      IF (IYN.EQ.'N' .OR. IYN.EQ.'n') IXTERM = 0
      IF (IXTERM.NE.1 .AND. IXTERM.NE.0) GO TO 20
      CALL QGINIT(TEMP_UNIT,OUT_UNIT,IXTERM,IN_UNIT,OUT_UNIT)
C
C====== Get max and min values of x and y
C
      CALL GETMAX(X,Y,NSTART,NEND,MINX,MAXX,MINY,MAXY,NVALS)
C
C====== Dialogue with user for titles,labels etc
C
 
      IJOIN = 0
      IPOST = 0
      WRITE (OUT_UNIT,6050)
      READ (IN_UNIT,6010) TITLE
      WRITE (OUT_UNIT,6060)
      READ (IN_UNIT,6010) X_LABEL
      WRITE (OUT_UNIT,6070)
      READ (IN_UNIT,6010) Y_LABEL
 
      WRITE (OUT_UNIT,6020) NVALS
      READ (IN_UNIT,6010) IYN
      IF (IYN.EQ.'Y' .OR. IYN.EQ.'y') IJOIN = 1
C
C====== Draw axes and labels
C
      CHSIZE = 20
      NTICK = 5
      NDECX = 4
      NDECY = 1
      MM_TO_PIXEL = 3.5
      CALL QGENTR
      CALL QGPSET(2)
      CALL QGTEXT(75.0*MM_TO_PIXEL,210.0*MM_TO_PIXEL,1.5*CHSIZE,TITLE)
      CALL QGTEXT(230*MM_TO_PIXEL,15*MM_TO_PIXEL,CHSIZE,X_LABEL)
      CALL QGTEXT(10*MM_TO_PIXEL,190*MM_TO_PIXEL,CHSIZE,Y_LABEL)
      CALL QGXAX1(20*MM_TO_PIXEL,35*MM_TO_PIXEL,250*MM_TO_PIXEL,MINX,
     +            MAXX,NTICK,NDECX,CHSIZE/2,CHSIZE)
      CALL QGYAX1(20*MM_TO_PIXEL,35*MM_TO_PIXEL,150*MM_TO_PIXEL,MINY,
     +            MAXY,NTICK,NDECY,CHSIZE/2,CHSIZE)
C
C====== Map graph area to X and Y space
C
      CALL QGMAPC(20*MM_TO_PIXEL,270*MM_TO_PIXEL,35*MM_TO_PIXEL,
     +            185*MM_TO_PIXEL,MINX,MAXX,MINY,MAXY)
C
C====== Set symbol type and plot points
C
      ISYMB = 2
      SYMBOL_SIZE = 3
C
C====== 2=diagonal cross, 3=square
C
      CALL QGMOVE(X(NSTART),Y(NSTART))
      DO 30 I = NSTART , NEND - 1
         IF (IJOIN.EQ.1 .AND. I.GT.NSTART) CALL QGDRAW(X(I),Y(I))
         CALL QGSYMB(X(I),Y(I),ISYMB,SYMBOL_SIZE)
         CALL QGMOVE(X(I),Y(I))
 30   CONTINUE
C
C====== Check if hard copy required
C
      CALL QGPSET(0)
      CALL QGPIXC
      CALL QGMOVE(20*MM_TO_PIXEL,15*MM_TO_PIXEL)
      WRITE (STRING,6030)
      CALL QGTEXT(20*MM_TO_PIXEL,15*MM_TO_PIXEL,CHSIZE,STRING)
      READ (IN_UNIT,6010) IYN
      IF (IYN.EQ.'Y' .OR. IYN.EQ.'y') IPOST = 1
      IF (IPOST.EQ.1) THEN
         WRITE (STRING,6040)
         CALL QGTEXT(20*MM_TO_PIXEL,10*MM_TO_PIXEL,CHSIZE,STRING)
         READ (IN_UNIT,6010) FNAME
         CALL QGTEXT(20*MM_TO_PIXEL,5*MM_TO_PIXEL,CHSIZE,FNAME)
         OPEN (10,FILE=FNAME,STATUS='new')
         CALL QGWRPS(10,1)
      END IF
      CALL QGALPH
      CALL QGCLOS
      RETURN
C
C====== Format statements
C
 6000 FORMAT (' Is this a unix xterm window (y/n): ',$)
 6010 FORMAT (A)
 6020 FORMAT (1X,I6,' points to plot, do you want to "join',
     +        ' the dots" (y/n): ',$)
 6030 FORMAT (' Do you want a hard copy plot (y/n): ',$)
 6040 FORMAT (' Name of output postscript file: ',$)
 6050 FORMAT (' Input title for plot(max 20 chars): ',$)
 6060 FORMAT (' Input label for x-axis (max 20 chars): ',$)
 6070 FORMAT (' Input label for y-axis(max 20 chars): ',$)
      END