      PROGRAM SAD
C
C VERSION
C     14-MAY-12  AD Allow for 14 boat divisions in Eights
C     26-JAN-05  AD Add rise.dat output
C     20-AUG-02  AD Remove non-standard FORTRAN features.
C     21-MAY-01  AD Add SPDFAC
C     21-APR-99  AD 
C
C DESCRIPTION
C     Statistical Analysis of Divisions.
C
      IMPLICIT NONE
C
      REAL RANGAU
C
C LOCAL CONSTANTS
      INTEGER       NPRB
        PARAMETER ( NPRB = 8 )
      INTEGER       NRUN
        PARAMETER ( NRUN = 10000 )
      INTEGER       MAXBNG   ! Max no bunglines in one divisions
        PARAMETER ( MAXBNG = 14 ) 
      INTEGER       MAXDIV
        PARAMETER ( MAXDIV = 4 )
      INTEGER       MAXCRW
        PARAMETER ( MAXCRW = (MAXBNG-1)*MAXDIV + 1 )
      REAL          SPDFAC   ! Difference in SAD scores equiv to 50% bump prob.
        PARAMETER ( SPDFAC = 0.4 )
C
C LOCAL VARIABLES
      INTEGER      IBNG, INEW, NBNG
      INTEGER      ICRW, JCRW, NCRW
      INTEGER      IDIV, NDIV
      INTEGER      IPOS, IRUN, ISAVE
      INTEGER      ISEED
      INTEGER      IDAY, NDAY
      INTEGER      CRWPOS(MAXCRW), POSCRW(MAXCRW)
      INTEGER      IPRB, PRBMAX
      INTEGER      LBNG
      INTEGER      BMPBNG, BPDBNG
      INTEGER      BMPPOS, BPDPOS
      INTEGER      NEWBNG(MAXBNG), NEWCRW(MAXBNG)
      INTEGER      NOUT
      REAL         AVGBMP(MAXCRW)
      REAL         AVGSPD(MAXCRW), DAYSPD(MAXBNG)
      REAL         DAYSD(MAXCRW), RUNSD(MAXCRW)
      REAL         RUNSPD(MAXCRW), BMPSPD(MAXCRW)
      REAL         SPDMAX
      REAL         RESULT(MAXCRW,MAXCRW)
      REAL         RMAX
      REAL         SUM
      LOGICAL      CREWIN(MAXBNG)
      LOGICAL      EIGHTS
      LOGICAL      MENDIV
      CHARACTER*2  MWSTR
      CHARACTER*8  STRING(NPRB)
      CHARACTER*80 INPFIL
      CHARACTER*80 DATFIL
      CHARACTER*80 HTMFIL
      CHARACTER*20 CREW(MAXCRW)
      CHARACTER*20 CRWFUL
C
C EXECUTABLE CODE -------------------------------------------------------------
C
      EIGHTS = .TRUE.
      IF ( EIGHTS ) THEN
        WRITE ( *, * ) 'Eights Rules'
        NBNG = 13
      ELSE
        WRITE ( *, * ) 'Torpids Rules'
        NBNG = 12
      END IF
      NOUT = 3*NBNG + 1           ! no.crews output to data files 37 or 40
c
      write ( *, * ) 'Spd Fac = ', SPDFAC

C Ask user for the number of divisions
C
      WRITE ( *, '($X,A)' ) 'No. Divisions: '
      READ ( *, * ) NDIV
      IF ( NDIV .GT. MAXDIV ) STOP 'NDIV .GT. MAXDIV'
C
C Ask for data file and read in crew stats.
C Last division should contain extra boat.
C
      NCRW = 0
      WRITE ( *, '($X,A)' ) 'Input data file: '
      READ ( *, '(A)' ) INPFIL
      IF ( INDEX ( INPFIL, 'sadm' ) .GT. 0 ) THEN
        MENDIV = .TRUE.
        WRITE ( *, * ) 'Assuming Men''s Divisions'
      ELSE
        MENDIV = .FALSE.
        WRITE ( *, * ) 'Assuming Women''s Divisions'
      END IF
      OPEN ( UNIT=1, FILE=INPFIL, STATUS='OLD' ) 
      DO IDIV = 1, NDIV
        DO ICRW = 1, NBNG
          NCRW = NCRW + 1
          READ ( 1, * ) 
     &      CREW(NCRW), AVGSPD(NCRW), RUNSD(NCRW), DAYSD(NCRW)
        END DO
        IF ( IDIV .EQ. NDIV ) THEN
          NCRW = NCRW + 1
          READ ( 1, * ) 
     &      CREW(NCRW), AVGSPD(NCRW), RUNSD(NCRW), DAYSD(NCRW)
        END IF           
      END DO
