C        PERIODIC REGRESSION AND HARMONIC ANALYSIS   MARCH 29, 1966
C        THIS IS A SIFTED VERSION OF BMD04R ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DIMENSION III(8),CCC(8)
      DOUBLE PRECISION A123,B123,C123,E123
       DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
      DATA A123,B123,C123,D123,E123/6HPROBLM,6HFINISH,6HCOVPRB,3HYES,
     X6HSPECTG/
C
  425 FORMAT(83H1BMD04R - PERIODIC REGRESSION AND HARMONIC ANALYSIS - VE
     1RSION OF MARCH 29, 1966    ,/
     240H HEALTH SCIENCES COMPUTING FACILITY,UCLA//
     3 14H PROBLEM CODE A6,/
     4 14H NO. OF REPS. I3,/
     518H NO. OF SUB-INTS. I4,/
     6 18H NO. OF HARMONICS I2,/)
C
 404  FORMAT(' A PROBLEM OR FINNISH CARD WAS EXPECTED.')
  412 FORMAT(A6,I1,8(I2,F6.0))
 413  FORMAT(18A4)
  501 FORMAT(2A6,I2,I3,A3,2I2,I1,A3,40X,2I2)
  509 FORMAT(1H I4,F11.4,F12.4,F11.4)
  510 FORMAT(1H010X,18HTABLE OF RESIDUALS/9X,19HFOR HARMONIC NUMBERI2//)
  511 FORMAT(15H  TIME   ACTUAL4X,20HPREDICTED   RESIDUAL//)
  512 FORMAT(2A6,A3,56X,I2)
 1017 FORMAT(' THE CARD THAT WAS READ-IN CONTAINS AT LEAST ONE ERROR.')
C
C
 4000 FORMAT(64H0ILLEGAL TRANSGENERATION CODE SPECIFIED. PROGRAM CANNOT   
     1PROCEED.)
 4001 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
 4201 FORMAT(' A COVARIANCE PROBLEM CARD WAS EXPECTED.')
 4202 FORMAT(' A SPECIAL TRANSGENERATION CARD6WAS EXPECTED.')
      NTAPE=5
	CALL USAGEB('BMD04R')
C
C     READ THE PROBLEM CARD
C
  888 READ (5,501)TODE,CODE,NFR,KNT,TRANS,KOV,ITAB,NHARM,XPLOT,MTAPE,MAT
     1Y
      IE=0
      IF(TODE.EQ.B123)GO TO 401
      IF(TODE.EQ.A123)GO TO 403
      WRITE (6,404)
 402  WRITE(6,1017)
  401 IF(NTAPE.EQ.5)GO TO 456
      REWIND NTAPE
C
C
C
  456 STOP
  403 CALL TPWD(MTAPE,NTAPE)
  459 IF(NFR*(NFR-20)) 604,402,402
 604  IF((KNT-NHARM-NHARM)*(KNT-401))405,402,402
  405 IF(NHARM*(NHARM-10)) 406,402,402
  406 IF(MATY.GT.0.AND.MATY.LE.10)GO TO 407
      WRITE(6,4001)
      MATY=1
C407  MATY=MATY*12
 407  MATY=MATY*18
      NDF1=NFR-1
      TOTOB=KNT*NFR
      WRITE (6,425)CODE,NFR,KNT,NHARM
      IF(TRANS.NE.D123)GO TO 4105
C
C     READ THE SPECIAL TRANSGENERATION CARDS
C
  411 READ (5,412) SPC,NNN,(III(I),CCC(I),I=1,NNN)
      IE=0
      IF(SPC.NE.E123)GO TO 402
 409  DO 410 I=1,NNN
      IF(III(I).EQ.0.OR.III(I).GE.11)GO TO 310
 410  CONTINUE
 4105 IF(IE)401,325,325
 310  WRITE (6,4000)
      IE=-11
      GO TO 410
C
C     READ THE VARIBLE FORMAT CARD
C
 325  READ (5,413)(TSUMX(I),I=1,MATY)
      WRITE(6,3251) (TSUMX(I), I=1,MATY)
 3251 FORMAT(' THE VARIABLE FORMAT IS  ',18A4/(1H ,24X,18A4))
      DO 414 I=1,NFR
C
C     READ THE INPUT DATA
C
  414 READ (NTAPE,TSUMX)(YDATA(I,J),J=1,KNT)
      IF(TRANS.NE.D123)GO TO 3
  416 DO 417 I=1,NFR
      DO 418 J=1,KNT
  418 EQY(J)=YDATA(I,J)
      CALL TRANS1(EQY,KNT,NNN,III,CCC,I)
      DO 417 J=1,KNT
  417 YDATA(I,J)=EQY(J)
    3 SUMT2=0.0
      SUMT=0.0
      FN=NFR
      TIME=KNT
      C=(2.0*3.14159265)/TIME
      DO101I=1,9
      DU(I)=0.0
      DV(I)=0.0
      AAY(I)=0.0
      BBY(I)=0.0
      DO101J=1,NFR
      AY(J,I)=0.0
      BY(J,I)=0.0
      RSUM(J)=0.0
  101 CONTINUE
      DO103J=1,KNT
  103 TSUM(J)=0.0
      DO21I=1,KNT
      DO20J=1,NFR
      RSUM(J)=RSUM(J)+YDATA(J,I)
      TSUM(I)=TSUM(I)+YDATA(J,I)
   20 SUMT2=SUMT2+YDATA(J,I)*YDATA(J,I)
   21 SUMT=SUMT+TSUM(I)
      ZEANY=SUMT/(FN*TIME)
      ICURVE=0
      IF(NHARM)66,66,67
   66 NHARM=1
      GOTO69
   67 ICURVE=NHARM
      NHARM=1
   69 N=NHARM
      H=N
      DO10I=1,KNT
      T=I-1
      U(I,N)=COS(C*T*H)
   10 V(I,N)=SIN(C*T*H)
      DO15I=1,KNT
      DU(N)=DU(N)+U(I,N)*U(I,N)
   15 DV(N)=DV(N)+V(I,N)*V(I,N)
      DO30I=1,NFR
      DO25J=1,KNT
      AY(I,N)=AY(I,N)+U(J,N)*YDATA(I,J)
   25 BY(I,N)=BY(I,N)+V(J,N)*YDATA(I,J)
      AY(I,N)=AY(I,N)/DU(N)
   30 BY(I,N)=BY(I,N)/DV(N)
      DO40J=1,KNT
      AAY(N)=AAY(N)+U(J,N)*TSUM(J)
   40 BBY(N)=BBY(N)+V(J,N)*TSUM(J)
      AAY(N)=AAY(N)/(FN*DU(N))
      BBY(N)=BBY(N)/(FN*DV(N))
      DO50J=1,KNT
      ADD=0.0
      DO45I=1,NHARM
   45 ADD=ADD+AAY(I)*U(J,I)+BBY(I)*V(J,I)
   50 EQY(J)=ZEANY+ADD
      IF(ITAB.NE.99)GO TO 80
   81 WRITE (6,510)NHARM
      WRITE (6,511)
      II=0
   82 II=II+1
      YBAR=TSUM(II)/FN
      DIFF=YBAR-EQY(II)
      IJK=II-1
      WRITE (6,509)IJK,YBAR,EQY(II),DIFF
      IF(II.LT.KNT)GO TO 82
   80 IF(NDF1.LE.0)GO TO 2000
 2001 CALL ANAVAR
      S6=SUMT2-CM
      S7=CM
      S1M=S1/DF1
      S2M(N)=S2(N)/DF2(N)
      S3M=S3/DF3
      S5M=S5/DF5
      FV1=S1M/S5M
      FV3=S3M/S5M
      S4M(N)=S4(N)/DF4(N)
      S8M(N)=S8(N)/DF8(N)
      S9M(N)=S9(N)/DF9(N)
      DO690 I=1,N
      FV4(I)=S4M(I)/S5M
      FV2(I)=S2M(I)/(S3M+S4M(I)-S5M)
      ESTD(I)=(S3M+S4M(I)-S5M)**2
      BOT=S3M**2/DF3+S4M(I)**2/DF4(I)+S5M**2/DF5
      ESTD(I)=ESTD(I)/BOT
      FV8(I)=S8M(I)/S5M
  690 FV9(I)=S9M(I)/S5M
 2000 IF(N-2) 2006,2007,2007
 2006 CALL MAXMIN
      AMP=2.0*SQRT(AAY(N)**2+BBY(N)**2)
      SEAMP=AMP/2.0
 2007 CALL PRITE
      IF(NDF1.LE.0)GO TO 2002
 2003 CALL VRITE
 2002 IF(XPLOT.NE.D123)GO TO 3001
 3000 CALL PLOT
 3001 IF(NHARM-9)90,100,100
   90 IF(KNT-2*N-1)100,100,91
   91 IF(ICURVE-NHARM)2004,100,2004
 2004 NHARM=NHARM+1
      GOTO69
 4021 WRITE(6,4201)
      GO TO 402
 4022 WRITE(6,4202)
      GO TO 402
  100 IF(NDF1.LE.0.OR.KOV.LE.0)GO TO 888
 1101 DO 102 KK=1,KOV
C
C     READ COVARIANCE PROBLEM CARD
C
      READ (5,512) ODE,CODEX,TRANS,MATX
      IE=0
      IF(ODE.NE.C123)GO TO 4021
      IF(TRANS.NE.D123)GO TO 419
 420  READ (5,412) SPC,NNN,(III(I),CCC(I),I=1,NNN)
      IF(SPC.NE.E123)GO TO 4022
 322  DO 350 I=1,NNN
      IF(III(I).GT.0.AND.III(I).LT.11)GO TO 350
 330  WRITE (6,4000)
      IE=-11
 350  CONTINUE
      IF(-IE)419,419,401
  419 IF(MATX.GT.0.AND.MATX.LE.10)GO TO 4195
      WRITE(6,4001)
      MATX=1
 4195 MATX=MATX*18
C
      READ (5,413)(EQY(I),I=1,MATX)
      WRITE(6,3252) (EQY(I), I=1,MATX)
 3252 FORMAT(' THE DATA INPUT FORMAT IS  ', 18A4/(1H ,26X,18A4))
      CALL COVAR(NTAPE,TRANS,D123,NNN,CODEX)
      CALLCRITE
  102 CONTINUE
      GO TO 888
      END
      SUBROUTINEANAVAR
C              SUBROUTINE ANAVAR FOR BMD04R          MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
C
      IF(N-2)90,100,100
   90 CM=SUMT**2/(FN*TIME)
      SS1=0.0
      DO3I=1,NFR
    3 SS1=SS1+RSUM(I)*RSUM(I)
      S1=SS1/TIME-CM
      NDF1=NFR-1
      DF1=NDF1
      NDF6=NFR*KNT-1
      DF6=NDF6
      NDF7=1
      DF7=1.0
      SS3=0.0
      DO4J=1,KNT
    4 SS3=SS3+TSUM(J)*TSUM(J)
  100 ADDS2=0.0
      S2(N)=(DU(N)*AAY(N)**2+DV(N)*BBY(N)**2)*FN
      NDF2(N)=2
      DF2(N)=2.0
      NDF4(N)=2*(NFR-1)
      DF4(N)=2.0*(FN-1.0)
      NDF8(N)=NFR-1
      DF8(N)=FN-1.0
      NDF9(N)=NFR-1
      DF9(N)=FN-1.0
      NDF3=KNT-2*N-1
      DF3=NDF3
      NDF5=(NFR-1)*NDF3
      DF5=NDF5
      DO5I=1,NHARM
    5 ADDS2=ADDS2+S2(I)
      S3=SS3/FN-CM-ADDS2
      S4(N)=0.0
      ADDS8(N)=0.0
      DO7I=1,NFR
      S4(N)=S4(N)+DU(N)*AY(I,N)**2+DV(N)*BY(I,N)**2
    7 ADDS8(N)=ADDS8(N)+SQRT(DU(N)*AY(I,N)**2+DV(N)*BY(I,N)**2)
      S4(N)=S4(N)-S2(N)
      S8(N)=2.0*S2(N)+S4(N)-2.0*SQRT(S2(N))*(ADDS8(N)/SQRT(FN))
      S9(N)=S4(N)-S8(N)
      ADDS4=0.0
      DO8I=1,NHARM
    8 ADDS4=ADDS4+S4(I)
      S5=SUMT2+CM-SS1/TIME-SS3/FN-ADDS4
      RETURN
      END
      SUBROUTINE COVAR(NTAPE,TRANS,D123,NNN,CODEX)
C           SUBROUTINE COVAR FOR BMD04R              MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
C
  513 FORMAT(18H1COVARIANCE CODE  A6)
      DO101I=1,NHARM
      AAX(I)=0.0
      BBX(I)=0.0
      DO101J=1,NFR
      AX(J,I)=0.0
      BX(J,I)=0.0
      RSUMX(J)=0.0
  101 CONTINUE
      DO 102 J=1,KNT
      TSUMX(J)=0.0
  102 CONTINUE
      SUMT2X=0.0
      SUMTX=0.0
      SUMXY=0.0
      DO 21 J=1,NFR
      READ(NTAPE,EQY)(XDATA(I),I=1,KNT)
      IF(TRANS.NE.D123)GO TO 103
      CALL TRANS1(XDATA,KNT,NNN,III,CCC,J)
  103 DO 20 I=1,KNT
      RSUMX(J)=RSUMX(J)+XDATA(I)
      TSUMX(I)=TSUMX(I)+XDATA(I)
      SUMT2X=SUMT2X+XDATA(I)*XDATA(I)
      SUMXY=SUMXY+YDATA(J,I)*XDATA(I)
      DO 20 N=1,NHARM
      AX(J,N)=AX(J,N)+U(I,N)*XDATA(I)
   20 BX(J,N)=BX(J,N)+V(I,N)*XDATA(I)
   21 CONTINUE
      WRITE (6,513)CODEX
      DO25J=1,KNT
   25 SUMTX=SUMTX+TSUMX(J)
      ZEANX=SUMTX/(FN*TIME)
      DO29N=1,NHARM
      DO30I=1,NFR
      AX(I,N)=AX(I,N)/DU(N)
   30 BX(I,N)=BX(I,N)/DV(N)
      DO40J=1,KNT
      AAX(N)=AAX(N)+U(J,N)*TSUMX(J)
   40 BBX(N)=BBX(N)+V(J,N)*TSUMX(J)
      AAX(N)=AAX(N)/(FN*DU(N))
      BBX(N)=BBX(N)/(FN*DV(N))
   29 CONTINUE
      DO50J=1,KNT
      ADD=0.0
      DO45I=1,NHARM
   45 ADD=ADD+AAX(I)*U(J,I)+BBX(I)*V(J,I)
   50 EQX(J)=ZEANX+ADD
      CALLCROSS
      CMX=SUMTX**2/(FN*TIME)
      SS1X=0.0
      DO3I=1,NFR
    3 SS1X=SS1X+RSUMX(I)*RSUMX(I)
      S1X=SS1X/TIME-CMX
      SS3X=0.0
      DO4J=1,KNT
    4 SS3X=SS3X+TSUMX(J)*TSUMX(J)
      ADDS2X=0.0
      ADDS4X=0.0
      DO22N=1,NHARM
      S4X(N)=0.0
      S2X(N)=(DU(N)*AAX(N)**2+DV(N)*BBX(N)**2)*FN
      DO7I=1,NFR
    7 S4X(N)=S4X(N)+DU(N)*AX(I,N)**2+DV(N)*BX(I,N)**2
      S4X(N)=S4X(N)-S2X(N)
   22 CONTINUE
      DO5I=1,NHARM
    5 ADDS2X=ADDS2X+S2X(I)
      S3X=SS3X/FN-CMX-ADDS2X
      DO8I=1,NHARM
    8 ADDS4X=ADDS4X+S4X(I)
      S5X=SUMT2X+CMX-SS1X/TIME-SS3X/FN-ADDS4X
      EXX1=S5X+S1X
      EXX3=S5X+S3X
      EXY1=S5XY+S1XY
      DO 26 N=1,NHARM
      EXX2(N)=S5X+S2X(N)
      EXX4(N)=S5X+S4X(N)
      EXY4(N)=S5XY+S4XY(N)
      EXY2(N)=S5XY+S2XY(N)
   26 CONTINUE
      EXY3=S5XY+S3XY
      NDFE1=NDF1+NDF5
      DO9I=1,NHARM
      NDFE2(I)=NDF2(I)+NDF5
    9 NDFE4(I)=NDF4(I)+NDF5
      NDFE3=NDF3+NDF5
      BE2=S5XY**2/S5X
      BS12=EXY1**2/EXX1-BE2
      BS32=EXY3**2/EXX3-BE2
      NDFERR=NDF5-1
      DFERR=NDFERR
      BE2M=(S5-BE2)/DFERR
      BS12M=(S1-BS12)/DF1
      BS32M=(S3-BS32)/DF3
      FC1=BS12M/BE2M
      FC3=BS32M/BE2M
      BYX=S5XY/S5X
      FC5=BE2/BE2M
      DO 31 N=1,NHARM
      BS22(N)=EXY2(N)**2/EXX2(N)-BE2
      BS42(N)=EXY4(N)**2/EXX4(N)-BE2
      BS42M(N)=(S4(N)-BS42(N))/DF4(N)
      BS22M(N)=(S2(N)-BS22(N))/DF2(N)
      FC2(N)=BS22M(N)/(BS32M+BS42M(N)-BE2M)
      FC4(N)=BS42M(N)/BE2M
      XNUM=(BS32M+BS42M(N)-BE2M)**2
      DEN=BS32M**2/DF3+BS42M(N)**2/DF4(N)+BE2M**2/DFERR
      DFESTC(N)=XNUM/DEN
   31 CONTINUE
      RETURN
      END
      SUBROUTINECRITE
C             SUBROUTINE CRITE FOR BMD04R            MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
C
    4 FORMAT(1H 45X22HANALYSIS OF COVARIANCE//)
    5 FORMAT(20H SOURCE OF VARIATION5X2HDF8X2HXX14X2HXY14X20HYY     *       
     1  BS2-BE24X2HDF7X6HRED.MS5X8HF RATIOS)
    6 FORMAT(22H BETWEEN REPLICATES   I5,F14.4,2F16.4,2H *F13.4,I5,F14.4
     1,F12.3)
    7 FORMAT(13H EFFECT OF (AI1,2H+BI1,5H)    I5,F14.4,2F16.4,2H *F13.4,
     1I5,F14.4,F12.3,1H*)
    8 FORMAT(22H SCATTER ABOUT CURVE  I5,F14.4,2F16.4,2H *F13.4,I5,F14.4
     1,F12.4)
    9 FORMAT(15H REPLICATE X (AI1,2H+BI1,3H)  I5,F14.4,2F16.4,2H *F13.4,
     1I5,F14.4,F12.4)
   10 FORMAT(6H ERROR16XI5,F14.4,2F16.4,2H *F13.4,I5,F14.4,F12.4)
   11 FORMAT(1H 74(1H*))
   12 FORMAT(22H ERROR+REPLICATES     I5,F14.4,F16.4)
   13 FORMAT(9H ERROR+(AI1,2H+BI1,9H)        I5,F14.4,F16.4,18X24H*EST D
     1F FOR DENOMINATOR I1,7H ABOVE=F6.2)
   14 FORMAT(14H ERROR+SCATTER8XI5,F14.4,F16.4,19X29HREGRESSION COEFFICI
     1ENT B(YX)=F8.4)
   15 FORMAT(16H ERROR+REPS X (AI1,2H+BI1,2H) I5,F14.4,F16.4)
   16 FORMAT(1H075X34HBE2=(XY).(XY)/XX    FROM ERROR ROW)
   17 FORMAT(1H 75X42HBS2=(XY).(XY)/XX    FROM SUMS IN LAST ROWS)
  100 FORMAT(1H 119(1H*))
  101 FORMAT(1H 75X,24HRED.MS=(YY-(BS2-BE2))/DF)
      NDFE1=NDF1+NDF5
      WRITE (6,4)
      WRITE (6,5)
      WRITE (6,100)
      WRITE (6,6)NDF1,S1X,S1XY,S1,BS12,NDF1,BS12M,FC1
      DO20I=1,NHARM
   20 WRITE (6,7)I,I,NDF2(I),S2X(I),S2XY(I),S2(I),BS22(I), NDF2(I),BS22M
     1(I),FC2(I)
      WRITE (6,8)NDF3,S3X,S3XY,S3,BS32,NDF3,BS32M,FC3
      DO25I=1,NHARM
   25 WRITE (6,9)I,I,NDF4(I),S4X(I),S4XY(I),S4(I),BS42(I),NDF4(I),BS42M(
     1I),FC4(I)
      WRITE (6,10)NDF5,S5X,S5XY,S5,BE2,NDFERR,BE2M,FC5
      WRITE (6,11)
      WRITE (6,12)NDFE1,EXX1,EXY1
      DO30I=1,NHARM
   30 WRITE (6,13)I,I,NDFE2(I),EXX2(I),EXY2(I),I,DFESTC(I)
      WRITE (6,14)NDFE3,EXX3,EXY3,BYX
      DO35I=1,NHARM
   35 WRITE (6,15)I,I,NDFE4(I),EXX4(I),EXY4(I)
      WRITE (6,100)
      WRITE (6,16)
      WRITE (6,17)
      WRITE (6,101)
      RETURN
      END
      SUBROUTINECROSS
C           SUBROUTINE CROSS FOR BMD04R              MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
C
      SS3XY=0.0
      SS1XY=0.0
      ADDXY=0.0
      CMXY=SUMT*SUMTX/(FN*TIME)
      DO10I=1,NFR
   10 SS1XY=SS1XY+RSUMX(I)*RSUM (I)
      S1XY=SS1XY/TIME-CMXY
      DO 15 J=1,KNT
   15 SS3XY=SS3XY+TSUMX(J)*TSUM (J)
      ADD4=0.0
      DO40N=1,NHARM
      S4XY(N)=0.0
      S2XY(N)=(AAX(N)*AAY(N)*DU(N)+BBX(N)*BBY(N)*DV(N))*FN
   40 CONTINUE
      DO20I=1,NHARM
   20 ADDXY=ADDXY+S2XY(I)
      S3XY=SS3XY/FN-CMXY-ADDXY
      DO45N=1,NHARM
      DO25I=1,NFR
   25 S4XY(N)=S4XY(N)+DU(N)*AX(I,N)*AY(I,N)+DV(N)*BX(I,N)*BY(I,N)
      S4XY(N)=S4XY(N)-S2XY(N)
   45 CONTINUE
      S6XY=SUMXY-CMXY
      S7XY=CMXY
      DO30I=1,NHARM
   30 ADD4=ADD4+S4XY(I)
      S5XY=SUMXY+CMXY-SS1XY/TIME-SS3XY/FN-ADD4
      RETURN
      END
      SUBROUTINEMAXMIN
C            SUBROUTINE MAXMIN FOR BMD04R            MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
C
      RANGE=2.0*SQRT(AAY(N)**2+BBY(N)**2)
      RATIO=ABS(BBY(N)/AAY(N))
      TP=ATAN(RATIO)
      IF(AAY(N))5,20,20
    5 IF(BBY(N))10,15,15
   10 THETA=3.14159265+TP
      GOTO35
   15 THETA=3.14159265-TP
      GOTO35
   20 IF(BBY(N))25,30,30
   25 THETA=6.2831853-TP
      GOTO35
   30 THETA=TP
   35 RMAX=TIME*THETA/6.2831853
C
      RMIN=RMAX+0.5*TIME
      RETURN
      END
      SUBROUTINE PLOT
C         SUBROUTINE PLOT FOR BMD04R                 MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DIMENSION ABCD(101)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
      DATA Q000HL,Q001HL,Q002HL,Q003HL/4H    ,4H.   ,4HO   ,4HP   /
C
      XMAX=-10.0**10
      XMIN=10.0**10
      DO 3 I=1,KNT
      EQX(I)=TSUM(I)/FN
      XMAX=AMAX1(XMAX,EQX(I),EQY(I))
      XMIN=AMIN1(XMIN,EQX(I),EQY(I))
    3 CONTINUE
C     DIMENSION ABCD(101)
      DO 1 I=1,6
 1    S2X(I)=(XMAX-XMIN)*(FLOAT(I)-1.0)*.2+XMIN
      WRITE (6,2)(S2X(I),I=1,6)
 2    FORMAT(6H- TIMEF14.4,5F20.4/1H07X,11(9X,1H.))
      DO 4 I=4,100
 4    ABCD(I)=(+Q000HL)
      K=2
      L=3
      DDD=100.0/(XMAX-XMIN)
      DO 5 I=1,KNT
      ABCD(K)=(+Q000HL)
      ABCD(L)=(+Q000HL)
      ABCD(1)=(+Q001HL)
      ABCD(101)=(+Q001HL)
      K=(EQX(I)-XMIN)*DDD+1.5
      L=(EQY(I)-XMIN)*DDD+1.5
      ABCD(K)=(+Q002HL)
      ABCD(L)=(+Q003HL)
      I1=I-1
 5    WRITE (6,6)I1,ABCD
 6    FORMAT(I5,12X,101A1)
      RETURN
      END
      SUBROUTINE PRITE
C              SUBROUTINE PRITE FOR BMD04R           MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
C
    3 FORMAT(20H0PARAMETER ESTIMATES)
    4 FORMAT(16H0CONSTANT(MEAN)=F13.4)
    5 FORMAT(13H0COEFFICIENTS)
    6 FORMAT(3H A(I1,2H)=F12.4,4H  B(I1,2H)=F12.4)
    7 FORMAT(17H0SEMI-AMPLITUDE =10X,F12.4)
    8 FORMAT(23H PHASE ANGLE (RADIANS)=4X,F12.4)
    9 FORMAT(8H RANGE =19X,F12.4)
   10 FORMAT(27H TIME OF MAXIMUM RESPONSE =F12.4)
   11 FORMAT(27H TIME OF MINIMUM RESPONSE =F12.4)
      WRITE (6,3)
      WRITE (6,4)ZEANY
      WRITE (6,5)
      DO20I=1,NHARM
   20 WRITE (6,6)I,AAY(I),I,BBY(I)
      IF(N-2)21,22,22
   21 WRITE (6,7)SEAMP
      WRITE (6,8)THETA
      WRITE (6,9)RANGE
      WRITE (6,10)RMAX
      WRITE (6,11)RMIN
   22 RETURN
      END
      SUBROUTINE TPWD(NT1,NT2)
C        SUBROUTINE TPWD FOR BMD04R                  MARCH 29, 1966
C
      IF(NT1)40,10,12
 10   NT1=5
 12   IF(NT1-NT2)14,19,14
   14 IF(NT2.EQ.5)GO TO 18
   15 REWIND NT2
   19 IF(NT1-5)18,24,18
 18   IF(NT1-6)22,40,22
 22   REWIND NT1
 24   NT2=NT1
 28   RETURN
 40   WRITE (6,49)
      CALL EXIT
      STOP
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      END
      SUBROUTINE TRANS1(X,NX,NTG,NCODE,CONS,IE)
C         SUBROUTINE TRANS1 FOR BMD04R               MARCH 29, 1966
      DIMENSION X(300),NCODE(8),CONS(8)
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
C
      DO 100 I=1,NTG
      JUMP=NCODE(I)
      C=CONS(I)
      DO 95 J=1,NX
      T=X(J)
      GO TO (1,2,3,4,5,6,7,8,9,10),JUMP
    1 IF(T)99,11,12
   11 D=0.0
      GO TO 90
   12 D=SQRT(T)
      GO TO 90
    2 IF(T) 99,13,14
   13 D=1.0
      GO TO 90
   14  D  =SQRT(T)+SQRT(T+1.0)
      GO TO 90
    3 IF(T) 99,99,15
   15 D=ALOG10(T)
      GO TO 90
    4  D  =EXP(T)
      GO TO 90
    5 IF(T) 99,11,17
   17 IF(T-1.0) 18,19,99
   18  D  =ASN(SQRT(T))
      GO TO 90
   19  D  =3.141593/2.0
      GO TO 90
    6 A=NX
      A=T/A
      IF(A) 99,20,21
   20  D  =ASN(SQRT(1.0/FLOAT(NX+1)))
      GO TO 90
   21 IF(A-1.0) 22,22,99
   22 A=NX+1
      B=T/A
       D  =ASN(SQRT(B))+ASN(SQRT(B+1.0/A))
      GO TO 90
    7 IF(T)24,99,24
   24  D  =1.0/T
      GO TO 90
    8 D=T+C
      GO TO 90
    9 D=T*C
      GO TO 90
   10 IF(T) 99,11,26
   26 D=T**C
   90 X(J)=D
   95 CONTINUE
  100  CONTINUE
      GO TO 200
 99   WRITE (6,4000)J,IE,JUMP
      GO TO 90
 4000 FORMAT(26H0THE VALUE OF SUB-INTERVALI4,14H FOR REPLICATEI3,58H VIO
     1LATED THE RESTRICTIONS FOR TRANSGENERATION OF THE TYPEI3,1H./52H T
     2HE PROGRAM CONTINUED LEAVING THIS VALUE UNCHANGED.)
  200 RETURN
      END
      SUBROUTINEVRITE
C           SUBROUTINE VRITE FOR BMD04R              MARCH 29, 1966
      COMMON  H      , N      , KNT    , T      , U      , V
      COMMON  DU     , DV     , NHARM  , NFR    , RSUM   , TSUM
      COMMON  YDATA  , SUMT2  , SUMT   , AY     , BY     , AAY
      COMMON  BBY    , FN     , ZEANY  , EQY    , CM     , S1
      COMMON  S2     , S3     , S4     , NDF2   , DF2    , NDF4
      COMMON  DF4    , NDF8   , DF8    , NDF9   , S6     , S7
      COMMON  S8     , S9     , S5     , NDF1   , DF1    , NDF3
      COMMON  DF3    , NDF5   , DF5    , NDF6   , DF6    , NDF7
      COMMON  DF7    , AAX    , BBX    , AX     , BX     , RSUMX
      COMMON  TSUMX  , SUMT2X , SUMTX  , EQX    , SUMXY  , XDATA
      COMMON  S4XY   , S2XY   , S3XY   , S5XY   , CMX    , S1X
      COMMON  S2X    , S3X    , S4X    , S5X    , S1XY   , EXX1
      COMMON  EXX2   , EXX3   , EXY1   , EXY2   , EXY3   , EXY4
      COMMON  NDFE2  , NDFE4  , NDFE3  , BE2    , BS12   , BS22
      COMMON  BS32   , BS42   , NDFERR , DFERR  , BE2M   , BS12M
      COMMON  BS22M  , BS32M  , BS42M  , FC1    , FC2    , FC3
      COMMON  FC4    , BYX    , RANGE  , THETA  , RMAX   , RMIN
      COMMON  S1M    , S2M    , S3M    , S4M    , S5M    , S8M
      COMMON  S9M    , FV1    , FV3    , FV4    , FV2    , ESTD
      COMMON  FV8    , FV9    , DFESTC , FC5    , DF9    , TIME
      COMMON  SEAMP  , EXX4   , ADDS8  , FOR
      DIMENSIONU(300,9),V(300,9),DU(9),DV(9),RSUM(20),TSUM(300),YDATA(15
     1,300),AY(20,9),BY(20,9),AAY(9),BBY(9),EQY(300),S2(9),S4(9),NDF2(9)
     2,DF2(9),NDF4(9),DF4(9),NDF8(9),DF8(9),NDF9(9),DF9(9),S8(9),S9(9),A
     3AX(9),BBX(9),AX(20,9),BX(20,9),RSUMX(20),TSUMX(300),EQX(300),XDATA
     4( 300  ),S2XY(9),S4XY(9),S2X(9),S4X(9),EXX2(9),EXX4(9),EXY2(9),EXY
     54(9),NDFE2(9),NDFE4(9),BS22(9),BS42(9),BS22M(9),BS42M(9),FC2(9),FC
     64(9),S2M(9),S4M(9),S8M(9),S9M(9),FV4(9),FV2(9),FV8(9),FV9(9),FMT(1
     72),DFESTC(9),ESTD(9),ADDS8(9),FOR(12)
      DOUBLE PRECISION TODE,CODE,SPC,ODE,CODEX
C
  509 FORMAT(1H134X20HANALYSIS OF VARIANCE//)
  510 FORMAT(55H SOURCE OF VARIATION   DF   SUM OF SQUARES  MEAN SQUARE6
     1X,25HF RATIOS   DF FOR F TESTS)
  511 FORMAT(20H BETWEEN REPLICATES I5,F15.4,F13.4,F14.4,5H    (I4,1H,I5
     1,1H))
  512 FORMAT(13H EFFECT OF (AI1,2H+BI1,3H)  I5,F15.4,F13.4,F14.4,5H    (
     1I4,1H,F8.2,1H))
  513 FORMAT(20H SCATTER ABOUT CURVEI5,F15.4,F13.4,F14.4,5H    (I4,1H,I5
     1,1H))
  514 FORMAT(15H REPLICATE X (AI1,2H+BI1,1H)I5,F15.4,F13.4,F14.4,5H    (
     1I4,1H,I5,1H))
  515 FORMAT(20H REPLICATE X SCATTERI5,F15.4,F13.4)
  516 FORMAT(1H 85(1H*))
  517 FORMAT(6H TOTAL14XI5,F15.4)
  518 FORMAT(20H CORRECTION TO MEAN I5,F15.4)
  519 FORMAT(19H REPLICATE X AMPLI I1,I5,F15.4,F13.4,F14.4,5H    (I4,1H,
     1I5,1H))
  520 FORMAT(19H REPLICATE X PHASE I1,I5,F15.4,F13.4,F14.4,5H    (I4,1H,
     1I5,1H))
      WRITE (6,509)
      WRITE (6,510)
      WRITE (6,516)
      WRITE (6,511)NDF1,S1,S1M,FV1,NDF1,NDF5
      DO3I=1,NHARM
    3 WRITE (6,512)I,I,NDF2(I),S2(I),S2M(I),FV2(I),NDF2(I), ESTD(I)
      WRITE (6,513)NDF3,S3,S3M,FV3,NDF3,NDF5
      DO4I=1,NHARM
    4 WRITE (6,514)I,I,NDF4(I),S4(I),S4M(I),FV4(I),NDF4(I),NDF5
      WRITE (6,515)NDF5,S5,S5M
      WRITE (6,516)
      WRITE (6,517)NDF6,S6
      WRITE (6,518)NDF7,S7
      WRITE (6,516)
      DO5I=1,NHARM
      WRITE (6,519)I,NDF8(I),S8(I),S8M(I),FV8(I),NDF8(I),NDF5
    5 WRITE (6,520)I,NDF9(I),S9(I),S9M(I),FV9(I),NDF8(I),NDF5
      WRITE (6,516)
      RETURN
      END

















































