	SUBROUTINE TETHYS(DAY,I,R)

C	SOUS PROGRAMME DE CALCUL DES COORDONNEES DIFFERENTIELLES 
C	RECTILIGNES EQUATORIALES B1950.0 DE MIMAS PAR RAPPORT A
C	LA PLANETE SATURNE
C	EN ENTREE :
C	*    DAY= DATE JULIENNE OBSERVATION EN JOUR FRACTIONNAIRE
C	*      I= INDICE DE LA THEORIE DU MOUVEMENT DU SATELLITE
C	ADOPTEE POUR CHAQUE ELEMENT ORBITAL ET PETITE VARIATION
C	        =1 G.STRUVE(1930)
C	        =2 Y.KOZAI (1957)
C	        =3 M.RAPAPORT(1977)PLANS DATE
C	        =4 I.G.CHUGUNOV (1984)
C	        =5 H.A.GARCIA (1972)
C	        =6 A.T.SINCLAIR (1977)
C	        =7 G.DOURNEAU (1987)
C	EN SORTIE :
C	*   R (K) = COORDONNEES RECTILIGNES EQUATORIALES B1950.0 
C	R(K=1,3)=DX,XY,DZ PAR RAPPORT A LA PLANETE SATURNE EN U.A. 

C	TRANSMISSIONS PAR COMMON A ETIQUETTE
C	*  COMMON/BLKTEVAL :

	IMPLICIT REAL*8(A-H,O-Z)
	CHARACTER*8 NOMSAT(8),NOM(8,24)
	DIMENSION EL(8,24),E0(8),AN(8),TETA0(8),GPI0(8),
	1 EX(8),AI(8),AA(8),TETAV(8),GPIV(8),A11(8),
	1 B1C(8),B1V(8),A13(8),A15(8),R(3),
	1 OMEG50C(8),ANGI50C(8),OMEG50V(8),ANGI50V(8),
	1 XTM(8),ACL(8),ACT(8)

	COMMON/BLKPI/PI,DPI,DEGRAD,RADEG
	COMMON/BLKDAT/DATE50,DATE2000,ANJBES
	COMMON/BLKTEVAL/E0,AN,TETA0,GPI0,EX,AI,AA,TETAV,GPIV,
	1 OMEG50C,ANGI50C,OMEG50V,ANGI50V,
	2 XTM,B1C,B1V,ACT,ACL

C	XTM Rapport des masses et demi-grands axes avec Mimas 
C	ACT Acceleration seculaire du noeud 
C	ACL Acceleration seculaire en longitude 

	COMMON/BLKMIVAL/EL1(8,24)
	EQUIVALENCE (EL(1,1),E0(1))

	DATA DS/9.538937D0/
	DATA TMK/216.11D0/

C	Calcul de la petite variation de A113 
C	correspondant a la variation de A111

	IF ((I-6)*(I-5).EQ.0) EL1(I,14)=-43.126D0
	B1V(I)=EL1(I,16)	!	ADOPTION DES VALEURS
	B1C(I)=EL1(I,15)	!	DES LIBRATIONS DE MIMAS.
	EL(8,16)=EL1(8,16)	!	FIXATION DES PETITES VARIATIONS
	EL(8,15)=EL1(8,15)	!	EGALES AUX VALEURS DE MIMAS
	  	
C	FIN DES DATA

	DT50=(DAY-DATE50)/36525.
	DJ=DAY-2411093.
	T=DT50*100.+1950.
	DT=DJ/365.25
	IT=1
1	OMEG50=(OMEG50C(I)+OMEG50V(I)*DT50)*DEGRAD
	ANGI50=(ANGI50C(I)+ANGI50V(I)*DT50)*DEGRAD
	CALL DANNO(DATE50,DATE50,OMEG50,ANGI50,
	1OMEG0,ANGI0,ANGN0,ANGJ0,POMEG0,DOM0)

C	FIXATION THEORIQUE DES ARGUMENTS DE LIBRATION (KOZAI,1957)

	IF ((I-4)*(I-5)*(I-6).NE.0)THEN
	COEFTM = - ( XTM(I) / 2 )
	A11(I) =  EL1(I,14)   * COEFTM
	A13(I) = (EL1(I,14)**3) * COEFTM
	1 / (3*(TMK*TMK) * (1-XTM(I))*(1-XTM(I)))
	A15(I) = (EL1(I,14)**5) * COEFTM
	1 / (5*(TMK**4)  * ((1-XTM(I))**4))
	ELSE IF ((I-5)*(I-6).EQ.0)THEN
	A11(I)=2.0616
	A13(I)=0.0335
	A15(I)=0.00094
	ELSE IF (I.EQ.4)THEN
	A11(I)=2.0072
	A13(I)=0.0341
	A15(I)=0.001
	END IF
	
	ARG1=B1V(I)*(T-B1C(I))
	ARG1=(DMOD(ARG1,360.D0))*DEGRAD

C	Acceleration seculaire en longitude ACL

	AL1=E0(I)+AN(I)*DJ + ACL(I)*(DT/100.)*(DT/100.)
	1+A11(I)*SIN(ARG1)+A13(I)*SIN(3.*ARG1)+A15(I)*SIN(5.*ARG1)
	AL1=(DMOD(AL1,360.D0))*DEGRAD
	GPI=GPI0(I)+GPIV(I)*DT
	GPI=(DMOD(GPI,360.D0))*DEGRAD
	ANOM=AL1-GPI
	GPI=GPI+DOM0

C	Acceleration seculaire du noeud ACT

	TETA=TETA0(I)+TETAV(I)*DT + ACT(I)*(DT/100.)*(DT/100.)
	TETA=(DMOD(TETA,360.D0))*DEGRAD+DOM0
	GPI=GPI-TETA
	AIR=AI(I)*DEGRAD
	TAA=AA(I)*DS*DEGRAD/3600.D0
	CALL EECR(TAA,EX(I),AIR,TETA,GPI,ANOM,R)
	CALL ROTA(R,0.D0,-ANGJ0,-ANGN0)
	RETURN
	END
