C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAU5 MYELOSUPPRESSION (WOJCIECH KRZYZANSKI)
C --- THIS IS EXAMPLE 3 OF "APPLYING STIFF ... DISTRIBUTED DELAYS"
C     version of 27 May 2025
C * * * * * * * * * * * * * * * * * * * * * * * * *
        include '../radau5.f'
        include '../decsol.f'
        include 'dc_sumexp.f'
        IMPLICIT REAL*8 (A-H,O-Z)
C --- PARAMETERS FOR RADAU5
        INTEGER, PARAMETER :: NDM=5000
        PARAMETER (LWORK=4*NDM*NDM+12*NDM+20,LIWORK=3*NDM+20)               
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(NDM),WORK(LWORK),IWORK(LIWORK),YEX(NDM)
        DIMENSION ATOL(NDM),RTOL(NDM)
        EXTERNAL FGAM1,JGAM1,JGAMSX1,FGAM2,JGAM2,JGAMSX2,SOLOUT,MAS
        REAL*4 start,finish
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT
        COMMON/EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha        
C
C --- CHOICE BETWEEN SETS OF PARAMETERS (IPROB = 1 0R 2)  
        IPROB = 2
        CALL PROBDATA(IPROB)
        IF (IPROB.EQ.1) MI=1
        IF (IPROB.EQ.2) MI=2
C --- FOR NEGATIVE ALPHA, AUGMENT IT BY 1 (mi = 2)
        ialph = 0
        IF (IPROB.EQ.2) ialph = 1
C --- DESIRED ACCURACY
        TOL=1.D-7
        EPSIL=TOL
C --- END OF INTEGRATION
        XEND=100.D0
        CALL HMNEXP (XEND,kappa,alpha,ialph,epsil,HEXP,MEXP,NEXP)
        WRITE(*,*) 'HEXP, MEXP, NEXP ', HEXP, MEXP, NEXP
        RPAR(6)=HEXP
        IPAR(1)=MEXP
        IPAR(2)=NEXP
C --- FULL LINEAR ALGEBRA (ILIN = 1) AND ADAPTED TO SUMEXP (ILIN = 2)
        DO ILIN = 1,2
          ND=4       
          IF (ILIN.EQ.1) N=ND+MI*(NEXP-MEXP)
          IF (ILIN.EQ.2) N=ND+2+MI*(NEXP-MEXP)
C --- ODE FORMULATION (IDAE = 0) OR DAE FORMULATION (IDAE = 1)
        DO IDAE = 0,1
          write (6,*) '  '
          IPAR(3)=IDAE
          X=0.0D0 
          Y=0.0D0
          Y(1)=(Kc/K)*W0
          Y(2)=W0
          Y(3)=A0
C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM? NO
          IMAS=1
          MLMAS=0
          MUMAS=0
C --- OUTPUT ROUTINE IS USED DURING INTEGRATION
          IOUT=0
C --- REQUIRED TOLERANCE
          RTOL=TOL
          ATOL=RTOL
          ITOL=1
          FACTOL = 1.0D2
          DO I=ND+1,N
            ATOL(I)=TOL*FACTOL
            RTOL(I)=TOL*FACTOL
          END DO
          FACTOL = 1.0D-2
          ATOL(ND)=TOL*FACTOL
          RTOL(ND)=TOL*FACTOL
          H=MAX(TOL,1.0D-5)
C --- SET DEFAULT VALUES 
          DO I=1,20
            IWORK(I)=0
            WORK(I)=0.D0
          END DO
          IF (ILIN.EQ.1) THEN
            IJAC=1
            MLJAC=N
