* ULTIMA MODIFICACION: 6/7/12 

* FORTRAN code for solving second-order oscillatory differential systems:
* y''= f(t,y), y(T0) = y0,y'(T0) = y'0, in [T0,TF] by using the high-order 
* methods EFRKN8,EFRKN10 and EFRKN12


* The program uses next subroutines:
*
*    F: It evaluates F(T,Y)
*    RKN1E: Solve algebraic equations of the EFSV1 method
*    RKN2E: Solve algebraic equations of the EFSV2 method
*    COEFICIENTES: This routine use  h in the calculation of the coefficients
 
        PARAMETER (MAXDIM=200)
        INTEGER M,IER,NITER,NJAC,IFUN,NTITER,NFUN,N
        INTEGER INDICE,OPCIMET
        INTEGER OPCISUBDV,IS1,IS2,IS3,IS4,IS5,IS6
		LOGICAL FALLO
        REAL*4 secnds !an intrinsic timer function in G77
        REAL*4 TI
        DOUBLE PRECISION  H,TOL,PI,T,X,T0,TF,TMAX,MEG,NORMA,ADIF
		DOUBLE PRECISION  H1,H2,H3,H4,H5,H6
		DOUBLE PRECISION  T1,T2,T3,T4,T5,T6
		DOUBLE PRECISION OMEGA1,OMEGA2,OMEGA3,OMEGA4 OMEGA5,OMEGA6
        DOUBLE PRECISION OMEGAG1,OMEGAG2,OMEGAG3,OMEGAG4,OMEGAG5
        DOUBLE PRECISION OMEGAG6
        DOUBLE PRECISION Y0(MAXDIM),Y(MAXDIM),YF(MAXDIM)
		DOUBLE PRECISION YS1(MAXDIM),YS2(MAXDIM),YS3(MAXDIM)
		DOUBLE PRECISION YS4(MAXDIM),YS5(MAXDIM),YS6(MAXDIM)
		DOUBLE PRECISION YFS1(MAXDIM),YFS2(MAXDIM),YFS3(MAXDIM)
		DOUBLE PRECISION YFS4(MAXDIM),YFS5(MAXDIM),YFS6(MAXDIM)
		DOUBLE PRECISION DIF(MAXDIM),DIFP(MAXDIM)
        DOUBLE PRECISION Y0P(MAXDIM),YP(MAXDIM),YFP(MAXDIM) 
        DOUBLE PRECISION YPS1(MAXDIM),YPS2(MAXDIM),YPS3(MAXDIM) 
		DOUBLE PRECISION YPS4(MAXDIM),YPS5(MAXDIM),YPS6(MAXDIM)
        DOUBLE PRECISION YFPS1(MAXDIM),YFPS2(MAXDIM),YFPS3(MAXDIM) 
		DOUBLE PRECISION YFPS4(MAXDIM),YFPS5(MAXDIM),YFPS6(MAXDIM)
		DOUBLE PRECISION F0(MAXDIM)
        DOUBLE PRECISION MYF,MYFP,AYF,AYFP,NORMAY,NORMAYP 
        DOUBLE PRECISION TYFMAX, TYFPMAX  
        DOUBLE PRECISION LANDA,KA,W
        DOUBLE PRECISION DIGITS,DIGITSP,MEGP,NORMAP,ADIFP,TMAXP
        DOUBLE PRECISION EST,EREL,EABS,TOLPV,HINI,RN,RMAX,RMIN,P  
		DOUBLE PRECISION EPOS, EVEL 
        DOUBLE PRECISION A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2 
        DOUBLE PRECISION ANINF
        EXTERNAL ANINF  
        COMMON /OPCIONES/ INDICE,OPCIMET
        COMMON /VARIPRO/ LANDA,KA,W
        COMMON/DATOS/ A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2    
         

*used constants describe them below


*        parameters
*       ************
* M: The IVP dimension 
* T0, TF: starting  and end  points of the  integration  interval, respectively
* H:  step-size 
* Y0: initial values vector
* IER: error code 
* IC:  number of steps 
* TMAX: time point whith the maximum global error (MEG) 
* MEG: maximum global error 
* YF:  solution vector in TF 
* YFP: the derivate of the solution vector in TF 
* DIF: difference vector between exact  and numerical solution in TF
* NFUN: Number of derivative functions evaluated 
* T: starting point of the interval (t move to X = T + H) 
* X: retrieved point to take a step, X = T + H 
* IFUN: number of  F evaluations by step 
* NJAC: number of steps  
* Y: Initial vector Y = Y (T) at each step  
* LANDA:fitting frequency for each problem 



         PRINT *,'Enter the number of problem to solve:'
         PRINT *,' '
         PRINT *,' 1.-NONLINEAR OSCILLATORY, 2.-TWO BODIES'
         READ *, INDICE