C
C Ask user for number of days of bumps racing to be predicted.
C
      WRITE ( *, '($X,A)' ) 'Number of Days: '
      READ ( *, * ) NDAY
C
      do icrw = 1, ncrw
        do jcrw = 1, ncrw
          result(jcrw,icrw) = 0.0
        enddo
      enddo
C
      ISEED = 112233 
      DO IRUN = 1, NRUN
C
C Initialise starting order at start of each new run, using random number *RUNSD
C
        DO ICRW = 1, NCRW
          POSCRW(ICRW) = ICRW
          CRWPOS(ICRW) = ICRW
          RUNSPD(ICRW) = AVGSPD(ICRW) + RANGAU(ISEED) * RUNSD(ICRW)
        END DO
C
C Loop over each requested day and division, starting with lowest division
C Copy finishing order from previous day to starting order
C Calculate new speed for each day using random number *DAYSD
C
        DO IDAY = 1, NDAY
          DO IDIV = NDIV, 1, -1
            DO IBNG = 1, NBNG+1
              IPOS = ( IDIV - 1 ) * NBNG + IBNG
              CREWIN(IBNG) = .TRUE.
              NEWBNG(IBNG) = IBNG
              ICRW = CRWPOS(IPOS)
              DAYSPD(IBNG) = RUNSPD(ICRW) + RANGAU(ISEED) * DAYSD(ICRW)
            END DO 
C
C Find quickest bump, swap crews around, remove from start order and repeat
C until no more bumps occur in division.
C
  200       CONTINUE
            SPDMAX = -100.0
C
C Find bungline number of top crew still racing (LBNG=1 always, in Torpids)
            LBNG = 1
            DO WHILE ( .NOT. CREWIN(LBNG) .AND. LBNG .LE. NBNG )
              LBNG = LBNG + 1
              IF ( .NOT. EIGHTS ) STOP 'F-SAD: Logical Error#1'
            END DO
C No. crews is NBNG+1 so if NBNG is even there should always be at least one
C crew still racing
            IF ( MOD ( NBNG, 2 ) .EQ. 0 .AND. LBNG .EQ. NBNG+2 ) 
     &        STOP 'F-SAD: Logical Error#2'
C With NBNG odd, it is possible that all crews have bumped out (LBNG=NBNG+2)
C            
            DO IBNG = LBNG+1, NBNG+1
              IF ( CREWIN(IBNG) ) THEN
                BMPSPD(IBNG) = 
     &            1.5*(DAYSPD(IBNG)-DAYSPD(LBNG)) / SPDFAC
     &            / ( FLOAT(IBNG-LBNG)*2.5 - 1.0 ) 
                IF ( BMPSPD(IBNG) .GT. SPDMAX ) THEN
                  SPDMAX = BMPSPD(IBNG)
                  BMPBNG = IBNG
                  BPDBNG = LBNG
                END IF
                LBNG = IBNG
              END IF
            END DO
C
C Found next bump, so swap crews around in finishing order and remove.
C 
            IF ( SPDMAX .GT. 1.0 ) THEN
              BMPPOS = ( IDIV - 1 ) * NBNG + BMPBNG
              BPDPOS = ( IDIV - 1 ) * NBNG + BPDBNG
              CREWIN(BMPBNG) = .FALSE.
              IF ( EIGHTS ) THEN
                CREWIN(BPDBNG) = .FALSE.
                NEWBNG(BPDBNG) = BMPBNG
                NEWBNG(BMPBNG) = BPDBNG
              ELSE 
                DO IBNG = NEWBNG(BPDBNG)+1, NEWBNG(BMPBNG)
                  NEWBNG(IBNG) = NEWBNG(IBNG) - 1
                  NEWBNG(BPDBNG) = NEWBNG(BPDBNG) + 1
                END DO
              END IF
              GOTO 200
            END IF