C --- CALL OF THE SUBROUTINE RADAU5
            IF (MI.EQ.1) THEN
              CALL cpu_time(start)
              CALL RADAU5(N,FGAM1,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAM1,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
              CALL cpu_time(finish)
            END IF
            IF (MI.EQ.2) THEN
              CALL cpu_time(start)
              CALL RADAU5(N,FGAM2,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAM2,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
              CALL cpu_time(finish)
            END IF
          END IF
          IF (ILIN.EQ.2) THEN
            IJAC=1
            MUJAC=0
            MLJAC=ND
C --- CALL OF THE SUBROUTINE RADAU5
            IF (MI.EQ.1) THEN
              CALL cpu_time(start)
              CALL RADAU5(N,FGAM1,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAMSX1,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
              CALL cpu_time(finish)
            END IF
            IF (MI.EQ.2) THEN
              CALL cpu_time(start)
              CALL RADAU5(N,FGAM2,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAMSX2,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
              CALL cpu_time(finish)
            END IF
          END IF
        IF (IPROB.EQ.1) THEN 
          YEX(1)=1.9067548914832040D0
          YEX(2)=6.3753272029666723D0
          YEX(3)=0.0D0
        END IF
         IF (IPROB.EQ.2) THEN 
          YEX(1)= 2.5190147371973723D0
          YEX(2)= 6.5179975042328913D0
          YEX(3)=0.0D0
        END IF
          ERR=0.0D0
          DO I=1,2
            ERR = ERR + ((Y(I)-YEX(I))/YEX(I))**2
          END DO
          ERR=SQRT(ERR)
C --- PRINT FINAL SOLUTION
          write (6,77) ILIN,IDAE,TOL,finish-start
 77       FORMAT(1X,'radau5, lin-alg = ',I1,'   IDAE = ',I1,
     &           '   TOL = ',E9.3,'    time = ',E8.2)
          WRITE (6,99) X,(Y(I),I=1,2),ERR
 99       FORMAT(1X,'X =',F10.3,'    Y =',2E20.12,'    ERR = ',E9.3)
          WRITE (6,91) (IWORK(J),J=14,20)
 91       FORMAT(' fcn=',I5,' jac=',I4,' step=',I4,' accpt=',I4,
     &        ' rejct=',I3,' dec=',I4,' sol=',I5)
        END DO
        END DO
        STOP
        END
C
        SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN)
C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5"
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),CONT(LRC)
        COMMON /INTERN/XOUT
        WRITE (9,99) X,Y(1),Y(2),Y(3), NR-1         
        IF (NR.EQ.1) THEN
           XOUT=0.0D0
        ELSE
 10        CONTINUE
           IF (X.GE.XOUT) THEN
              WRITE (6,99) XOUT,CONTR5(1,XOUT,CONT,LRC),
     1                          CONTR5(2,XOUT,CONT,LRC),
     2                          CONTR5(3,XOUT,CONT,LRC), NR-1        
              XOUT=XOUT+1.0D0
              GOTO 10
           END IF
        END IF
 99     FORMAT(1X,'X =',F5.1,'    Y =',3E18.10,'    NSTEP = ',I4)
        RETURN
        END
C
        SUBROUTINE FGAM1(N,X,Y,F,RPAR,IPAR)
C --- RIGHT-HAND SIDE OF EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),F(N)
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT
        COMMON /EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha        
        PI=4*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        IDAE=IPAR(3)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        APPINT=0.D0
        GINC=0.0D0
        II=5
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI = DEXP(IND*HEXP) + kappa
          CI     = DEXP(alpha*IND*HEXP)
          F(II)  = - GAMMAI*Y(II) + Y(1)
          APPINT =  APPINT + CI*Y(II)
          GINC =  GINC + CI*(1-EXP(-GAMMAI*X))/GAMMAI
          II=II+1
        END DO
        FACT=HEXP*kappa**(1.D0-alpha)*SIN(PI*alpha)/PI
        GINC=GINC*FACT
        APPINT=APPINT*FACT
c
c        GINC = gammds ( kappa*X, 1.D0-alpha, ifault )
c
        C = Y(3)/V
        F(1) = (Kp*(W0/Y(2))**gamv- Ks*C - K)*Y(1)
        GINC = 1.D0-GINC
        C4=GINC*(Kc/K)*W0
        F(2) = K*(Y(4)+C4) - Kc*Y(2)
        IF (IDAE.EQ.0) F(3) = -Vmax*Y(3)/(Km+C)
        IF (IDAE.EQ.1) THEN
          X0 = 0.0D0
          FA = - (Y(3)-A0)/(V*Km) - Vmax*(X-X0)/Km
          F(3) = Y(3) - A0*EXP(FA)
        END IF
        F(4) = APPINT - Y(4)
        RETURN
        END  
C
        SUBROUTINE JGAM1(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN OF THE EQUATION (FULL OPTION)
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),DFY(LDFY,N)
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT
        COMMON /EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha
        PI=4*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        IDAE=IPAR(3)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        C = Y(3)/V
        DFY=0.0D0
        DFY(1,1)=Kp*(W0/Y(2))**gamv- Ks*C - K
        DFY(1,2)=-gamv*Kp*Y(1)*(W0/Y(2))**gamv/Y(2)
        DFY(1,3)=-Ks*Y(1)/V
        DFY(2,2)=-Kc
        DFY(2,4)=K
        IF (IDAE.EQ.0) DFY(3,3)=-Km*Vmax/(Km+C)**2
        IF (IDAE.EQ.1) THEN
          X0 = 0.0D0
          FA = - (Y(3)-A0)/(V*Km) - Vmax*(X-X0)/Km
          DFY(3,3)=1.0D0 + A0*EXP(FA)/(V*Km)
        END IF
        DFY(4,4)=-1.0D0
        II=5
        FAC=HEXP*kappa**(1.D0-alpha)*SIN(PI*alpha)/PI
        DO I=1,LEXP
          IND = MEXP + (I-1)
          GAMMAI     = DEXP(IND*HEXP)+kappa
          CI         = DEXP(alpha*IND*HEXP)
          DFY(4,II)  = FAC*CI
          DFY(II,1)  = 1.0D0
          DFY(II,II) = - GAMMAI
          II=II+1
        END DO
        RETURN
        RETURN
        END
C
        SUBROUTINE JGAMSX1(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN OF THE EQUATION (SUMEXP OPTION)
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR     
        DIMENSION Y(N),DFY(LDFY,N)
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT
        COMMON /EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha
        PI=4*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        IDAE=IPAR(3)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        C = Y(3)/V
        DFY = 0.0D0
        DFY(1,1)=Kp*(W0/Y(2))**gamv- Ks*C - K
        DFY(1,2)=-gamv*Kp*Y(1)*(W0/Y(2))**gamv/Y(2)
        DFY(1,3)=-Ks*Y(1)/V
        DFY(2,2)=-Kc
        DFY(2,4)=K
        IF (IDAE.EQ.0) DFY(3,3)=-Km*Vmax/(Km+C)**2
        IF (IDAE.EQ.1) THEN
          X0 = 0.0D0
          FA = - (Y(3)-A0)/(V*Km) - Vmax*(X-X0)/Km
          DFY(3,3)=1.0D0 + A0*EXP(FA)/(V*Km)
        END IF
        DFY(4,4)=-1.0D0
        FAC=HEXP*kappa**(1.D0-alpha)*SIN(PI*alpha)/PI
        DFY(4,5)=1.0D0
        DFY(1,6)=1.0D0
        II=7
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI =  DEXP(IND*HEXP)+kappa
          CI =      DEXP(alpha*IND*HEXP)
          DFY(1,II) = CI*FAC
          DFY(2,II) = - GAMMAI
          DFY(3,II) = 0.0D0
          II=II+1
        END DO
        RETURN
        END
C
        SUBROUTINE FGAM2(N,X,Y,F,RPAR,IPAR)
C --- RIGHT-HAND SIDE OF EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),F(N)
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT
        COMMON /EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha        
        PI=4*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        IDAE=IPAR(3)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        APPINT=0.D0
        II=5 
        GINC=0.0D0
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI  = DEXP(IND*HEXP)+kappa
          CI      = DEXP(alpha*IND*HEXP)
          F(II)   = - GAMMAI*Y(II) + Y(1)
          F(II+1) = - GAMMAI*Y(II+1) + Y(II)
          APPINT  = APPINT+CI*Y(II+1)
          GINC =  GINC + CI*(1-EXP(-GAMMAI*X)*(1+GAMMAI*X))/GAMMAI**2
          II=II+2
        END DO
        FACT=HEXP*kappa**(2.D0-alpha)*SIN(PI*alpha)/(PI*(1-alpha))
        GINC=GINC*FACT
        APPINT=APPINT*FACT
c
c        GINC = gammds ( kappa*X, 1.D0-alpha0, ifault )
c
        C = Y(3)/V
        F(1) = (Kp*(W0/Y(2))**gamv - Ks*C - K)*Y(1)
        GINC = 1.D0-GINC
        C4=GINC*(Kc/K)*W0
        F(2) = K*(Y(4)+C4) - Kc*Y(2)
        IF (IDAE.EQ.0) F(3) = -Vmax*Y(3)/(Km+C)
        IF (IDAE.EQ.1) THEN
          X0 = 0.0D0
          FA = - (Y(3)-A0)/(V*Km) - Vmax*(X-X0)/Km
          F(3) = Y(3) - A0*EXP(FA)
        END IF  
        F(4) = APPINT - Y(4)
        RETURN
        END  
C        
        SUBROUTINE JGAM2(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN OF THE EQUATION (FULL OPTION)
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR
        DIMENSION Y(N),DFY(LDFY,N)
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT
        COMMON /EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha
        PI=4*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        IDAE=IPAR(3)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        C = Y(3)/V
        DFY=0.0D0
        DFY(1,1)=Kp*(W0/Y(2))**gamv - Ks*C - K
        DFY(1,2)=-gamv*Kp*Y(1)*(W0/Y(2))**gamv/Y(2)
        DFY(1,3)=-Ks*Y(1)/V
        DFY(2,2)=-Kc
        DFY(2,4)=K
        IF (IDAE.EQ.0) DFY(3,3)=-Km*Vmax/(Km+C)**2
        IF (IDAE.EQ.1) THEN
          X0 = 0.0D0
          FA = - (Y(3)-A0)/(V*Km) - Vmax*(X-X0)/Km
          DFY(3,3)=1.0D0 + A0*EXP(FA)/(V*Km)
        END IF
        DFY(4,4)=-1.0D0
        FAC=HEXP*kappa**(2.D0-alpha)*SIN(PI*alpha)/(PI*(1-alpha))
        II=5
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI =  DEXP(IND*HEXP)+kappa
          CI =      DEXP(alpha*IND*HEXP)
          DFY(4,II+1) = FAC*CI
          DFY(II,1) = 1.0D0
          DFY(II,II) = - GAMMAI
          DFY(II+1,II+1) = - GAMMAI
          DFY(II+1,II) = 1.0D0
          II=II+2
        END DO
        RETURN
        END
C
        SUBROUTINE JGAMSX2(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN OF THE EQUATION (SUMEXP OPTION)
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR     
        DIMENSION Y(N),DFY(LDFY,N)
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT
        COMMON /EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha
        PI=4*ATAN(1.0D0)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        IDAE=IPAR(3)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        C = Y(3)/V
        DFY = 0.0D0
        DFY(1,1)=Kp*(W0/Y(2))**gamv- Ks*C - K
        DFY(1,2)=-gamv*Kp*Y(1)*(W0/Y(2))**gamv/Y(2)
        DFY(1,3)=-Ks*Y(1)/V
        DFY(2,2)=-Kc
        DFY(2,4)=K
        IF (IDAE.EQ.0) DFY(3,3)=-Km*Vmax/(Km+C)**2
        IF (IDAE.EQ.1) THEN
          X0 = 0.0D0
          FA = - (Y(3)-A0)/(V*Km) - Vmax*(X-X0)/Km
          DFY(3,3)=1.0D0 + A0*EXP(FA)/(V*Km)
        END IF
        DFY(4,4)=-1.0D0
        FAC=HEXP*kappa**(2.D0-alpha)*SIN(PI*alpha)/(PI*(1-alpha))
        DFY(4,5)=1.0D0
        DFY(1,6)=1.0D0
        II=7
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI =  DEXP(IND*HEXP)+kappa
          CI =      DEXP(alpha*IND*HEXP)
          DFY(1,II) = 0.0D0
          DFY(1,II+1) = CI*FAC
          DFY(2,II) = - GAMMAI
          DFY(2,II+1) = - GAMMAI
          DFY(3,II) = 1.0D0
          DFY(3,II+1) = 0.0D0
          II=II+2
        END DO
        RETURN
        END
C
        SUBROUTINE PROBDATA(IPROB)
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8 kappa,K,Ks,Kp,Km,Kc,MTT,NU
        COMMON /EQPAR/kappa,K,Ks,Kp,Km,Kc,gamv,Vmax,V,MTT,A0,W0,alpha
c
        IF(IPROB.EQ.1) THEN
          alpha = 0.036D0
          NU = 1.D0-alpha 
          MTT = 47.5D0
          K = NU/MTT
          kappa=K
          Kp = K
          Kc = K
          gamv = 0.664D0
          Ks = 0.0328D0
          Vmax = 77.2D0
          Km = 16.9D0
          V = 1.35D0
          W0 = 14.4D0
          A0 = 127.D0   
        END IF     
c
        IF(IPROB.EQ.2) THEN
          alpha = -0.46D0
          NU = 1.D0-alpha 
          MTT = 55.6D0
          K = NU/MTT
          kappa=K
          Kp = K
          Kc = K
          gamv = 0.507D0
          Ks = 0.0213D0
          Vmax = 100.0D0
          Km = 22.0D0
          V = 1.03D0
          W0 = 14.4D0
          A0 = 127.D0   
        END IF     
c
         RETURN
         END
C
         SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR)
         IMPLICIT REAL*8 (A-H,O-Z)
         DIMENSION AM(LMAS,N)
         REAL*8, dimension(20) :: RPAR
         INTEGER, dimension(20) :: IPAR
         IDAE = IPAR(3)
         DO I=1,N
           AM(1,I)=1.0D0
         END DO
         IF (IDAE.EQ.1) AM(1,3) = 0.0D0
         AM(1,4) = 0.0D0
         RETURN
         END
C        
        SUBROUTINE HMNEXP (XEND,kappa,alpha,ialph,epsil,HEXP,MEXP,NEXP)
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8 kappa,alpha,epsil
C --- COMPUTES THE VALUES OF HEXP,MEXP,NEXP FOR THE GAMMA DISTRIBUTION
        Tl=1.0D-15
        Tr=XEND
        Vl=((kappa*Tl)**(-alpha)*exp(-kappa*Tl)/gamma(1-alpha))/epsil
        Vr=((kappa*Tr)**(-alpha)*exp(-kappa*Tr)/gamma(1-alpha))/epsil
  33    CONTINUE
        T=(Tr+Tl)/2
        VT=((kappa*T)**(-alpha)*exp(-kappa*T)/gamma(1-alpha))/epsil
        IF (VT.LE.1.D0) THEN
            Tr=T
        ELSE
            Tl=T
        END IF
        T=(Tl+Tr)/2
        IF (ABS(Tr-Tl).LE.1.D-15) GOTO 33
        T=MIN(XEND,T)
        deltlog = LOG(GAMMA(2.D0-alpha)*epsil)/(1.D0-alpha)
C ---
        alpha = alpha + ialph
C ---
        pi=4.D0*ATAN(1.D0)
        a=(PI/2.D0)*(1.D0+alpha/((alpha+1.D0)*LOG(epsil)));
        HEXP=2.D0*PI*a/LOG(1.D0+(2.D0/epsil)*COS(a)**(-alpha));
        xllog = LOG(GAMMA(alpha+1.D0)*epsil)/alpha
        xs = -LOG(GAMMA(alpha)*epsil)
        NEXP = INT((LOG(xs)-deltlog+LOG(kappa))/HEXP) + 1
        MEXP = INT((xllog - LOG(T))/HEXP)-1
        RETURN
        END