* We Introduce  the interval of integration of the IVP
* and some parameters of the problem. 

        IF (INDICE .EQ. 1) THEN !Nonlinear oscillatory 
            M = 2         !M is the number of equations of the problem
            T0=0.D0
			TF=100.D0
            W=50.D0
            KA= 0.5D0 
  
        PRINT *,'The test problem  we integrate is:'
        PRINT *,'Y1"=-0.5*((w^2+k^2+1)*Y1 + (w^2-k^2-1)*Y2) + 0.5*k^2*(
     &Y1-Y2)^3'
        PRINT *,'Y2"=-0.5*((w^2-k^2-1)*Y1 + (w^2+k^2+1)*Y2) - 0.5*k^2*(
     &Y1-Y2)^3'
        PRINT *,' We move from T0 = 0. up TF=100.'

        PRINT *,' '
 
		ELSE IF (INDICE .EQ. 2) THEN !two bodies

            M = 2        !M is the number of equations of the problem
            T0=0.D0
			TF=20.D0
			W=1.D0

        PRINT *,'The test problem  we integrate is:'
        PRINT *,'Y1"=-Y1/sqrt((Y1^2+YI^2)^3)'
        PRINT *,'Y2"=-Y2/sqrt((Y1^2+YI^2)^3)'
        PRINT *,' We move from T0=0. up TF=20.'
        PRINT *,' '

 		END IF

* Select the method

        PRINT *,'Enter the method:'
        PRINT *,'(1-SV-I  Method  Storner/Verlet 1 stage)'
        PRINT *,'(2-SV-II  Method  Storner/Verlet 2 stages)'
        READ *, OPCIMET
* We calculate the number of steps to give

        NJAC = (TF-T0)/H + 0.001

* Read the index that it determines the number substeps and the combination 
* linear coefficient omega 

        PRINT *,'Enter the substages of progress:'

        PRINT *,'1-(1-2-3-4 subdiv.)/(1-2-3 subdiv.) 8(6)'
		PRINT *,'2-(1-2-3-4-5 subdiv.)/(1-2-3-4 subdiv.) 10(8)'
		PRINT *,'3-(1-2-3-4-5-6 subdiv.)/(1-2-3-4-5 subdiv.) 12(10)'
		READ *, OPCISUBDV

		OMEGA1 = 0.D0
        OMEGA2 = 0.D0
		OMEGA3 = 0.D0
		OMEGA4 = 0.D0
		OMEGA5 = 0.D0
		OMEGA6 = 0.D0

		OMEGAG1 = 0.D0
        OMEGAG2 = 0.D0
		OMEGAG3 = 0.D0
		OMEGAG4 = 0.D0
		OMEGAG5 = 0.D0
		OMEGAG6 = 0.D0

		IS1=1
		IS2=1
		IS3=1
		IS4=1
		IS5=1
		IS6=1

	    IF (OPCISUBDV .EQ. 1) THEN
		OMEGA1 = - 1.D0/360.D0
        OMEGA2 = 16.D0/45.D0
		OMEGA3 = -729.D0/280.D0
		OMEGA4 = 1024.D0/315.D0

		OMEGAG1 = 1.D0/24.D0
        OMEGAG2 = -16.D0/15.D0
		OMEGAG3 = 81.D0/40.D0
        OMEGAG4 = 0.D0 

		IS1=1
		IS2=2
		IS3=3
		IS4=4

		ELSE IF (OPCISUBDV .EQ.2) THEN
		OMEGA1 =  1.D0/8640.D0
        OMEGA2 =- 64.D0/945.D0
		OMEGA3 = 6561.D0/4480.D0
		OMEGA4 =- 16384.D0/2835.D0
		OMEGA5 = 390625.D0/72576.D0

		OMEGAG1 = - 1.D0/360.D0
        OMEGAG2 = 16.D0/45.D0
		OMEGAG3 = -729.D0/280.D0
		OMEGAG4 = 1024.D0/315.D0
		OMEGAG5 = 0.D0

		IS1=1
		IS2=2
		IS3=3
		IS4=4
		IS5=5
		ELSE IF (OPCISUBDV .EQ.3) THEN
		OMEGA1 = -1.D0/302400.D0
        OMEGA2 = 8.D0/945.D0
		OMEGA3 = -2187.D0/4480.D0
		OMEGA4 = 65536.D0/14175.D0
		OMEGA5 = -9765625.D0/798336.D0
		OMEGA6 = 17496.D0/1925.D0

		OMEGAG1 =  1.D0/8640.D0
        OMEGAG2 =- 64.D0/945.D0
		OMEGAG3 = 6561.D0/4480.D0
		OMEGAG4 =- 16384.D0/2835.D0
		OMEGAG5 = 390625.D0/72576.D0
	    OMEGAG6= 0.D0

		IS1=1
		IS2=2
		IS3=3
		IS4=4
		IS5=5
		IS6=6

		END IF

*The order of the methods is P, and it is used to calculate the HINI

        IF(OPCISUBDV.EQ.1)THEN 
        P=8
		ELSE IF (OPCISUBDV.EQ.2)THEN
		P=10
		ELSE IF (OPCISUBDV.EQ.3)THEN
		P=12
        END IF


