C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAU5: EXAMPLE 1 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 FGAMI,JGAMI,JGAMIEX,MAS,SOLOUT
        REAL*4 start,finish
        REAL*8 kappa
C --- PARAMETERS OF THE PROBLEM
        alpha=0.5D0
c        alpha= 0.1D0
        ialph=0
        kappa = 0.25D0
        RPAR(1)=alpha
        RPAR(2)=kappa
C --- DESIRED ACCURACY
        TOL=1.0D-9
        EPSIL=TOL
C --- END OF INTEGRATION
        XEND=50.D0
        CALL HMNEXP (XEND,kappa,alpha,ialph,epsil,HEXP,MEXP,NEXP)
        write (6,*) '  '
        WRITE(*,*) 'HEXP,MEXP,NEXP ', HEXP, MEXP, NEXP
        RPAR(6)=HEXP
        IPAR(1)=MEXP
        IPAR(2)=NEXP
C --- SAME TOL FOR ALL COMNPONENTS (ITOL = 0) COMPONENTWISE TOL (ITOL = 1)
        DO ITOL = 0,0
C --- FULL LINEAR ALGEBRA (ILIN = 1) AND ADAPTED TO SUMEXP (ILIN = 2)
        DO ILIN = 1,2
           write (6,*) '  '
           MI=1
           ND=2
           N=ND+MI*(NEXP-MEXP)
           IF (ILIN.EQ.2) N=N+2
C --- INITIAL VALUES
           X=0.0D0 
           Y=0.0D0
C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM? NO
           IMAS=1
           MLMAS=0
           MUMAS=0
C --- OUTPUT ROUTINE IS USED DURING INTEGRATION? NO
           IOUT=0
C --- REQUIRED TOLERANCE
           RTOL=TOL
           ATOL=TOL
           IF (ITOL.EQ.1) THEN
             FACTOL = SQRT(TOL)/10.0D0
             DO I=ND+1,N
               ATOL(I)=FACTOL
               RTOL(I)=FACTOL
             END DO
           END IF
           H=0.1D0
C --- SET DEFAULT VALUES 
           DO I=1,20
             IWORK(I)=0
             WORK(I)=0.D0
           END DO
           IF (ILIN.EQ.1) THEN
             IJAC=1
c             IJAC=0
             MLJAC=N
C --- CALL OF THE SUBROUTINE RADAU5
             CALL cpu_time(start)
             CALL RADAU5(N,FGAMI,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAMI,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
             CALL cpu_time(finish)
           END IF
           IF (ILIN.EQ.2) THEN
             IJAC=1
             MUJAC=0
             MLJAC=2    !  1 y-component + 1 integral
C --- CALL OF THE SUBROUTINE RADAU5
             CALL cpu_time(start)
             CALL RADAU5(N,FGAMI,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JGAMIEX,IJAC,MLJAC,MUJAC,
     &                  MAS,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
             CALL cpu_time(finish)
           END IF
           YEX(1)=XEND/2
           ERR=0.0D0
           DO I=1,1
             ERR = ERR + ((Y(I)-YEX(I))/YEX(I))**2
           END DO
           ERR=SQRT(ERR)
C --- PRINT FINAL SOLUTION
           write (6,77) ILIN,TOL,finish-start
 77        FORMAT(1X,'radau5, lin-alg = ',I1,'   TOL = ',E9.3,
     &           '    time = ',E8.2)
           WRITE (6,99) X,Y(1),ERR
 99        FORMAT(1X,'X =',F13.3,'    Y =',E20.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), 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), NR-1        
              XOUT=XOUT+1.0D0
              GOTO 10
           END IF
        END IF
 99     FORMAT(1X,'X =',F10.1,'    Y =',E18.10,'    NSTEP = ',I4)
        RETURN
        END
C
        SUBROUTINE FGAMI(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
        alpha = RPAR(1)
        kappa = RPAR(2)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)    
        PI=4.D0*ATAN(1.D0)
        APPINT=0.0D0
        DO I=3,LEXP+2
           IND=MEXP+(I-3)
           GAMMAI =  EXP(IND*HEXP) + kappa
           CI     =  EXP(alpha*IND*HEXP)
           F(I)   = - GAMMAI*Y(I) + Y(1)
           APPINT =  APPINT + CI*Y(I)
        END DO
        APPINT=APPINT*HEXP*(kappa**(1.D0-alpha))*SIN(alpha*PI)/PI
        F(1) = (1.D0 - Y(1))*ERF(SQRT(X)/2.D0)
     &          - SQRT(X/PI)*EXP(-X/4.D0) + Y(2) + 0.5D0
        F(2) = APPINT - Y(2)
        RETURN
        END
C
        SUBROUTINE JGAMI(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN with full linear algebra
        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
        alpha = RPAR(1)
        kappa = RPAR(2)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)    
        PI=4.D0*ATAN(1.D0)
        DFY = 0.0D0    
        DFY(1,1)=-ERF(SQRT(X)/2.D0)
        DFY(1,2)=1.0D0
        DFY(2,2)=-1.0D0
        FAC=HEXP*kappa**(1.D0-alpha)*SIN(alpha*PI)/PI
        II=3
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI =  EXP(IND*HEXP) + kappa
          CI =      EXP(alpha*IND*HEXP)
          DFY(2,II) = CI*FAC
          DFY(II,II) = - GAMMAI
          DFY(II,1) = 1.0D0
          II=II+1
        END DO
        RETURN
        END
C
        SUBROUTINE JGAMIEX(N,X,Y,DFY,LDFY,RPAR,IPAR)
C --- JACOBIAN for the option "sumexp" in the linear algebra
        IMPLICIT REAL*8 (A-H,O-Z)
        REAL*8 kappa
        REAL*8, dimension(20) :: RPAR
        INTEGER, dimension(20) :: IPAR       
        DIMENSION Y(N),DFY(LDFY,N)
        alpha = RPAR(1)
        kappa = RPAR(2)
        MEXP=IPAR(1)
        NEXP=IPAR(2)
        LEXP=NEXP-MEXP
        HEXP=RPAR(6)        
        PI=4.D0*ATAN(1.D0)
        DFY=0.0D0
        DFY(1,1)=-ERF(SQRT(X)/2.D0)
        DFY(1,2)=1.0D0
        DFY(2,2)=-1.0D0
        DFY(2,3)=1.0D0
        DFY(1,4)=1.0D0
        FAC=HEXP*kappa**(1.D0-alpha)*SIN(alpha*PI)/PI
        II=5
        DO I=1,LEXP
          IND=MEXP+(I-1)
          GAMMAI =  EXP(IND*HEXP) + kappa
          CI =      EXP(alpha*IND*HEXP)
          DFY(1,II) = CI*FAC
          DFY(2,II) = - GAMMAI
          II=II+1
        END DO
        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
         DO I=1,N
           AM(1,I)=1.0D0
         END DO
         AM(1,2) = 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
