      DIMENSION T(251),C(251),TX(250),CY(250)
      CHARACTER XS*10,YS*10,TN*14,MA*1
      NS=50
      WRITE(*,*) '###########################################'
      WRITE(*,*) '#   MOMENT CALCULATION FOR PLASMA LEVEL   #'
      WRITE(*,*) '#   COPYRIGHT MARCH 22(1996) BY K.YAMAOKA #'
      WRITE(*,*) '#   FACULTY OF PHARMACEUTICAL SCIENCES    #'
      WRITE(*,*) '#   KYOTO UNIVERSITY                      #'
      WRITE(*,*) '###########################################'
      WRITE(*,*)
      CALL FMA1(N,TX,CY,TN)
      N1=N+1
      DO 10 I=1,N
      T(I+1)=TX(I)
      C(I+1)=CY(I)
   10 CONTINUE
      T0=0.0
      T9=-1E+10
      Y0=0.0
      Y9=-1E+10
      DO 30 I=1,N
      IF(T9.LT.TX(I)) T9=TX(I)
      IF(Y9.LT.CY(I)) Y9=CY(I)
   30 CONTINUE
      XS='TIME'
      YS='CP'
      WRITE(*,*)
      CALL PLOT2(0,0,XT,YT,T0,T9,Y0,Y9,XS,YS,N,TX,CY)
      CALL TMP
      IF(C(2).EQ.0.0) THEN
      DO 35 I=3,N1
      TX(I-2)=T(I)
      CY(I-2)=C(I)
   35 CONTINUE
      DO 36 I=2,N1
      T(I-1)=T(I)
      C(I-1)=C(I)
   36 CONTINUE
      N1=N1-1
      N=N-1
      GOTO 37
      ENDIF
      A1=(C(3)-C(2))/(T(3)-T(2))
      B1=C(2)-A1*T(2)
      IF(A1.GT.0.0) THEN
      T(1)=-B1/A1
      IF(T(1).LT.0.0) T(1)=0.0
      C(1)=0.0
      WRITE(*,'('' LAG TIME='',F10.5)') T(1)
      ELSE
      T(1)=0.0
      C(1)=B1
      WRITE(*,'('' INITAL CONCENTRATION='',F10.5)') C(1)
      END IF
   37 S0=0.0
      S1=0.0
      S2=0.0
      DO 20 I=2,N1
      S0=S0+0.5*(C(I)+C(I-1))*(T(I)-T(I-1))
      S1=S1+0.5*(T(I)*C(I)+T(I-1)*C(I-1))*(T(I)-T(I-1))
      S2=S2+0.5*(T(I)*T(I)*C(I)
     &+T(I-1)*T(I-1)*C(I-1))*(T(I)-T(I-1))
   20 CONTINUE
      AUC=S0
      RMT=S1/S0
      VRT=S2/S0-RMT*RMT
      WRITE(*,'('' **  '',A,''<< FOR FINITE TIME>> **'')') TN
      WRITE(*,'('' AUC='',F10.5)') AUC
      WRITE(*,'('' MRT='',F10.5)') RMT
      WRITE(*,'('' VRT='',F10.5)') VRT
      CALL TMP
      Y0=1E+10
      Y9=-1E+10
      DO 40 I=1,N
      CY(I)=ALOG(CY(I))
      IF(Y0.GT.CY(I)) Y0=CY(I)
      IF(Y9.LT.CY(I)) Y9=CY(I)
   40 CONTINUE
      YS='LOG CP'
      WRITE(*,*)
      CALL PLOT2(0,0,XT,YT,T0,T9,Y0,Y9,XS,YS,N,TX,CY)
   47 WRITE(*,*) 'NUMBER FOR EXTRAPOLATION(0:STOP,1:AUTO)?'
      READ(*,*) NP
      IF(NP.EQ.0) STOP
      IF(NP.EQ.1) GOTO 46
      DO 42 I=1,NP
      TX(I)=T(N1+I-NP)
      CY(I)=ALOG(C(N1+I-NP))
   42 CONTINUE
      CALL LINE(NP,TX,CY,A,B)
      AA=A
      BB=B
      NN=NP
      B=EXP(BB)
      A=-AA
      P0=S0+B/A*EXP(-A*T(N1))
      P1=S1+(B/A/A+B*T(N1)/A)*EXP(-A*T(N1))
      P2=S2+(2*B/A/A/A+2*B*T(N1)/A/A+B*T(N1)*T(N1)/A)*EXP(-A*T(N1))
      WRITE(*,*)
      WRITE(*,'('' **  '',A,''<< FOR INFINITE TIME>> **'')') TN
      WRITE(*,'('' ** EXTRAPOLATION FUNCTION NP='',I3,''**'')') NN
      WRITE(*,'(1X,F10.5,''*EXP(-'',F10.5,''*T)'')') B,A
      AUC=P0
      RMT=P1/P0
      VRT=P2/P0-RMT*RMT
      WRITE(*,'('' AUC='',F10.5)') AUC
      WRITE(*,'('' MRT='',F10.5)') RMT
      WRITE(*,'('' VRT='',F10.5)') VRT
      GOTO 47
   46 CM=-1E+10
      DO 52 I=1,N1
      IF(CM.GE.C(I)) GOTO 52
      IM=I
      CM=C(I)
   52 CONTINUE
      IM=IM-1
      IF(IM.GT.N1-3) THEN
      WRITE(*,*) '** NOT ENOUGH DATA **'
      GOTO 100
      ENDIF
      AM=1E+10
      DO 70 J=IM,N1-3
      NP=N1-J
      DO 60 I=1,NP
      TX(I)=T(N1+I-NP)
      CY(I)=ALOG(C(N1+I-NP))
   60 CONTINUE
      CALL LINE(NP,TX,CY,A,B)
      SS=0.0
      DO 62 I=1,NP
      SS=SS+(CY(I)-A*TX(I)-B)*(CY(I)-A*TX(I)-B)
   62 CONTINUE
      IF(SS.EQ.0.0) THEN
      AA=A
      BB=B
      NN=NP
      GOTO 72
      END IF
      AIC=FLOAT(NP)*ALOG(SS)
      IF(AM.LE.AIC) GOTO 70
      AM=AIC
      AA=A
      BB=B
      NN=NP
   70 CONTINUE
   72 B=EXP(BB)
      A=-AA
      P0=S0+B/A*EXP(-A*T(N1))
      P1=S1+(B/A/A+B*T(N1)/A)*EXP(-A*T(N1))
      P2=S2+(2*B/A/A/A+2*B*T(N1)/A/A+B*T(N1)*T(N1)/A)*EXP(-A*T(N1))
      WRITE(*,*)
      WRITE(*,'('' **  '',A,''<< FOR INFINITE TIME>> **'')') TN
      WRITE(*,'('' ** EXTRAPOLATION FUNCTION NP='',I3,''**'')') NN
      WRITE(*,'(1X,F10.5,''*EXP(-'',F10.5,''*T)'')') B,A
      AUC=P0
      RMT=P1/P0
      VRT=P2/P0-RMT*RMT
      WRITE(*,'('' AUC='',F10.5)') AUC
      WRITE(*,'('' MRT='',F10.5)') RMT
      WRITE(*,'('' VRT='',F10.5)') VRT
      GOTO 47 
  100 END
      SUBROUTINE PLOT2(NS,NE,XT,YT,X0,X9,Y0,Y9,XS,YS,NP,TX,CY)
      DIMENSION XT(50),YT(50,5),XG(6),TX(250),CY(250)
      CHARACTER XS*10,YS*10,IP(51,17)
      XM=50.0
      YM=16.0
      DO 100 I=1,51
      DO 100 J=1,17
      IP(I,J)=' '
  100 CONTINUE
      DO 110 I=1,NS
      IX=INT(XM*(XT(I)-X0)/(X9-X0)+0.5)+1
      DO 110 J=1,NE
      IY=INT(YM*(Y9-YT(I,J))/(Y9-Y0)+0.5)+1
      IP(IX,IY)=CHAR(J+48)
  110 CONTINUE
      DO 107 I=1,NP
      IX=INT(XM*(TX(I)-X0)/(X9-X0)+0.5)+1
      IY=INT(YM*(Y9-CY(I))/(Y9-Y0)+0.5)+1
      IP(IX,IY)=CHAR(79)
  107 CONTINUE
      WRITE(*,'(5X,A10)') YS
      DO 130 K=1,17
      JJ=K-4*INT(K/4)
      IF(JJ.EQ.1) THEN
      YY=(Y9-Y0)*(16-K+1)/16.0+Y0
      IF((YY.GE.10000.0).OR.((YY.LT.0.1).AND.(YY.NE.0.0))) THEN
      WRITE(*,'(1X,E12.4,''+'',51A1)') YY,(IP(I,K),I=1,51)
      ELSE
      WRITE(*,'(5X,F8.3,''+'',51A1)') YY,(IP(I,K),I=1,51)
      END IF
      ELSE
      WRITE(*,'(13X,''I'',51A1)') (IP(I,K),I=1,51)
      END IF
  130 CONTINUE
      WRITE(*,'(13X,''++'',5A10)') ('---------+',I=1,5)
      DO 140 I=1,6
      XG(I)=(X9-X0)*(I-1)/5.0+X0
  140 CONTINUE
      WRITE(*,'(7X,6F10.2)') (XG(I),I=1,6)
      WRITE(*,'(35X,A10)') XS
      END
      SUBROUTINE TMP
      CHARACTER TP
      WRITE(*,*) 'DO YOU CONTINUE(Y/N)?'
      READ(*,'(A)') TP
      IF(TP.EQ.'N') STOP
      END
      SUBROUTINE LINE(N,X,Y,A,B)
      DIMENSION X(50),Y(50)
      XS=0.0
      YS=0.0
      XX=0.0
      YY=0.0
      XY=0.0
      DO 10 I=1,N
      XS=XS+X(I)
      YS=YS+Y(I)
      XX=XX+X(I)*X(I)
      YY=YY+Y(I)*Y(I)
      XY=XY+X(I)*Y(I)
   10 CONTINUE
      XS=XS/FLOAT(N)
      YS=YS/FLOAT(N)
      XY=XY/FLOAT(N)
      XX=XX/FLOAT(N)
      YY=YY/FLOAT(N)
      A=(XY-XS*YS)/(XX-XS*XS)
      B=YS-A*XS
      END