* Absolute and relative errors

		EABS=0.D0
        PRINT*, 'ENTER THE RELATIVE ERROR'
        READ*, EREL

* We calculate the initial values

        CALL EXACTA(M,T0,Y0,Y0P)
        CALL F(M,T0,Y0,F0)
        TOLPV = MAX(1.D0,ANINF(Y0,M))*EREL+EABS


        IF(ANINF(Y0,M).LE. 1.D-10)THEN
         HINI=0.1
        ELSE
         HINI = (TOLPV/ANINF(F0,M))**(1.d0/DBLE(P)) 
        END IF

        H=MIN(HINI,(TF-T0)/10)
        WRITE(*,*) 'HINI ',H
 
        CALL COEFICIENTES(H)
		
        CALL EXACTA(M,T0,Y0,Y0P)

        DO 1 I=1,M

          Y(I)=Y0(I)
          YP(I)=Y0P(I)
		  YS1(I) = Y0(I)
		  YS2(I) = Y0(I)
		  YS3(I) = Y0(I)
		  YS4(I) = Y0(I)
		  YS5(I) = Y0(I)
		  YS6(I) = Y0(I)
		  YPS1(I) = Y0P(I)
 		  YPS2(I) = Y0P(I)
 		  YPS3(I) = Y0P(I)
 		  YPS4(I) = Y0P(I)
 		  YPS5(I) = Y0P(I)
 		  YPS6(I) = Y0P(I)
		
1       CONTINUE
 
        MEG = 0.D0
        MEGP= 0.D0
        NTITER = 0
        NFUN = 0
        IC=0
        NPF=0
        RMAX=3.D0
        RMIN=0.1D0


********************************************************
* We start the main loop                               *
********************************************************

        TI=SECNDS(0.0)
		T=T0

        DO WHILE(T.LT.TF)

			IF(IS1.GT. 0) H1=H/DBLE(IS1)		
	
            CALL COEFICIENTES(H1)
			
			DO 21 IS=0,IS1-1
				T1= T+ IS*H1
				IF (OPCIMET .EQ. 1) THEN
				CALL RKN1E (M,H1,T1,YS1,YPS1,YFS1,YFPS1)
				ELSE IF (OPCIMET .EQ. 2) THEN
			    CALL RKN2E (M,H1,T1,YS1,YPS1,YFS1,YFPS1)
				END IF
				DO 211 K=1,M
				  YS1(K)= YFS1(K)
				  YPS1(K) = YFPS1(K)
211				CONTINUE
				
21			CONTINUE
		
			
			IF(IS2.GT. 0) H2=H/DBLE(IS2)
            CALL COEFICIENTES(H2)

			DO 22 IS=0, IS2-1
				T2 = T + IS*H2
				IF(OPCIMET.EQ.1) THEN
					CALL RKN1E (M,H2,T2,YS2,YPS2,YFS2,YFPS2)
				ELSE IF (OPCIMET .EQ. 2) THEN				
					CALL RKN2E (M,H2,T2,YS2,YPS2,YFS2,YFPS2)
				END IF
				DO 221 K=1,M
				  YS2(K)= YFS2(K)
				  YPS2(K) = YFPS2(K)
221				CONTINUE
22			CONTINUE

				
			IF(IS3.GT. 0) H3=H/DBLE(IS3)
            CALL COEFICIENTES(H3)
			DO 23 IS=0, IS3-1
				T3 = T + IS*H3
				IF(OPCIMET.EQ.1) THEN
					CALL RKN1E (M,H3,T3,YS3,YPS3,YFS3,YFPS3)
				ELSE IF (OPCIMET .EQ. 2) THEN				
					CALL RKN2E (M,H3,T3,YS3,YPS3,YFS3,YFPS3)
				END IF
				DO 231 K=1,M
				  YS3(K)= YFS3(K)
				  YPS3(K) = YFPS3(K)
231				CONTINUE
23			CONTINUE
		

			IF(IS4.GT. 0) H4=H/DBLE(IS4)
            CALL COEFICIENTES(H4)
			DO 24 IS=0, IS4-1
				T4 = T + IS*H4
				IF(OPCIMET.EQ.1) THEN
					CALL RKN1E (M,H4,T4,YS4,YPS4,YFS4,YFPS4)
				ELSE IF (OPCIMET .EQ. 2) THEN				
					CALL RKN2E (M,H4,T4,YS4,YPS4,YFS4,YFPS4)
				END IF
				DO 241 K=1,M
				  YS4(K)= YFS4(K)
				  YPS4(K) = YFPS4(K)