C
C End of division: put each crew into new finishing position
C
            DO IBNG = 1, NBNG+1        ! Loop over original bungline#
              INEW = NEWBNG(IBNG)  ! New bungline# 
              IPOS = ( IDIV - 1 ) * NBNG + IBNG
              ICRW = CRWPOS(IPOS)  ! Crew# on original bungline 
              NEWCRW(INEW) = ICRW  ! Crew# on new bungline
              POSCRW(ICRW) = IPOS + INEW - IBNG ! New pos of orig.crew
            END DO
            DO IBNG = 1, NBNG+1
              IPOS = ( IDIV - 1 ) * NBNG + IBNG
              CRWPOS(IPOS) = NEWCRW(IBNG)  ! Crew# on new bunglines
            END DO
          END DO
        END DO
C
C At the end of each run note crew's finishing position
C
        DO ICRW = 1, NCRW
          DO IPOS = 1, NCRW
            IF ( CRWPOS(IPOS) .EQ. ICRW ) THEN
              RESULT(ICRW,IPOS) = RESULT(ICRW,IPOS) + 1.0
              GOTO 300
            END IF
          END DO
  300     CONTINUE
        END DO
      END DO
C
C Work out mean rise/fall of each crew
      DO ICRW = 1, NCRW
        SUM = 0.0
        DO IPOS = 1, NCRW
          SUM = SUM + RESULT(ICRW,IPOS)*(ICRW-IPOS)
        END DO
        AVGBMP(ICRW) = SUM / NRUN
      END DO
C      
C Print out results 
C
      IF ( MENDIV ) THEN
        MWSTR = '_m'
      ELSE 
        MWSTR = '_w'
      END IF 
      DATFIL = 'sad_pred' // MWSTR // '.dat'
      HTMFIL = 'sad_pred' // MWSTR // '.html'
      WRITE ( *, * ) 'Writing results to file ', DATFIL(1:30)
      OPEN ( UNIT=2, FILE=DATFIL, STATUS='UNKNOWN' ) 
      OPEN ( UNIT=3, FILE=HTMFIL, STATUS='UNKNOWN' ) 
      WRITE ( 3, * ) '<TABLE>'
      WRITE ( 2, * ) 'Finishing Positions'
      WRITE ( 3, * ) '<TR><TH>Pos <TH>Start Crew<TH> Score <TH> Av.Rise'
      WRITE ( 3, * ) '  <TH COLSPAN=8> Finishing Crew'
      DO IPOS = 1, NCRW
        IF ( MOD ( IPOS, NBNG ) .EQ. 1 ) THEN
          IDIV = 1 + INT(IPOS/NBNG)
          WRITE ( 2, '(X,A,I1)' ) 'Division ', IDIV
          IF ( IDIV .LE. 4 ) 
     &  WRITE ( 3, * ) '<TR><TH ALIGN=LEFT COLSPAN=4> Division', IDIV
          IBNG = 0
        END IF
        IBNG = IBNG + 1
        DO IPRB = 1, NPRB
          STRING(IPRB) = ' '
        END DO
        DO IPRB = 1, NPRB
          PRBMAX = 0
          DO ICRW = 1, NCRW
           IF ( RESULT(ICRW,IPOS) .GT. PRBMAX ) THEN
              JCRW = ICRW
              PRBMAX = RESULT(ICRW,IPOS)
            END IF
          END DO
          IF ( PRBMAX .GT. 0 ) THEN
            WRITE ( STRING(IPRB), '(A4,I3,A)' ) CREW(JCRW)(1:8), 
     &        ( NINT ( 100.0*RESULT(JCRW,IPOS)/FLOAT(NRUN) ) ),'%'
            RESULT(JCRW,IPOS) = -RESULT(JCRW,IPOS)
          ELSE
            GOTO 400
          END IF
        END DO      
  400   CONTINUE
        WRITE ( 2, '(X,I2,a3,a4,F4.1,SPF5.1,A,8(2X,A8))' ) IBNG, '. (', 
     &    CREW(IPOS)(1:4), AVGSPD(IPOS), AVGBMP(IPOS), ')',
     &    (STRING(IPRB),IPRB=1,NPRB)
        IF ( IPOS .LE. NOUT ) THEN
          CALL CRWNAM ( CREW(IPOS)(1:4), CRWFUL )
          WRITE ( 3, '(A,I2,A)' ) 
     &      '<TR><TD ALIGN=RIGHT>', IBNG, '<TD><A HREF=./sad/' //
     &      CREW(IPOS)(1:4) // MWSTR // '.html>' // CRWFUL // '</A>'
          WRITE ( 3, '(A,F4.1,A,SPF5.1)' ) 
     &      '  <TD ALIGN=CENTER>', AVGSPD(IPOS), 
     &        '<TD ALIGN=CENTER>', AVGBMP(IPOS)
          WRITE ( 3, * ) 
     &      ('  <TD ALIGN=RIGHT> ', STRING(IPRB), IPRB = 1, NPRB )
        END IF
      END DO
      CLOSE ( 2 )
      WRITE ( 3, * ) '</TABLE>'
      CLOSE ( 3 )
