{ for the numerical integration
   method presented by professor Everhart E. Celestial Mech., 1974, v.10 }

 UNIT uneveras ;

 INTERFACE

 Uses
   unsagcot ; { for type TVect3 }

 Var
   StepRada  : Array[1..8] of Extended ;
   DimCoefR  : Array[1..8] of Extended ;
   DimCoefV  : Array[1..8] of Extended ;
   CoefEver  : Array[1..8] of TVect3 ;
   EverAlfa  : Array[1..8] of TVect3 ;
   MemAlpha  : Byte ;
   AlphaMom  : Array[1..2] of Array[1..4] of Extended ;

 Procedure EverInit ;

 Procedure EverCoefNullo ;

 Procedure ClcVectFromPolinom ( DStep  : Extended ;
                                ZeroR , ZeroV  : TVect3 ;
                                Var
                                  VectR , VectV  : TVect3 ) ;

 Procedure PosVelEverGrator ( StepFull         : Extended ;
                              Var
                                TBeg           : Extended ;
                              Var
                                ZeroR , ZeroV  : TVect3 ;
                              Var
                                TEnd           : Extended ;
                              Var
                                VectR , VectV  : TVect3 ) ;

 IMPLEMENTATION

 Uses
   unequras ; { for the procedure ClcRightHandCoor }

 Var
   DifT  : Array[1..8] of Extended ;
   DifD  : Array[1..8] of Extended ;
   Coef  : Array[1..8] of Array[1..8] of Extended ;

 procedure EverInit ; { initial value to some variables }
 begin
   StepRada[1]:=0.000000000000000000;
   StepRada[2]:=0.056262560526922147;
   StepRada[3]:=0.180240691736892365;
   StepRada[4]:=0.352624717113169637;
   StepRada[5]:=0.547153626330555383;
   StepRada[6]:=0.734210177215410532;
   StepRada[7]:=0.885320946839095768;
   StepRada[8]:=0.977520613561287501;
   DimCoefR[1]:=0.500000000000000000;
   DimCoefR[2]:=0.166666666666666667;
   DimCoefR[3]:=0.083333333333333333;
   DimCoefR[4]:=0.050000000000000000;
   DimCoefR[5]:=0.033333333333333333;
   DimCoefR[6]:=0.023809523809523810;
   DimCoefR[7]:=0.017857142857142857;
   DimCoefR[8]:=0.013888888888888889;
   DimCoefV[1]:=1.000000000000000000;
   DimCoefV[2]:=0.500000000000000000;
   DimCoefV[3]:=0.333333333333333333;
   DimCoefV[4]:=0.250000000000000000;
   DimCoefV[5]:=0.200000000000000000;
   DimCoefV[6]:=0.166666666666666667;
   DimCoefV[7]:=0.142857142857142857;
   DimCoefV[8]:=0.125000000000000000;
 end;

 Procedure EverCoefNullo ;  { initial moment nullo to start }
 Var
   i,j  : Byte ;
 Begin
   For i:=1 To 8 Do
     For j:=1 To 3 Do
       Begin
         CoefEver[i][j]:=0.0;
         EverAlfa[i][j]:=0.0;
       End;
   MemAlpha:=0;
   For i:=1 To 2 Do
     For j:=1 To 4 Do
       AlphaMom[i][j]:=0.0;
 End;

 { to obtain vector
      on the moment  t  from the coefficients of the polinoms
      constructed in integrator process presented by professor Everhart }


 Procedure ClcVectFromPolinom ( DStep  : Extended ;
                                ZeroR , ZeroV  : TVect3 ;
                                Var
                                  VectR , VectV  : TVect3 ) ;
 Var  { DStep in day from nullo point }
   i,j   : Byte ;
   CurR  : TVect3 ;
   CurV  : TVect3 ;
 Begin
   For j:=1 To 3 Do
     Begin
       CurV[j]:=DimCoefV[8]*CoefEver[8][j];
       CurR[j]:=DimCoefR[8]*CoefEver[8][j];
     End;
   For i:=7 DownTo 1 Do
     For j:=1 To 3 Do
       Begin
         CurV[j]:=DimCoefV[i]*CoefEver[i][j]+CurV[j]*DStep;
         CurR[j]:=DimCoefR[i]*CoefEver[i][j]+CurR[j]*DStep;
       End;
   For j:=1 To 3 Do
     Begin
       VectV[j]:=ZeroV[j]+CurV[j]*DStep;
       VectR[j]:=ZeroR[j]+(ZeroV[j]+CurR[j]*DStep)*DStep;
     End;
 End;

 Procedure FromPrevious ( Var  TBeg           : Extended ;
                          Var  ZeroR , ZeroV  : TVect3 ;
                          Var  TEnd           : Extended ;
                          Var  VectR , VectV  : TVect3 ) ;
 Var
   i  : Byte ;
 Begin
   TBeg:=TEnd; { to change condition from previous step }
   For i:=1 To 3 Do
     Begin
       ZeroR[i]:=VectR[i];
       ZeroV[i]:=VectV[i];
     End;
 End;

 Procedure CoeFromRecur ( StepFull  : Extended ) ;
 Var
   i,j  : Byte ;
 Begin
   DifT[1]:=0; { StepFull in day }
   For i:=2 To 8 Do
     Begin
       DifT[i]:=StepRada[i]*StepFull;
       DifD[i]:=DifT[i];
     End;
   For i:=1 To 8 Do
     For j:=1 To 8 Do  Coef[i][j]:=0;
   { the simple coefficients from recurrence formula }
   For i:=1 To 8 Do  Coef[i][i]:=1;
   For i:=2 To 8 Do
     Coef[i][1]:=-DifT[i]*Coef[i-1][1];
   For i:=3 To 8 Do
     For j:=2 To i-1 Do
       Coef[i][j]:=Coef[i-1][j-1]-DifT[i]*Coef[i-1][j];
 End;

 Procedure CoeFromAlpha ( icur  : Byte );
 Var
   ia,j,k  : Byte ;
 Begin
   For j:=1 To 3 Do
     For ia:=2 To icur Do   { to improve the value of coefficients }
       Begin
         CoefEver[ia][j]:=EverAlfa[ia][j];
         For k:=ia+1 To 8 Do
           CoefEver[ia][j]:=CoefEver[ia][j]
                           +Coef[k-1][ia-1]*EverAlfa[k][j];
       End;
 End;

 Procedure ToChangeAlpha ( TBeg  : Extended ) ;
 Var
   j  : Byte ;
   d  : Extended ;
 Begin
   If  MemAlpha < 2  Then  Exit ; { no action no alpha values in memory }
   d:=(TBeg-AlphaMom[2][4])/(AlphaMom[2][4]-AlphaMom[1][4]);
   For j:=1 To 3 Do { linear interpolation }
     EverAlfa[2][j]:=AlphaMom[2][j]+(AlphaMom[2][j]-AlphaMom[1][j])*d;
   CoeFromAlpha(2);  { to improve the value of the second coefficients }
 End;

 Procedure AlphaValueToMemory ( TBeg  : Extended ) ;
 Var
   j  : Byte ;
 Begin
   MemAlpha:=MemAlpha+1;
   If  MemAlpha = 3
     Then
       Begin
         MemAlpha:=2;
         For j:=1 To 3 Do
           AlphaMom[1][j]:=AlphaMom[2][j];
         AlphaMom[1][4]:=AlphaMom[2][4];
       End;
   For j:=1 To 3 Do
     AlphaMom[MemAlpha][j]:=EverAlfa[2][j];
   AlphaMom[MemAlpha][4]:=TBeg;
 End;

 { to integrate the full function
      by the use of the method presented by professor Everhart }

 Procedure PosVelEverGrator ( StepFull         : Extended ;
                              Var
                                TBeg           : Extended ;
                              Var
                                ZeroR , ZeroV  : TVect3 ;
                              Var
                                TEnd           : Extended ;
                              Var
                                VectR , VectV  : TVect3 ) ;
 Var
   i,j,ia,it  : Byte ;
   tc,dst,dt  : Extended ;
   FunB,FunC  : TVect3 ;
   RCur,VCur  : TVect3 ;
 Begin

   FromPrevious(TBeg,ZeroR,ZeroV,TEnd,VectR,VectV);{ to new condition }

   If  Abs(StepFull) < 1.0e-12  Then  Exit; { no step needs }

   CoeFromRecur(StepFull); { StepFull in day }

   { value of the forces at the first point  t0 }
   ClcRightHandCoor(TBeg,ZeroR,FunB); { TBeg in mod. julian day }

   { initial values of the coefficients of the polinoms for forces }
   For j:=1 To 3 Do
     CoefEver[1,j]:=FunB[j]; { the first coefficients }

   ToChangeAlpha(TBeg); { if there are two previous alpha values }

   it:=0;
 Repeat
   it:=it+1;
   For i:=2 To 8 Do
     Begin

       dst:=DifT[i];
       ClcVectFromPolinom(dst,ZeroR,ZeroV,RCur,VCur);

       tc:=TBeg+DifD[i]; { step in day TCur in mod. julian day }

       { value of the forces at the next point i }
       ClcRightHandCoor(tc,RCur,FunC); { from UnForces }

       dt:=1.0/DifT[i];
       For j:=1 To 3 Do
         EverAlfa[i][j]:=(FunC[j]-FunB[j])*dt;

       { to get the next alpha value }
       For ia:=3 To i Do
         Begin
           dt:=1.0/(DifT[i]-DifT[ia-1]);
           For j:=1 To 3 Do
             EverAlfa[i][j]:=(EverAlfa[i][j]-EverAlfa[ia-1][j])*dt;
         End;

       CoeFromAlpha(i);  { to improve the value of coefficients }

     End;
 Until it = 4 ;

   { coefficients of the polinoms have been calculated
     it is possible to obtain the vectors at the last point }
   TEnd:=TBeg+StepFull; { StepFull in day }
   ClcVectFromPolinom(StepFull,ZeroR,ZeroV,VectR,VectV);

   AlphaValueToMemory(TBeg) ;

 End;

 END.