C   ====== FILE I/O ======= 
      SUBROUTINE FMA1(N,X,Y,FN)
      DIMENSION X(50),Y(50)
      CHARACTER FN*14,IC*1
  100 WRITE(*,*) 'E)XIT, I)NPUT, R)EAD, W)RITE, S)TOP ?'
      READ(*,'(A)') IC
      IF(IC.EQ.'S') STOP
      IF(IC.EQ.'E'.AND.N.GT.0) RETURN
      IF(IC.EQ.'I') GOTO 1100
      IF(IC.EQ.'R') GOTO 1200
      IF(IC.EQ.'W') GOTO 1300
      GOTO 100
 1100 WRITE(*,'('' NUMBER OF POINTS ?'')')
      READ(*,*) N
      DO 1127 I=1,N
      WRITE(*,1125) I,I
 1125 FORMAT(' T',I2,', CP',I2,'?')
      READ(*,*) X(I),Y(I)
 1127 CONTINUE
      GOTO 100
 1200 WRITE(*,*) 'FILE NAME TO READ ?'
      READ(*,'(A)') FN
      OPEN(UNIT=5,FILE=FN)
      READ(5,*) N
      WRITE(*,*)
      DO 1210 I=1,N
      READ(5,*) X(I),Y(I)
      WRITE(*,1225) I,X(I),I,Y(I)
 1225 FORMAT(' T',I2,'=',F10.5,', CP',I2,'=',F10.5)
 1210 CONTINUE
      CLOSE(UNIT=5)
      GOTO 100
 1300 WRITE(*,*) 'FILE NAME TO WRITE ?'
      READ(*,'(A)') FN
      OPEN(UNIT=6,FILE=FN,STATUS='NEW')
      WRITE(6,*) N
      DO 1327 I=1,N
      WRITE(6,'(X,F12.5,'','',F12.5)')X(I),Y(I)
 1327 CONTINUE 
      CLOSE(UNIT=6)
      GOTO 100
      END
      
