***********************************************************************
      SUBROUTINE GAUSS_15(X,TS,TF,STEP,ERR,N,NOR,NI,NS,NBS,NF,FUN)
      ! Gauss-Everhart Integrator (Radau & Lobatto Spacings)
      ! Written by Avdyushev V. (2006)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (IPS0=100,EPS=1.D-15,SGM=3.16D0)
      REAL*8, ALLOCATABLE  :: A(:,:),B(:,:),C(:,:),D(:,:),E(:,:)
      REAL*8, ALLOCATABLE  :: SMS(:,:),SPC(:),BL(:,:)
      REAL*8, DIMENSION(N) :: X,Y,F,F0,P,XK
      LOGICAL LS,VS,INC,CNV
      COMMON/SPACING/SPACING(56)
      EXTERNAL FUN

      IF(TS.EQ.TF) RETURN
      IF(NOR.LT.2.OR.NOR.GT.15) 
     *STOP 'Inadmissible Order!'

      K=NOR/2 
      ALLOCATE(SPC(0:K),SMS(K,K))
      ALLOCATE(A(K,N),B(K,N),BL(K,N))
      ALLOCATE(C(0:K,0:K),D(0:K,0:K),E(0:K,0:K))
      SPC(0)=0.D0; L=K*(NOR-K-1); SPC(1:K)=SPACING(L+1:L+K)

      DO I=1,K; DO J=1,I-1
      SMS(I,J)=1.D0/(SPC(I)-SPC(J))
      END DO; END DO

      C=0.D0; C(0,0)=1.D0
      D=0.D0; D(0,0)=1.D0
      DO J=1,K; DO I=J,K
      C(I,J)=C(I-1,J-1)-SPC(I-1)*C(I-1,J)
      D(I,J)=D(I-1,J-1)+SPC(J  )*D(I-1,J)
      END DO; END DO

      E=0.D0; E(0:K,0)=1.D0
      DO I=1,K; DO J=1,I
      E(I,J)=E(I-1,J-1)+E(I-1,J)
      END DO; END DO

      DO I=1,K
      FLT=DFLOAT(I+1)
      C(:,I)=C(:,I)/FLT
      D(I,:)=D(I,:)*FLT
      E(I,:)=E(I,:)*FLT
      END DO

      DEG=1.D0/DFLOAT(K+1)
      RMX=SGM**DEG; RMN=1.D0/RMX
      DIR=(TF-TS)/DABS(TF-TS); H_=DABS(STEP)
      VS=ERR.NE.0.D0; CNV=NI.LE.0; INC=.TRUE. 
      T=TS; NS=0; NBS=0; NST=0; NF=0; IPS=IPS0
    3 H=H_*DIR; R=1.D0; LS=.FALSE.
      B=0.D0; XK=0.D0

      IF(H_.EQ.0.D0) THEN
      IF(.NOT.VS)
     *STOP 'Zero Constant Step!'
      CALL FUN(TS,X,F0); NF=NF+1
      H_=EPS; B_=0.D0
      DO WHILE(B_.EQ.0.D0)
      H_=H_*10.D0; H=H_*DIR; Y=X+H*F0
      IF(DABS(TF-TS).LE.H_) GOTO 1
      CALL FUN(TS+H,Y,F); NF=NF+1
      DO J=1,N; B_=B_+(F(J)-F0(J))**2; END DO 
      END DO
      H_=DSQRT(2.D0*ERR*H_/DSQRT(B_))
      H=H_*DIR
      END IF

    1 IF(DABS(TF-T).LE.H_) THEN
      STEP=H_; H=TF-T; H_=DABS(H) 
      R=R*H_/STEP; LS=.TRUE.; END IF

      IF(NS.LE.1) BL=B

      Q=1.D0
      DO L=1,K
      P=0.D0
      DO M=L,K
      P=P+E(M,L)*B(M,:)
      END DO 
      B(L,:)=B(L,:)-BL(L,:)
      Q=Q*R; BL(L,:)=P*Q/DFLOAT(L+1)
      B(L,:)=BL(L,:)+B(L,:)
      END DO

      DO L=1,K; P=0.D0; DO M=L,K
      P=P+D(M,L)*B(M,:)
      END DO; A(L,:)=P; END DO

      CALL FUN(T,X,F0); NF=NF+1

      DO IT=1,IPS

      DO I=1,K

      S=SPC(I); P=0.D0
      DO J=K,1,-1; P=S*(B(J,:)+P); END DO
      Y=X+H*S*(F0+P)

      CALL FUN(T+H*S,Y,F); NF=NF+1

      P=(F-F0)/S
      DO J=1,I-1; P=(P-A(J,:))*SMS(I,J); END DO
      DO J=1,I; B(J,:)=B(J,:)+C(I,J)*(P-A(I,:)); END DO
      A(I,:)=P

      END DO

      IF(ALL(DABS(Y-XK).LE.DABS(EPS*Y))) EXIT
      XK=Y

      END DO

      IF(IT.GT.IPS0) THEN
      IF(CNV.AND.NS.NE.0) NBS=NBS+1
      END IF

      IF(VS) THEN
      B_=0.D0 
      DO J=1,N
      B_=B_+B(K,J)**2
      END DO 
      IF(B_.NE.0.D0) THEN
      R=(ERR/DSQRT(B_)/H_)**DEG
      ELSE; R=RMX; END IF
      IF(NS.EQ.0) THEN 
      INC=.NOT.LS.AND.INC
      IF(INC.AND.R.GT.RMX.OR.R.LT.RMN) THEN
      NST=NST+1; IF(NST.GT.100.AND.R.GT.RMX) GOTO 4
      H_=H_*R; GOTO 3 
      END IF; END IF
    4 IF(R.GT.RMX) R=RMX
      END IF

      P=0.D0
      DO J=K,1,-1; P=B(J,:)+P; END DO
      X=X+H*(F0+P); T=T+H; NS=NS+1 
      H=H*R; H_=DABS(H); XK=0.D0

      IF(CNV) THEN; IPS=IPS0
      ELSE; IPS=NI; END IF

      IF(LS) GOTO 2; GOTO 1

    2 DEALLOCATE(A,B,C,D,E,BL,SMS,SPC)

      END SUBROUTINE GAUSS_15