C
C Sort by average change in position
      DATFIL = 'rise' // MWSTR // '.dat'
      OPEN (UNIT=2,FILE=DATFIL,STATUS='UNKNOWN')
      DO IPOS = 1, min ( NCRW, NOUT )
        RMAX = -100
        DO ICRW = 1, min ( NCRW, NOUT )      ! top 3 divisions only
          IF ( AVGBMP(ICRW) .GT. RMAX ) THEN
            JCRW = ICRW
            RMAX = AVGBMP(ICRW)
          END IF 
        END DO
        WRITE ( 2, '(I2,X,A4,SPF5.1,SI3)' ) 
     &    IPOS, CREW(JCRW)(1:4), AVGBMP(JCRW), JCRW
        AVGBMP(JCRW) = -200
      END DO
C
      END
      SUBROUTINE CRWNAM ( SHORT, LONG )
C
      IMPLICIT NONE
C
      CHARACTER*4 SHORT
      CHARACTER*(*) LONG
C
C EXECUTABLE CODE ------------------------------------------------------------
      LONG = SHORT
      IF ( SHORT .EQ. 'Ball' ) LONG = 'Balliol' 
      IF ( SHORT .EQ. 'Bal2' ) LONG = 'Balliol II' 
      IF ( SHORT .EQ. 'Bras' ) LONG = 'Brasenose' 
      IF ( SHORT .EQ. 'Chri' ) LONG = 'Christ Church'
      IF ( SHORT .EQ. 'ChC2' ) LONG = 'Christ Church II'
      IF ( SHORT .EQ. 'CCh2' ) LONG = 'Christ Church II'
      IF ( SHORT .EQ. 'Corp' ) LONG = 'Corpus Christi' 
      IF ( SHORT .EQ. 'Exet' ) LONG = 'Exeter'
      IF ( SHORT .EQ. 'Exe2' ) LONG = 'Exeter II'
      IF ( SHORT .EQ. 'Hert' ) LONG = 'Hertford' 
      IF ( SHORT .EQ. 'Her2' ) LONG = 'Hertford II' 
      IF ( SHORT .EQ. 'Jesu' ) LONG = 'Jesus' 
      IF ( SHORT .EQ. 'Kebl' ) LONG = 'Keble' 
      IF ( SHORT .EQ. 'Keb2' ) LONG = 'Keble II' 
      IF ( SHORT .EQ. 'Lady' ) LONG = 'L.M.H.' 
      IF ( SHORT .EQ. 'Lina' ) LONG = 'Linacre' 
      IF ( SHORT .EQ. 'Lna2' ) LONG = 'Linacre II' 
      IF ( SHORT .EQ. 'Linc' ) LONG = 'Lincoln' 
      IF ( SHORT .EQ. 'Lnc2' ) LONG = 'Lincoln II' 
      IF ( SHORT .EQ. 'Magd' ) LONG = 'Magdalen' 
      IF ( SHORT .EQ. 'Mag2' ) LONG = 'Magdalen II' 
      IF ( SHORT .EQ. 'Mans' ) LONG = 'Mansfield' 
      IF ( SHORT .EQ. 'Mert' ) LONG = 'Merton'
      IF ( SHORT .EQ. 'Mer2' ) LONG = 'Merton II'
      IF ( SHORT .EQ. 'NewC' ) LONG = 'New College' 
      IF ( SHORT .EQ. 'New2' ) LONG = 'New College II' 
      IF ( SHORT .EQ. 'Orie' ) LONG = 'Oriel'
      IF ( SHORT .EQ. 'Ori2' ) LONG = 'Oriel II'
      IF ( SHORT .EQ. 'Osle' ) LONG = 'Osler House'
      IF ( SHORT .EQ. 'OslG' ) LONG = 'Osler-Green'
      IF ( SHORT .EQ. 'Pemb' ) LONG = 'Pembroke'
      IF ( SHORT .EQ. 'Pem2' ) LONG = 'Pembroke II'
      IF ( SHORT .EQ. 'Quee' ) LONG = 'Queen''s'
      IF ( SHORT .EQ. 'SAnn' ) LONG = 'St Anne''s'
      IF ( SHORT .EQ. 'SAnt' ) LONG = 'St Antony''s'
      IF ( SHORT .EQ. 'SCat' ) LONG = 'St Catherine''s' 
      IF ( SHORT .EQ. 'SCa2' ) LONG = 'St Catherine''s II' 
      IF ( SHORT .EQ. 'SEdm' ) LONG = 'St Edmund Hall' 
      IF ( SHORT .EQ. 'SEH2' ) LONG = 'St Edmund Hall II' 
      IF ( SHORT .EQ. 'SHil' ) LONG = 'St Hilda''s' 
      IF ( SHORT .EQ. 'SHi2' ) LONG = 'St Hilda''s II' 
      IF ( SHORT .EQ. 'SHug' ) LONG = 'St Hugh''s' 
      IF ( SHORT .EQ. 'SJoh' ) LONG = 'St John''s' 
      IF ( SHORT .EQ. 'SJo2' ) LONG = 'St John''s II' 
      IF ( SHORT .EQ. 'SPet' ) LONG = 'St Peter''s' 
      IF ( SHORT .EQ. 'Some' ) LONG = 'Somerville' 
      IF ( SHORT .EQ. 'Trin' ) LONG = 'Trinity' 
      IF ( SHORT .EQ. 'Tri2' ) LONG = 'Trinity II' 
      IF ( SHORT .EQ. 'Univ' ) LONG = 'University'
      IF ( SHORT .EQ. 'Uni2' ) LONG = 'University II'
      IF ( SHORT .EQ. 'Wadh' ) LONG = 'Wadham'
      IF ( SHORT .EQ. 'Wad2' ) LONG = 'Wadham II'
      IF ( SHORT .EQ. 'Wolf' ) LONG = 'Wolfson'
      IF ( SHORT .EQ. 'Wol2' ) LONG = 'Wolfson II'
      IF ( SHORT .EQ. 'Worc' ) LONG = 'Worcester'
      IF ( SHORT .EQ. 'Wor2' ) LONG = 'Worcester II'
      END

      REAL FUNCTION RANGAU ( ISEED )