241				CONTINUE
24			CONTINUE
		

		    IF(IS5.GT. 0) H5=H/DBLE(IS5)
            CALL COEFICIENTES(H5)

			DO 25 IS=0, IS5-1
				T5 = T + IS*H5
				IF(OPCIMET.EQ.1) THEN
					CALL RKN1E (M,H5,T5,YS5,YPS5,YFS5,YFPS5)
				ELSE IF (OPCIMET .EQ. 2) THEN				
					CALL RKN2E (M,H5,T5,YS5,YPS5,YFS5,YFPS5)
				END IF
				DO 251 K=1,M
				  YS5(K)= YFS5(K)
				  YPS5(K) = YFPS5(K)
251				CONTINUE
25			CONTINUE

		
			IF(IS6.GT. 0) H6=H/DBLE(IS6)
            CALL COEFICIENTES(H6)

			DO 26 IS=0, IS6-1
				T6 = T + IS*H6
				IF(OPCIMET.EQ.1) THEN
					CALL RKN1E (M,H6,T6,YS6,YPS6,YFS6,YFPS6)
				ELSE IF (OPCIMET .EQ. 2) THEN				
					CALL RKN2E (M,H6,T6,YS6,YPS6,YFS6,YFPS6)
				END IF
				DO 261 K=1,M
				  YS6(K)= YFS6(K)
				  YPS6(K) = YFPS6(K)
261				CONTINUE
26			CONTINUE

*			LINEAR COMBINATION OF SOLUTIONS

			DO 7 K=1,M
                YF(K) = OMEGA1*YS1(K)+OMEGA2*YS2(K)+OMEGA3*YS3(K)+
     &                  OMEGA4*YS4(K)+OMEGA5*YS5(K)+OMEGA6*YS6(K)
                YFP(K) = OMEGA1*YPS1(K)+OMEGA2*YPS2(K)+OMEGA3*YPS3(K)+
     &					OMEGA4*YPS4(K)+OMEGA5*YPS5(K)+OMEGA6*YPS6(K)
7       	CONTINUE

*         ESTIMATION OF THE ERROR IN POSITION AND SPEED

			EPOS=0.D0
			EVEL=0.D0

            DO 71 K=1,M

             EPOS= EPOS+ ((OMEGAG1-OMEGA1)*YS1(K) + 
     &             (OMEGAG2-OMEGA2)*YS2(K) + (OMEGAG3-OMEGA3)*YS3(K)+
     &             (OMEGAG4-OMEGA4)*YS4(K) + (OMEGAG5-OMEGA5)*YS5(K)+
     &             (OMEGAG6-OMEGA6)*YS6(K))**2

             EVEL= EVEL+ ((OMEGAG1-OMEGA1)*YPS1(K) + 
     &             (OMEGAG2-OMEGA2)*YPS2(K) + (OMEGAG3-OMEGA3)*YPS3(K)+
     &             (OMEGAG4-OMEGA4)*YPS4(K) + (OMEGAG5-OMEGA5)*YPS5(K)+
     &             (OMEGAG6-OMEGA6)*YPS6(K))**2


71          CONTINUE

           EPOS=DSQRT(EPOS)
           EVEL=DSQRT(EVEL)
           EST=MAX(EPOS,EVEL)
		   TOLPV=EABS+EREL*MAX(1.D0,ANINF(YF,M))

* CALCULATION OF FACTOR OF INCREASE OR REDUCTION STEP

          IF( EST.GE. 1.D-16) THEN
            RN= 0.9D0*(TOLPV/EST)**(1.D0/P)
          ELSE
             RN=RMAX
          END IF


* DECIDE if the step is acceptable, in THIS case we calculate the 
* new step and if not, reduce step. 


        IF ( EST.LT.TOLPV) THEN
           IC=IC+1 

		   IF(OPCISUB.EQ.1)THEN
		
			IS5=0
			IS6=0
		   ELSE IF(OPCISUB.EQ.2)THEN
		
			IS6=0
		   END IF

* number of function evaluations


			IF(OPCIMET.EQ.1)THEN
			NFUN=NFUN+ MAX(IS1,IS2,IS3,IS4,IS5,IS6)
			ELSE
			NFUN=NFUN+ MAX(IS1,IS2,IS3,IS4,IS5,IS6)+1 
			END IF

* We calculate the difference between  EXACT solution and the NUMERICAL solution 
* DIF = Y (X) - YF, as well as the maximum global error (MEG)

        X = T+H
		
        CALL EXACTA(M,X,Y,YP)

        DO 8 K=1,M
                DIF(K) = Y(K)-YF(K)
                DIFP(K) = YP(K)-YFP(K)
8       CONTINUE


        NORMA = 0.D0
        NORMAP= 0.D0
        NORMAY= 0.D0
        NORMAYP=0.D0

        DO 12 K=1,M
                ADIF = ABS(DIF(K))
                ADIFP=ABS(DIFP(K))
				AYF=ABS(YF(K))
				AYFP=ABS(YFP(K))
                IF (NORMA .LT. ADIF) THEN
                        NORMA = ADIF
                ENDIF
                IF (NORMAP.LT.ADIFP) THEN
                        NORMAP= ADIFP
                ENDIF
                IF (NORMAY .LT. AYF) THEN
                        NORMAY = AYF
                ENDIF
                IF (NORMAYP.LT.AYFP) THEN
                        NORMAYP= AYFP
                ENDIF