***********************************************************************
      BLOCK DATA
      REAL*8, DIMENSION(56) :: SPACING
      COMMON/SPACING/SPACING
      DATA SPACING/
     *1.00000000000000000000000000000000,
     *0.66666666666666666666666666666667,
     *0.50000000000000000000000000000000,
     *1.00000000000000000000000000000000,
     *0.35505102572168219018027159252941,
     *0.84494897427831780981972840747059,
     *0.27639320225002103035908263312687,
     *0.72360679774997896964091736687313,
     *1.00000000000000000000000000000000,
     *0.21234053823915294397475811012400,
     *0.59053313555926528913507374793117,
     *0.91141204048729605260445385623054,
     *0.17267316464601142810085377187657,
     *0.50000000000000000000000000000000,
     *0.82732683535398857189914622812343,
     *1.00000000000000000000000000000000,
     *0.13975986434378055215208708112488,
     *0.41640956763108317994330233133708,
     *0.72315698636187617231995400231437,
     *0.94289580388548231780687880744588,
     *0.11747233803526765357449851302033,
     *0.35738424175967745184292450297956,
     *0.64261575824032254815707549702044,
     *0.88252766196473234642550148697967,
     *1.00000000000000000000000000000000,
     *0.09853508579882642612349889788775,
     *0.30453572664636390548538517627883,
     *0.56202518975261385599498747999477,
     *0.80198658212639182746420786320470,
     *0.96019014294853125765919330990667,
     *0.08488805186071653506398389301626,
     *0.26557560326464289309811405904562,
     *0.50000000000000000000000000000000,
     *0.73442439673535710690188594095438,
     *0.91511194813928346493601610698373,
     *1.00000000000000000000000000000000,
     *0.07305432868025888514812603418031,
     *0.23076613796994549908311663988435,
     *0.44132848122844986791860665819448,
     *0.66301530971884570090294702791922,
     *0.85192140033151570815002314750402,
     *0.97068357284021510802794972308684,
     *0.06412992574519669233127711938966,
     *0.20414990928342884892774463430102,
     *0.39535039104876056561567136982732,
     *0.60464960895123943438432863017268,
     *0.79585009071657115107225536569898,
     *0.93587007425480330766872288061033,
     *1.00000000000000000000000000000000,
     *0.05626256053692214646565219103231,
     *0.18024069173689236498757994280918,
     *0.35262471711316963737390777017124,
     *0.54715362633055538300144855765235,
     *0.73421017721541053152321060830661,
     *0.88532094683909576809035976293249,
     *0.97752061356128750189117450042915/
      END BLOCK DATA
