***********************************************************************

      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