C
C VERSION
C     29-JUL-91  RJW Add test on zero divide
C     14-DEC-89  RJW  We all make mistakes - for EXTERNAL read INTRINSIC
C     13-DEC-89  RJW  Declare RAN as REAL EXTRENAL
C     13-OCT-89  AD
C
C DESCRIPTION
C
C     Generate random numbers with a Normal Distribution.
C     General Purpose Function.
C
C     The numbers are generated with a Gaussian Distribution, Mean=0, SD=1,
C     i.e. the probability of a value x varies as exp(-x**2/2).
C     The function is a simplified form of GASDEV, in Numerical Recipes 
C     (Press et al.) with the substitution of the VAX RAN function for the 
C     Numerical Recipes RAN1 function, and doesn't save the other random number
C     for the next call. 
C
C     The use of the VAX RAN function changes the requirement of ISEED to be a 
C     large, odd integer (just used internally for RAN function).
C
      IMPLICIT NONE
C
C ARGUMENTS
      INTEGER   ISEED   ! I/O Initially a large, odd integer. Updated each call.
C
C     INTRINSIC RAN     ! VAX VMS Random Number Generator
      REAL      RAND
C
C LOCAL VARIABLES
C
      REAL V1    ! Random number in the range [0,1]
      REAL V2    ! Random number in the range [0,1]
      REAL RAD   ! Random 'Radius'**2 in the range [0,1]
C
C- EXECUTABLE CODE -------------------------------------------------------------
C
 100  CONTINUE
      V1 = 2.0 * RAND (0) - 1.0           ! Convert range [0,1] to [-1,1]
c      iseed = iseed + 3
      V2 = 2.0 * RAND (0) - 1.0  
c      iseed = iseed + 3
      RAD = V1**2 + V2**2
      IF ( RAD .GE. 1.0 .OR. RAD .LE. 0.0 ) GOTO 100 ! Does not lie within unit circle
      RANGAU = V2 * SQRT ( -2.0 * LOG (RAD) / RAD )
C
      END