***********************************************************************
      SUBROUTINE FUN(T,X,F)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(4),F(4)
      R=DSQRT(X(1)**2+X(2)**2)
      F(1:2)= X(3:4)
      F(3:4)=-X(1:2)/R**3
      END SUBROUTINE FUN
***********************************************************************
      SUBROUTINE LG_16(X,TS,TF,STEP,ERR,N,NOR,NI,NS,NBS,NF,FUN)
      ! Gauss-Everhart Integrator (Legendre Spacings)
      ! Written by Avdyushev V. (2006)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (IPS0=100,EPS=1.D-15,SGM=3.16D0)
      REAL*8, ALLOCATABLE  :: A(:,:),B(:,:),C(:,:),D(:,:),E(:,:)
      REAL*8, ALLOCATABLE  :: SMS(:,:),SPC(:),BL(:,:)
      REAL*8, DIMENSION(N) :: X,Y,F,F0,P,XK
      LOGICAL LS,VS,INC,CNV
      COMMON/SPACING/SPACING(36)
      EXTERNAL FUN

      IF(TS.EQ.TF) RETURN
      IF(NOR.LT.2.OR.NOR.GT.16) 
     *STOP 'Inadmissible Order!'
      IF(NOR/2*2-NOR.NE.0) 
     *STOP 'Odd Order!'

      K=NOR/2-1
      ALLOCATE(SPC(-1:K),SMS(0:K,0:K))
      ALLOCATE(A(0:K,N),B(0:K,N),BL(0:K,N))
      ALLOCATE(C(-1:K,-1:K),D(-1:K,-1:K),E(-1:K,-1:K))
      SPC(-1)=0.D0; L=K*(K+1)/2+1; SPC(0:K)=SPACING(L:L+K)

      DO I=0,K; DO J=0,I-1
      SMS(I,J)=1.D0/(SPC(I)-SPC(J))
      END DO; END DO

      C=0.D0; C(-1,-1)=1.D0
      D=0.D0; D(-1,-1)=1.D0
      DO J=0,K; DO I=J,K
      C(I,J)=C(I-1,J-1)-SPC(I-1)*C(I-1,J)
      D(I,J)=D(I-1,J-1)+SPC(J  )*D(I-1,J)
      END DO; END DO

      E=0.D0; E(0:K,0)=1.D0
      DO I=1,K; DO J=1,I
      E(I,J)=E(I-1,J-1)+E(I-1,J)
      END DO; END DO

      DO I=1,K
      FLT=DFLOAT(I+1)
      C(:,I)=C(:,I)/FLT
      D(I,:)=D(I,:)*FLT
      E(I,:)=E(I,:)*FLT
      END DO

      DEG=1.D0/DFLOAT(K+1)
      RMX=SGM**DEG; RMN=1.D0/RMX
      DIR=(TF-TS)/DABS(TF-TS); H_=DABS(STEP)
      VS=ERR.NE.0.D0; CNV=NI.LE.0; INC=.TRUE. 
      T=TS; NS=0; NBS=0; NST=0; NF=0; IPS=IPS0
    3 H=H_*DIR; R=1.D0; LS=.FALSE.
      B=0.D0; XK=0.D0

      IF(H_.EQ.0.D0) THEN
      IF(.NOT.VS)
     *STOP 'Zero Constant Step!'
      CALL FUN(TS,X,F0); NF=NF+1
      H_=EPS; B_=0.D0
      DO WHILE(B_.EQ.0.D0)
      H_=H_*10.D0; H=H_*DIR; Y=X+H*F0
      IF(DABS(TF-TS).LE.H_) GOTO 1
      CALL FUN(TS+H,Y,F); NF=NF+1
      DO J=1,N; B_=B_+(F(J)-F0(J))**2; END DO 
      END DO
      H_=DSQRT(2.D0*ERR*H_/DSQRT(B_))
      H=H_*DIR
      END IF

    1 IF(DABS(TF-T).LE.H_) THEN
      STEP=H_; H=TF-T; H_=DABS(H) 
      R=R*H_/STEP; LS=.TRUE.; END IF

      IF(NS.LE.1) BL=B

      Q=1.D0
      DO L=0,K
      P=0.D0
      DO M=L,K
      P=P+E(M,L)*B(M,:)
      END DO 
      B(L,:)=B(L,:)-BL(L,:)
      BL(L,:)=P*Q/DFLOAT(L+1)
      B(L,:)=BL(L,:)+B(L,:)
      Q=Q*R 
      END DO

      DO L=0,K; P=0.D0; DO M=L,K
      P=P+D(M,L)*B(M,:)
      END DO; A(L,:)=P; END DO

      DO IT=1,IPS

      DO I=0,K

      S=SPC(I); P=0.D0
      DO J=K,0,-1; P=S*(B(J,:)+P); END DO
      Y=X+H*P

      CALL FUN(T+H*S,Y,F); NF=NF+1

      P=F
      DO J=0,I-1; P=(P-A(J,:))*SMS(I,J); END DO
      DO J=0,I; B(J,:)=B(J,:)+C(I,J)*(P-A(I,:)); END DO
      A(I,:)=P

      END DO

      IF(ALL(DABS(Y-XK).LE.DABS(EPS*Y))) EXIT
      XK=Y

      END DO

      IF(IT.GT.IPS0) THEN
      IF(CNV.AND.NS.NE.0) NBS=NBS+1
      END IF

      IF(VS) THEN
      B_=0.D0 
      DO J=1,N
      B_=B_+B(K,J)**2
      END DO 
      IF(B_.NE.0.D0) THEN
      R=(ERR/DSQRT(B_)/H_)**DEG
      ELSE; R=RMX; END IF
      IF(NS.EQ.0) THEN 
      INC=.NOT.LS.AND.INC
      IF(INC.AND.R.GT.RMX.OR.R.LT.RMN) THEN
      NST=NST+1; IF(NST.GT.100.AND.R.GT.RMX) GOTO 4
      H_=H_*R; GOTO 3 
      END IF; END IF
    4 IF(R.GT.RMX) R=RMX
      END IF

      P=0.D0
      DO J=K,0,-1; P=B(J,:)+P; END DO
      X=X+H*P; T=T+H; NS=NS+1 
      H=H*R; H_=DABS(H); XK=0.D0

      IF(CNV) THEN; IPS=IPS0
      ELSE; IPS=NI; END IF

      IF(LS) GOTO 2; GOTO 1

    2 DEALLOCATE(A,B,C,D,E,BL,SMS,SPC)

      END SUBROUTINE LG_16