12      CONTINUE

		IF(DABS(X-TF).LT.(H/1.05D0)) THEN
		PRINT*,X,' ERRTF=',ADIF, '  ERRTFV=',ADIFP
		END IF

        IF (MEG .LT. NORMA) THEN
                        TMAX = X
                        MEG = NORMA
        ENDIF
        IF (MEGP .LT. NORMAP) THEN
                        TMAXP = X
                        MEGP = NORMAP
        ENDIF 

        IF (MYF .LT. NORMAY) THEN
                        TYFMAX = X
                        MYF = NORMAY
        ENDIF
        IF (MYFP .LT. NORMAYP) THEN
                        TYFPMAX = X
                        MYFP = NORMAYP
        ENDIF 


        DO 16 K=1,M
                Y(K) = YF(K)
                YP(K) = YFP(K)
16      CONTINUE


           IF (FALLO) THEN
               RMAX=1.D0
           END IF
           RN=MIN(RN,RMAX)
           IF (FALLO) THEN
               FALLO=.FALSE.
               RMAX=3.D0
           END IF
           RN=MAX(RMIN,RN)
           H=RN*H
           IF(X+H.GT.TF) THEN
              H=TF-X
           ELSE IF(X+2*H.GT.TF) THEN
              H=(TF-X)/2.D0
           END IF
           T=X
        ELSE !REJECTED STEP

           FALLO=.TRUE.
           NPF=NPF+1 !NUMBER of REJECTED STEP

           RN=MAX(RMIN,RN)
           H=RN*H   !REDUCE THE STEP

        END IF

        DO 17 I=1,M 
				YS1(I) = Y(I)
		  		YS2(I) = Y(I)
		  		YS3(I) = Y(I)
		  		YS4(I) = Y(I)
		  		YS5(I) = Y(I)
		  		YS6(I) = Y(I)
		  		YPS1(I) = YP(I)
 		  		YPS2(I) = YP(I)
 		  		YPS3(I) = YP(I)
 		  		YPS4(I) = YP(I)
 		  		YPS5(I) = YP(I)
 		  		YPS6(I) = YP(I)

17      CONTINUE
        CALL COEFICIENTES(H) 

		END DO

		TI=SECNDS(0.0)-TI
        DIGITS=DLOG10(MEG)
        DIGITSP=DLOG10(MEGP)
      
***********************************************************
* write the results 
***********************************************************

30      WRITE(*,*) '  '
	  WRITE(*,*) '  '

      WRITE (*,*) 'MEG = ',  MEG,' DIGITS = ',DIGITS
      WRITE(*,*) 'MEGP= ', MEGP, ' DIGITSP = ', DIGITSP
      WRITE (*,*) 'cpu = ',  TI
      WRITE(*,*) '  '
	  WRITE(*,*) '  '
      WRITE(*,*) 'IC= ',IC,' NFUN= ', NFUN,' NPF= ',NPF
      WRITE(*,*) '  '
	  WRITE(*,*) '  '

        IF(OPCISUBDV.EQ.1)THEN !usando SVI
        WRITE(*,*) 'NFUNTOTAL = ', (IC+NPF)*4, 'DIGITS = ', DIGITS
		ELSE IF (OPCISUBDV.EQ.2)THEN
        WRITE(*,*) 'NFUNTOTAL = ', (IC+NPF)*5, 'DIGITS = ', DIGITS
		ELSE IF (OPCISUBDV.EQ.3)THEN
        WRITE(*,*) 'NFUNTOTAL = ', (IC+NPF)*6, 'DIGITS = ', DIGITS
        END IF

	  WRITE(*,*) '  '
        DO 40 K=1,M
                WRITE (*,*) 'YF(',K,')=',YF(K),'  YFP(K)=',YFP(K)
40      CONTINUE

   		  WRITE(*,*) '  '
     
        PAUSE
        STOP
        END
*****************************************************************

        FUNCTION ANINF(X,M)
        DOUBLE PRECISION ANINF, X(M)
        ANINF = DABS(X(1))
        DO 1 I=2,M
1         ANINF = MAX( ANINF, DABS(X(I)))
        RETURN
        END


*****************************************************************
*  Subroutine for the exact solution computation
*****************************************************************


        SUBROUTINE EXACTA(M,X,Y,YP)
        PARAMETER(MAXDIM=200)
        INTEGER INDICE,OPCIMET,M,NITER
        DOUBLE PRECISION X,Y(MAXDIM),YP(MAXDIM)
        DOUBLE PRECISION LANDA,KA,W
		DOUBLE PRECISION TETAO,PI,XX,DKK,DKC2,SN,CN,DN
		DOUBLE PRECISION SN1,CN1,DN1,DKC21
		DOUBLE PRECISION TOL,IERR,EPS,U
        COMMON /OPCIONES/ INDICE,OPCIMET
        COMMON /VARIPRO/ LANDA,KA,W

        DO 10 I=1,M
          Y(I) =0.D0
