C              ANALYSIS OF VARIANCE FOR ONE-WAY DESIGN JUNE 15, 1966
C        THIS IS A SIFTED VERSION OF BMD01V ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
      DOUBLE PRECISION PROBLM,CODE,FINISH,SAMSIZ,SPECTG,PROBNO
      DIMENSION  DATE(2), X(9000), NC(5000), FMT(180), FT(12),
     1 SDEV(12),XBAR(12),NUSE(20),LCODE(72),CONST(72),L1(8),D2(8)
      COMMON  X      , NCI    , NCUE   , FMT    , NTAPE  , LUMP
      COMMON  LCODE  , CONST  , NH     , FM     , NTG    , CASE
C
  200 FORMAT('1BMD01V - ANALYSIS OF VARIANCE FOR ONE-WAY DESIGN - REVISE
     1D JUNE 24, 1969'/
     341H HEALTH SCIENCES COMPUTING FACILITY, UCLA//)
C
      DATA YES,XNO,PROBLM,FINISH,SAMSIZ,SPECTG/3HYES,2HNO,6HPROBLM,
     16HFINISH,6HSAMSIZ,6HSPECTG/
	CALL USAGEB('BMD01V')
      DO  666  I = 1, 20
  666 NUSE(I) = 0
      NEVER = 0
      WRITE (6,200)
 1111 READ (5,300)CODE,PROBNO,NP,OUT,NTG,BIN,WIND,NTAPE,IVF
      IF(CODE.EQ.FINISH)GO TO 1914
      IF(CODE.EQ.PROBLM) GO TO 2
  608 WRITE(6,6080)
 6080 FORMAT(87H0PROBLM OR FINISH IS PUNCHED WRONG ON THE PROPER CARD OR
     1 EITHER CARD IS OUT OF SEQUENCE)
      STOP
    2 IF((NP-1)*(NP-5001)) 6001,6002,6002
 6002 WRITE (6,6102)
 6102 FORMAT(72H0ERROR ON PROBLM CARD, CHECK PROGRAM LIMITS ON NUMBER OF
     1 TREATMENT GROUP)
      STOP
 6005 WRITE(6,6015)
 6015 FORMAT(58H0 WRONG PROBLEM CARD, THE NUMBER OF SPECTG CARD IS TOO B
     1IG)
      STOP
  609 WRITE(6,6019)
 6019 FORMAT(55H00 AND 5 CAN NOT BE SPECIFIED AS INPUT BINARY TAPE UNIT)
      STOP
  610 WRITE(6,6110)
 6110 FORMAT(33H0TOTAL NUMBER OF CASES IS TOO BIG)
      STOP
  406 WRITE(6,4016)J,NI
 4016 FORMAT(11H0THE NUMBER,I3,15H SPECTG CODE ON,I3,23H SPECTG CARD IS 
     1ILLEGAL)
      STOP
  605 WRITE(6,6050)
 6050 FORMAT(31H0ILLEGAL NUMBER OF SPECTG CARDS)
      STOP
  613 WRITE(6,6013)I
 6013 FORMAT(' THE NUMBER OF CASES IN GROUP',I5,' IS OUTSIDE THE RANGE O
     *F BETWEEN 1 AND 9000, PROBLEM WILL BE IGNORED')
      STOP
 6001 IF(NTG-9) 6003,6003,6005
 6003 NCUE=1
      IF(BTWEEN(NP,2,500).LE.0) GO TO 6002
      IF(BIN.NE.YES)GO TO 1113
 1112 NCUE=2
 1113 IF  (BTWEEN(NP, 2, 5000)*BTWEEN(IVF, 0, 10)*BTWEEN(NCUE, 1, 2)*
     1 (BTWEEN(NTAPE, 0, 5) + BTWEEN(NTAPE, 7, 15)))  607, 607, 3
    3 NUTS=NTAPE*(NTAPE-5)
      GO TO  (4, 5), NCUE
    5 IF  (NUTS)  7, 609, 7
    4 IF(IVF.GT.0.AND.IVF.LE.10)GO TO 6
      IVF=1
      WRITE(6,4000)
    6 IF  (NUTS)  7, 8, 7
    8 NTAPE = 5
      GO TO  9
    7 NUSE(NTAPE) = 1
      IF(WIND.EQ.XNO)GO TO 9
  667 REWIND  NTAPE
    9 IF  (NEVER)  10, 10, 11
   11 WRITE (6,200)
   10 NEVER = NEVER + 1
      WRITE (6,400)PROBNO, NP
      NCARDS=(NP+12)/13
      I2=0
      DO 17 J=1,NCARDS
      I1=I2+1
      I2=I2+13
      IF(NP-I2)15,16,16
   15 I2=NP
   16 READ (5,700)CODE,(NC(I),I=I1,I2)
      IF(CODE.EQ.SAMSIZ)GO TO 17
      WRITE(6,7000)
7000  FORMAT(55H0SAMSIZ IS SPELLED WRONG OR SAMSIZ CARD OUT OF SEQUENCE)
      CALL EXIT
  165 NCARDS=-NCARDS
   17 CONTINUE
      IF(-NCARDS)18,18,6002
   18 CASE=0.0
      DO 19 I=1,NP
      IF(BTWEEN(NC(I),1,9000))613,613,181
  181 CASE1=NC(I)
      CASE=CASE+CASE1
   19 CONTINUE
      IF(CASE-(10.0**8))195,195,610
  195 IF(NTG)605,499,401
  401 NH=0
      DO 445 NI=1,NTG
      READ(5,425)CODE,M,(L1(I),D2(I),I=1,M)
      IF(CODE.EQ.SPECTG)GO TO 404
  403 WRITE (6,405)
      STOP
  404 DO 6004  J=1,M
      IF(L1(J)*(L1(J)-11)) 6004,406,406
 6004 CONTINUE
      DO 415 NG=1,M
      NG1=NH+NG
      LCODE(NG1)=L1(NG)
  415 CONST(NG1)=D2(NG)
  445 NH=NH+M
  499 GO TO (12,13),NCUE
   12 WRITE (6,500)IVF, NTAPE
      WRITE(6,750)
  750 FORMAT(20H0THE VARIABLE FORMAT)
      IVF = IVF*18
      READ (5,100)(FMT(I), I = 1, IVF)
      WRITE(6,101)(FMT(I),I=1,IVF)
 101  FORMAT(1X,18A4)
      GO TO  14
   13 READ (5,601)LUMP
      WRITE (6,600)LUMP, NTAPE
   14 IF(NTG)605,439,426
  426 WRITE (6,427)
      DO 428 NG=1,NH
  428 WRITE (6,431)LCODE(NG),CONST(NG)
  439 SPSQ=0.0
      SMSQ = 0.0
      GT = 0.0
      NP1 = NP - 1
      NCI = NC(1)
      TOBS = NCI
      CALL  READIN
      DO  20  I = 1, NCI
   20 GT = GT + X(I)
      Y = GT/TOBS
      DO  21  J = 1, NCI
   21 SPSQ = SPSQ + (X(J) - Y)**2
      NGO = 1
      IF(OUT.EQ.YES)GO TO 31
   30 WRITE (6,801)
      GO TO  40
   31 SDV = SQRT( SPSQ/(TOBS - 1.0) )
      CALL  DECFND(Y, LY)
      CALL  DECFND(SDV, LSDV)
      IF  (NP - 12)  32, 32, 33
   33 WRITE (6,900)
      LY= MAX0(0, 7 - LY)
      LSDV= MAX0(0, 7 - LSDV)
      WRITE (6,901)NGO, NCI, Y, SDV
      NGO = 3
      GO TO  40
   32 LY= MAX0(LY, LSDV)
      LSDV= MAX0(0, 5 - LY)
      XBAR(1) = Y
      SDEV(1) = SDV
      NGO = 2
   40 DO  50  I = 2, NP
      NCI = NC(I)
      XNI = NCI
      TOBS = TOBS + XNI
      CALL  READIN
      AVE = 0.0
      DO  51  J = 1, NCI
   51 AVE = AVE + X(J)
      GT = GT + AVE
      AVE = AVE/XNI
      SDV = 0.0
      DO  52  J = 1, NCI
   52 SDV = (X(J) - AVE)**2 + SDV
      SPSQ = SPSQ + SDV
      SMSQ = SMSQ + XNI*(AVE - Y)**2
      SDV = SQRT( SDV/(XNI - 1.0) )
      GO TO  (50, 61, 62), NGO
   61 XBAR(I) = AVE
      SDEV(I) = SDV
      GO TO  50
   62 WRITE (6,901)I, NCI, AVE, SDV
   50 CONTINUE
      GO TO  (70, 71, 70), NGO
   71 WRITE (6,1000)(I, I = 1, NP)
 1000 FORMAT(////22H TREATMENT GROUP       , 12(I7, 2X) )
      WRITE (6,1001)(NC(I), I = 1, NP)
      WRITE (6,1002)(XBAR(I), I = 1, NP)
      WRITE (6,1003)(SDEV(I), I = 1, NP)
   70 WRITE (6,1100)
      GT = GT/TOBS
      SMSQ = SMSQ - TOBS*(GT - Y)**2
      GT = SMSQ + SPSQ
      TOBS = TOBS - 1.0
      TLEVEL = NP1
      DF = TOBS - TLEVEL
      AVE = SMSQ/TLEVEL
      SDV = SPSQ/DF
      FRATIO = AVE/SDV
      WRITE (6,1200)SMSQ, NP1, AVE, FRATIO
      INTDF=DF +.0000001
      INTOBS=TOBS+.0000001
      WRITE (6,1500)SPSQ, INTDF   , SDV
      WRITE (6,1600)GT, INTOBS
      GO TO  1111
  607 WRITE (6,6070)
      STOP
 1914 DO  73  I = 1, 20
      IF  (NUSE(I))  72, 73, 72
   72 IF(I.NE.5)REWIND I
   73 CONTINUE
      WRITE (6,1900)
 1915 STOP
 100  FORMAT(18A4)
  300 FORMAT(2A6,I4,A3,I1,43X,A3,A2,2I2)
  400 FORMAT(///// 15H PROBLEM CODE   , A6, /  20H NUMBER OF TREATMENT ,
     1 9H GROUPS  , I4  )
  405 FORMAT(51H0ERROR ON TRANS-GENERATION CARD, SPECTG IS MISSPELT)
  425 FORMAT(A6,I1,8(I2,F6.0))
  427 FORMAT(1H023HSPECIAL TRANSGENERATION/1H 3X4HCODE4X8HCONSTANT)
  431 FORMAT(1H 5XI2,1XF11.5)
  500 FORMAT(34H NUMBER OF VARIABLE FORMAT CARDS   , I2  /
     1 18H DATA INPUT TAPE   , I2 )
  600 FORMAT(23H BINARY RECORD LENGTH  , I3 / 18H DATA INPUT TAPE  , I2)
  601 FORMAT(I3)
  700 FORMAT(A6, 1X, 13I5)
  801 FORMAT(////41H LISTING OF TREATMENT MEANS NOT REQUESTED   )
  900 FORMAT(//// 49H TREATMENT    SAMPLE SIZE      MEAN      STANDARD
     1 10H DEVIATION   //)
  901 FORMAT    (1X,I6,8X,I6,F16.5,        F17.5)
 1001 FORMAT(/ 22H SAMPLE SIZE           , 12(I7, 2X)   )
 1002 FORMAT    (/22H MEAN                 , 12F9.4)
 1003 FORMAT    (/22H STANDARD DEVIATION   , 12F9.4)
 1100 FORMAT(/////  33X, 21H ANALYSIS OF VARIANCE     ///
     1 19X, 15H SUM OF SQUARES, 6X, 2HDF, 5X, 12H MEAN SQUARE , 5X,
     2 8H F RATIO     /  )
 1200 FORMAT(15H0BETWEEN GROUPS      , F17.4, I10, 2F15.4  )
 1300 FORMAT(2F11.0)
 1400 FORMAT(2A5, 1X, 2A5)
 1500 FORMAT(15H0WITHIN GROUPS    , F17.4, I10, F15.4)
 1600 FORMAT(// 1H , 5X, 5HTOTAL, F21.4, I10)
 1900 FORMAT(40H1COMPUTATION HALTED BY  FINISH  CARD.   )
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
 6070 FORMAT(123H0PROBLEM CARD IS SET UP WRONG, EITHER INPUT TAPE UNIT I
     1S SPECIFIED WRONG OR THE NUMBER OF F-TYPE VARIABLE FORMAT IS TOO B
     2IG)
      END
C        FUNCTION BTWEEN FOR BMD01V                   JUNE 15, 1966
C        BTWEEN = 1.0 IF N IS BETWEN I AND J
      FUNCTION  BTWEEN(N, I, J)
      BTWEEN=0.0
      IF(N.GE.I.AND.N.LE.J)BTWEEN=1.0
      RETURN
      END
C        SUBROUTINE DECFND FOR BMD01V                 JUNE 15, 1966
C        X HAS A MAGNITUDE BETWEEN 10**(LEFT-1) AND 10**(LEFT)
      SUBROUTINE  DECFND(X, LEFT)
      A=10.
      DO 1 LEFT=1,19
      IF(A.GT.X) GO TO 2
    1 A=A*10.
    2 RETURN
      END
C        SUBROUTINE DIVALG FOR BMD01V                 JUNE 15, 1966
C        DIVISION ALGORITHM
      SUBROUTINE  DIVALG(NA, NB, NQ, NR, NCODE)
      NN = NA
      IF  (NCODE)  1, 1, 2
    2 NN = NN - 1
    1 NQ = NN/NB
      NR = NA - NB*NQ
      RETURN
      END
C        SUBROUTINE READIN FOR BMD01V                 JUNE 15, 1966
      SUBROUTINE  READIN
      DIMENSION X(9000),LCODE(72),CONST(72),FMT(180)
      COMMON  X      , NCI    , NCUE   , FMT    , NTAPE  , LUMP
      COMMON  LCODE  , CONST  , NH     , FM     , NTG
      GO TO  (1, 2), NCUE
    1 READ (NTAPE,FMT)(X(I), I = 1, NCI)
    6 IF(NTG)7,8,7
    7 CALL TRNGEN
    8 RETURN
    2 CALL  DIVALG(NCI, LUMP, NDO, LEFT, 1)
      N2 = 0
      IF  (NDO)  3, 3, 4
    4 DO  5  J = 1, NDO
      N1 = N2 + 1
      N2 = N2 + LUMP
    5 READ(NTAPE,9876)(X(I),I=N1,N2)
    3 N1 = N2 + 1
      READ(NTAPE,9876)(X(I),I=N1,NCI)
      GO TO 6
 9876 FORMAT(20A4)
      END
C        SUBROUTINE TRNGEN FOR BMD01V                 JUNE 15, 1966
      SUBROUTINE TRNGEN
      DIMENSION X(9000),LCODE(72),CONST(72),FMT(180)
      COMMON  X      , NCI    , NCUE   , FMT    , NTAPE  , LUMP
      COMMON  LCODE  , CONST  , NH     , FM     , NTG    , SAMP
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
      DO 210 I=1,NH
      FM=CONST(I)
      JUMP=LCODE(I)
      DO 85 J=1,NCI
      D=X(J)
      GO TO(10,20,30,40,50,60,70,80,90,110),JUMP
   10 IF(D)99,11,12
   11 D=0.0
      GO TO 8
   12 D=SQRT(D)
      GO TO 8
   20 IF(D)99,21,22
   21 D=1.0
      GO TO 8
   22 D=SQRT(D)+SQRT(D+1.0)
      GO TO 8
   30 IF(D)99,99,31
   31 D=ALOG10(D)
      GO TO 8
   40 D=EXP(D)
      GO TO 8
   50 IF(-D)52,11,99
   52 IF(D-1.0)53,54,99
   54 D=3.14159265/2.0
      GO TO 8
   53 D=ASN(SQRT(D))
      GO TO 8
   60 A=D/(SAMP+1.0)
      B=A+1.0/(SAMP+1.0)
      IF(A)99,61,62
   61 IF(-B)64,11,99
   64 D=ASN(SQRT(B))
      GO TO 8
   62 IF(B)99,65,66
   65 D=ASN(SQRT(A))
      GO TO 8
   66 D=ASN(SQRT(A))+ASN(SQRT(B))
      GO TO 8
   70 IF(D)71,99,71
   71 D=1.0/D
      GO TO 8
   80 D=D+FM
      GO TO 8
   90 D=D*FM
      GO TO 8
  110 IF(D)111,11,111
  111 D=D**FM
    8 X(J)=D
   85 CONTINUE
  210 CONTINUE
      GO TO 1000
   99 WRITE(6,105)I,D
      STOP
  105 FORMAT(23H0TRANS-GENERATION ERROR//10H PASS NO.=I3,9H X VALUE=F10.
     15)
 1000 RETURN
      END










