***********************************************************************
      BLOCK DATA
      REAL*8, DIMENSION(36) :: SPACING
      COMMON/SPACING/SPACING
      DATA SPACING/
     /0.50000000000000000000000000000000,
     /0.21132486540518711774542560974902,
     /0.78867513459481288225457439025098,
     /0.11270166537925831148207346002176,
     /0.50000000000000000000000000000000,
     /0.88729833462074168851792653997824,
     /0.06943184420297371238802675555359,
     /0.33000947820757186759866712044838,
     /0.66999052179242813240133287955162,
     /0.93056815579702628761197324444640,
     /0.04691007703066800360118656085030,
     /0.23076534494715845448184278964990,
     /0.50000000000000000000000000000000,
     /0.76923465505284154551815721035010,
     /0.95308992296933199639881343914970,
     /0.03376524289842398609384922275300,
     /0.16939530676686774316930020249005,
     /0.38069040695840154568474913915964,
     /0.61930959304159845431525086084036,
     /0.83060469323313225683069979750995,
     /0.96623475710157601390615077724700,
     /0.02544604382862073773690515797607,
     /0.12923440720030278006806761335960,
     /0.29707742431130141654669679396151,
     /0.50000000000000000000000000000000,
     /0.70292257568869858345330320603848,
     /0.87076559279969721993193238664039,
     /0.97455395617137926226309484202393,
     /0.01985507175123188415821956571526,
     /0.10166676129318663020422303176208,
     /0.23723379504183550709113047540538,
     /0.40828267875217509753026192881991,
     /0.59171732124782490246973807118009,
     /0.76276620495816449290886952459462,
     /0.89833323870681336979577696823792,
     /0.98014492824876811584178043428474/
      END BLOCK DATA