10        YP(I) =0.D0


        IF(INDICE.EQ.1) THEN        !NONLINEAR OSCILLATORY 

		PI=4.D0*DATAN(1.D0)
		DKC21=1.D0-(KA*KA)
        CALL SNCNDN(X,DKC21,SN1,CN1,DN1)


        Y(1)=1.D0/DSQRT(2.D0)*(DCOS(PI/4.D0+W*X)-SN1)
		Y(2)=1.D0/DSQRT(2.D0)*(DCOS(PI/4.D0+W*X)+SN1)

        YP(1)=1.D0/DSQRT(2.D0)*(-W*DSIN(PI/4.D0+W*X)-
     &       DSQRT(1.D0-KA**2*SN1**2)*CN1)
        YP(2)=1.D0/DSQRT(2.D0)*(-W*DSIN(PI/4.D0+W*X)+
     &       DSQRT(1.D0-KA**2*SN1**2)*CN1)


        ELSE IF(INDICE.EQ.2) THEN    !2 BODIES
*********************************************************


        TOL =1.D-14
	    IERR=0
		NITER=0

        EPS=1.D-3 
		CALL FKEPLER(X,EPS,TOL,IERR,NITER,U)

		IF(IERR.NE.0) THEN
			PRINT*, 'ERROR IN THE RESOLUTION OF KEPLER'
			PRINT*, IERR
			GOTO 100
		END IF

100	    CONTINUE

         Y(1) = DCOS(U)-EPS
         Y(2) = DSQRT(1-EPS**2)*DSIN(U)
         YP(1) = -DSIN(U)/(1.D0-EPS*DCOS(U))
         YP(2) =DSQRT(1-EPS**2)*DCOS(U)/(1.D0-EPS*DCOS(U))

        END IF
        RETURN
        END

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

		SUBROUTINE FKEPLER(T,E,TOL,IERR,NITER,U)
*********************************************************
* Auxillary subroutine that solves by a Newton improved method
* the equation of kepler. It is included autoinicialization.

* arguments 
* t: point where the equation is solved 
* e: eccentricity or epsilon appearing in the equation 
* tol: tolerance 
* ierr: INTEGER code, if ierr is not 0, there is an error in resolution 
* niter: number of iterations that have been used for 
* reach the solution (max. 100) 
* U: value of the solution of the  kepler equation 
*    t= U(t)-e*sin(U(t))

		IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        PARAMETER(NITERMAX=100)
		INTEGER NITER
		DOUBLE PRECISION TOL,IERR,EPS,U

C WE CALCULATE A VALUE U0 SUITABLE FOR NEWTON
		DOSPI = 8.D0*DATAN(1.D0)        
        TSTART=DMOD(T,DOSPI)
        SIG = DSIGN(1.D0,DSIN(TSTART))
        U0= TSTART +SIG*0.85D0*E

c newton improved implementation

        ESIN= E*DSIN(U0)
        F0 = U0-ESIN-TSTART
        F=F0
        DO WHILE((DABS(F) .GE. TOL) .AND. (NITER.LE. NITERMAX))
           NITER=NITER+1
           ECOS= E*DCOS(U0)
           FP= 1.D0-ECOS
           FPP= ESIN
           FPPP=ECOS
           DELTAU0= -F/FP
           DELTAU0= -F/(FP+DELTAU0*FPP/2.D0)
           DELTAU0= -F/(FP+DELTAU0*FPP/2.D0 +(DELTAU0**2.D0)*FPPP/6.D0)
		   U0= U0+DELTAU0
           ESIN= E*DSIN(U0)
           F = U0-ESIN-TSTART
		   
        END DO
        IF (NITER . GE. NITERMAX) THEN
          IERR=10        
        END IF
        U= U0  

        RETURN
        END


*****************************************************************
        SUBROUTINE COEFICIENTES(H)
*****************************************************************
        INTEGER OPCIMET,INDICE
        DOUBLE PRECISION H, C0
        DOUBLE PRECISION  A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2 
        DOUBLE PRECISION P
        DOUBLE PRECISION S,T,U
        DOUBLE PRECISION LANDA,KA,W
		DOUBLE PRECISION XNU
        COMMON /OPCIONES/ INDICE,OPCIMET
        COMMON /VARIPRO/ LANDA,KA,W
        COMMON/DATOS/ A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2 

**** sv-1=1   SV-2=2  

        XNU= W*H
 
        IF (OPCIMET.EQ.1) THEN      ! SV-1
   
        PI = 4.D0*DATAN(1.D0)

        C0 = DCOS(XNU/2.0D0)
        GAM2 = 2.0D0*DSIN(XNU/2.D0)/(XNU*DCOS(XNU/2.D0))
		GAM1 = 2.D0*DSIN(XNU/2.D0)/XNU

        C1 =0.5D0
        C2 = 0.D0

        A21=0.D0
		

        D1 = 2.D0*DSIN(XNU/2.D0)/XNU  !YP
        D2 = 0.D0

        B1 = 0.5D0*GAM2*D1 !Y
        B2 = 0.D0

         
        ELSE IF (OPCIMET.EQ.2) THEN          !SV-2
   
        PI = 4.D0*DATAN(1.D0)

		GAM1 = DSIN(XNU)/XNU
        C1 = 0.D0
        C2 =1.D0 
		D1 = GAM1/(1.D0+DCOS(XNU))
		D2=D1

		B1 = GAM1*D1
        B2 = 0.0D0

		A21 = B1

        END IF
        RETURN
        END



*****************************************************************
        SUBROUTINE RKN1E (M,H,T,Y,YP,YF,YFP)
*********************************************************
* This subroutine solves algebraic equations of the EFRKN 
* one-stage explicit method (EFSV1)


* parameters:

* M: dimension of the IVP 
* H: step size (from T to T + H) 
* T: starting point of the interval 
* initial vector Y: (Y =Y (T)) 
* YF: solution  vector (YF = Y(T+H))
* YFP: derivative solution  vector  (YFP= YP(T+H)) 
* B1, D1, C1:  RKN1 coefficients
* PI: constant PI 
* T1: intermediate points used in RKN1 
* Y1: intermediate vectors used in RKN1 
* F1: vector corresponding to the evaluate  F in Y1 
* SUM: constant

* Declaration of variables

        PARAMETER(MAXDIM=200)
        INTEGER M,NITER,IER,IFUN
        INTEGER INDICE,OPCIMET
        INTEGER ORDER1(MAXDIM),ORDER2(MAXDIM)
        DOUBLE PRECISION H,T,TOL,PI, C0
        DOUBLE PRECISION XJAC1(MAXDIM,MAXDIM),XJAC2(MAXDIM,MAXDIM)
        DOUBLE PRECISION T1,T2,YNORM,YNORM1,YNORM2
        DOUBLE PRECISION Y(MAXDIM),YF(MAXDIM),Y1(MAXDIM),Y2(MAXDIM)
        DOUBLE PRECISION YP(MAXDIM),YFP(MAXDIM)
        DOUBLE PRECISION F1(MAXDIM),F2(MAXDIM),INCRE(MAXDIM)
        DOUBLE PRECISION G(MAXDIM)
        DOUBLE PRECISION XJAC(MAXDIM,MAXDIM)
        DOUBLE PRECISION LANDA,KA,W
        DOUBLE PRECISION  A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2 
        COMMON /OPCIONES/ INDICE,OPCIMET
        COMMON /VARIPRO/ LANDA,KA,W
        COMMON/DATOS/ A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2 


        PI = 4.D0*DATAN(1.D0)

		T1=T+ 0.5D0*H
        C0 = DCOS(W*H*0.5D0) 
         
        DO 2 K=1,M
                Y1(K) = C0*Y(K)+ H*C1*GAM1*YP(K)
2      CONTINUE

        CALL F(M,T1,Y1,F1)

            DO 24 K=1,M

             YFP(K) = YP(K) + H*D1*F1(K) 
			 YF(K)=Y(K)+H*GAM2*YP(K)+(H**2)*B1*F1(K)
 
24          CONTINUE


        RETURN
        END




*****************************************************************
        SUBROUTINE RKN2E (M,H,T,Y,YP,YF,YFP)
*****************************************************************

* This subroutine solves algebraic equations of the EFRKN 
* two-stage explicit method (EFSV2)


*  parameters 

* M: dimension of the IVP 
* H:step size (from T to T + H) 
* T: starting point of the interval 
* initial vector Y: (Y =Y (T)) 
* YF: solution  vector (YF = Y(T+H))
* YFP: solution derivative  vector (YFP= YP(T+H)) 
* A21,B1,B2,D1,D2,C1,C2: ERKN2 coefficients
* PI: constant PI 
* T1,T2:  intermediate points used in RKN2 
* Y1,Y2: intermediate vectors used in RKN2
* F1,F2: vectors corresponding to  evaluate  F in Y1  and Y2
* SUM: constant 

* Declaration of variables
*
        PARAMETER(MAXDIM=200)
        INTEGER M,NITER,IER,IFUN
        INTEGER INDICE,OPCIMET
        INTEGER ORDER1(MAXDIM),ORDER2(MAXDIM)
        DOUBLE PRECISION H,T,TOL,PI
        DOUBLE PRECISION XJAC1(MAXDIM,MAXDIM),XJAC2(MAXDIM,MAXDIM)
        DOUBLE PRECISION T1,T2,YNORM,YNORM1,YNORM2
        DOUBLE PRECISION Y(MAXDIM),YF(MAXDIM),Y1(MAXDIM),Y2(MAXDIM)
        DOUBLE PRECISION YP(MAXDIM),YFP(MAXDIM)
        DOUBLE PRECISION F1(MAXDIM),F2(MAXDIM),INCRE(MAXDIM)
        DOUBLE PRECISION G(MAXDIM)
        DOUBLE PRECISION XJAC(MAXDIM,MAXDIM)
        DOUBLE PRECISION LANDA,KA,W
        DOUBLE PRECISION  A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2 
        COMMON /OPCIONES/ INDICE,OPCIMET
        COMMON /VARIPRO/ LANDA,KA,W
        COMMON/DATOS/  A21,B1,B2,C1,C2,D1,D2,GAM1,GAM2 

        PI = 4.D0*DATAN(1.D0)

		T1 = T
        T2=  T+ H

        DO 2 K=1,M
                Y1(K) = Y(K)+ H*C1*YP(K)
2      CONTINUE

        CALL F(M,T1,Y1,F1)

        DO 3 K=1,M
                Y2(K) = Y(K)+ H*C2*GAM1*YP(K)+ (H**2)*A21*F1(K)
3      CONTINUE

        CALL F(M,T2,Y2,F2)

           DO 24 K=1,M

             YFP(K) = YP(K) + H*(D1*F1(K) + D2*F2(K))
			 YF(K)=Y(K)+H*GAM1*YP(K)+(H**2)*(B1*F1(K)+B2*F2(K))
 
24          CONTINUE


        RETURN
        END



*********************************************************************
        SUBROUTINE F(M,T,Y,YPP)
*********************************************************************


* This subroutine computes the second derivative function F(T,Y) where
*  T is a point  of the interval of integration

 

        PARAMETER(MAXDIM=200)

        INTEGER INDICE,OPCIMET,M
        DOUBLE PRECISION YPP(MAXDIM),Y(MAXDIM),AM(MAXDIM,MAXDIM)
        DOUBLE PRECISION LANDA,KA,W
		DOUBLE PRECISION TETAO,EPS,PI,T
        DOUBLE PRECISION XINC2,XINC4,COCI(MAXDIM),ALFA,EF(MAXDIM)
        COMMON /OPCIONES/ INDICE,OPCIMET
        COMMON /VARIPRO/ LANDA,KA,W

        IF(INDICE.EQ.1) THEN      !NOLINEAR OSCILLATORY

		YPP(1)=-(1.D0/2.D0)*((1.D0+KA**2+W**2)*Y(1)-(1.D0+KA**2-W**2)*
     &        Y(2))+(1.D0/2.D0)*KA**2*(Y(1)-Y(2))**3

		YPP(2)=-(1.D0/2.D0)*(-(1.D0+KA**2-W**2)*Y(1)+(1.D0+KA**2+W**2)*
     &        Y(2))-(1.D0/2.D0)*KA**2*(Y(1)-Y(2))**3
        ELSE IF(INDICE.EQ.2)  THEN     !2 bodies
       
		YPP(1)=-Y(1)/DSQRT((Y(1)**2+Y(2)**2)**3)
        YPP(2)=-Y(2)/DSQRT((Y(1)**2+Y(2)**2)**3)

        END IF
        RETURN
        END
******************************************************************
*                                                                                                                                                                                                                                                          
* SUBROUTINE that CALCULATES the SINE, COSINE and DN
* ELLIPTIC. PROBLEM OF DUFFING.
*                                                                                                                                                                                                                                                          
******************************************************************
        SUBROUTINE SNCNDN(UU,EMMC,SN,CN,DN)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        PARAMETER(CA=.0003D0)
        LOGICAL BO
        DIMENSION EM(13),EN(13)
        EMC=EMMC
        U=UU
        IF(EMC.NE. 0.D0) THEN
          BO=(EMC.LT.0.D0)
          IF(BO) THEN
              D=1.D0-EMC
              EMC=-EMC/D
              D=DSQRT(D)
              U=D*U
          ENDIF
          A=1.D0
          DN=1.D0
          DO 11 I=1,13
              L=I
              EM(I)=A
              EMC=DSQRT(EMC)
              EN(I)=EMC
              C=0.5D0*(A+EMC)
              IF (DABS(A-EMC).LE. CA*A) GOTO 1
              EMC=A*EMC
              A=C
11        CONTINUE
1         U=C*U
          SN=DSIN(U)
          CN=DCOS(U)
          IF(dabs(SN).le.1.D-16) GOTO 2
          A=CN/SN
          C=A*C
          DO 12 II=L,1,-1
              B=EM(II)
              A=C*A
              C=DN*C
              DN=(EN(II)+A)/(B+A)
              A=C/B
12        CONTINUE
          A=1.D0/DSQRT(C**2+1.D0)
          IF( SN.LT.0.D0) THEN
              SN= -A
          ELSE
              SN=A
          END IF
          CN=C*SN
2         IF(BO) THEN
              A=DN
              DN=CN
              CN=A
              SN=SN/D
          ENDIF
        ELSE
          CN=1.D0/DCOSH(U)
          DN=CN
          SN=DTANH(U)
        PRINT *,'SN',SN
        ENDIF
        RETURN
        END
*****************************************************************
******************************************************************



