c Voila le fichier qui contient toute les subroutines necessaire
c pour calculer la matrice de transformation utilisant la NRO.
c Il faut fournir en entree le temps (3 variable differentes: tps, tumjd, tu)
c les coordonnes du pole terrestre xp, yp, les ecarts aux coordonnees du
c pole celeste dX0, dY0 (ex dpsi, deps fournie par la C04 / a UAI80,
c ils corresponde maintenant aux ecarts observes de X et Y par rapport
c a ceux calcule par le modele de Nicole qui est equivalent a IAU2000 que
c j'appelle X,Y UIA2000) plus la clef iclef qui n'a pas changee.  

c Calcul de la matrice d'orientation de la Terre a une epoque fixee
c en utilisant la serie combinee C04 : transformation de coordonnees 
c
c
c fonction matricielle orientation.f 
c [cjour cmois can cheure cmin csec cpm cut cdP ctides] 
c 

	implicit none
	
	integer n,i,j,tides
	parameter (n=41)
	real*8 jjm1,mjd,x1,y1,UT1_UTC1,LOD1,dX1,dY1,tdj,tmjd1,tmjd2	
	real*8 sx1,sy1,sUT1_UTC1,sLOD1,sdX1,sdY1
	real*8 jjm(n),x(n),y(n),ut1_utc(n),lod(n),dX(n),dY(n)
	real*8 spline,interplin,cor_x,cor_y,cor_ut1,cor_lod
	real*8 tai_utc,s,pm,ut,nut,tartuf
	
	real*8 tps, tumjd, tu,tmjd, tmjd_min, tmjd0, tmjd_utc,xp, yp, dX0, dY0, xq(3,3),
     .  dxqdxp(3,3), dxqdyp(3,3), dxqdut1(3,3),dxqdx(3,3),dxqdy(3,3)
	
	real*8 DJ
	external DJ
	integer iclef
	integer mois,jour,an,nh,nm
        character*80 can,cmois,cjour,cheure,cmin,csec,
     &  cpm,cut,cnut,ctides,fileinput,
     &  fileinputp,kom,kom1
		
        fileinput='/home/ftp/iers/eop/eopc04_05/eopc04_IAU2000.62-now' !  serie 9000 depuis 1962
	
        call getarg(1,can)
        
	if (Len_trim(can).eq.0) then
	write(6,*) '<<   Transf. de coordonnees du repere terrestre international Oxyz >>'
	write(6,*) '<<   au repere celeste international OXYZ                          >>'
	write(6,*) 'Usage : rotation                     '
        write(6,*) '     [an][mois][jour][h][mn][s][pm][ut][nut][tides]'
	write(6,*) ' [an] ...... date : annee  '
        write(6,*) ' [mois] .... date : mois   '
        write(6,*) ' [jour]  ... date : jour   '       
        write(6,*) ' [h]  ...... date : heure  '     
        write(6,*) ' [m] ......  date : minute '     
        write(6,*) ' [s]  ...... date : seconde'	  	
	write(6,*) ' [pm] ...... mouvement du pole   :  oui :1 /non :0'
	write(6,*) ' [ut] .......UT1 - UTC           :  oui :1 /non :0'
	write(6,*) ' [nut] ..... ecarts pole celeste :  oui :1 /non :0'
	write(6,*) ' [tides] ... corr. de marees diurnes et semi-diurnes :  oui :1 /non :0'
	stop
	endif
	
	call getarg(2,cmois)
        call getarg(3,cjour)
        call getarg(4,cheure)
        call getarg(5,cmin)
        call getarg(6,csec)
        call getarg(7,cpm)
        call getarg(8,cut)        
	call getarg(9,cnut)
	call getarg(10,ctides)

	read(can,*)    an	
	read(cmois,*)  mois
        read(cjour,*)  jour
	read(cheure,*) nh
	read(cmin,*)   nm
	read(csec,*)   s
	read(cpm,*) pm
	read(cut,*) ut
	read(cnut,*) nut
	read(ctides,*) tides
 	
c conversion en jour julien 
	
	
	tmjd = DJ(an,mois,jour,0,0,0.d0)- 2400000.5  ! date julienne modifiee UTC a 0hUTC 
	tmjd_utc = tmjd ! date MJD UTC
	tmjd0 = dint(tmjd)   ! date mjd de reference (permet de minimiser les erreurs sur le calcul des temps par soustraction de cette date)                           
	tmjd = tmjd - tmjd0
	tmjd = tmjd + dble(nh)/24.D0 + dble(nm)/60.D0/24.D0+ s/86400.D0 + tartuf(tmjd_utc)/86400.d0  ! DATE TAI
		
	tmjd_min = DJ(1962,1,1,0,0,0.d0) - 2400000.5
       
	j=0
	
	if (tmjd.ge.tmjd_min-tmjd0) then 
	 
c interpolation des EOP de la C04 a la date tdj

	tmjd1 = tmjd-10.d0
	tmjd2 = tmjd+10.d0
c	write(6,*) tdj1,tdj2	
	
c selection des EOP entre tmjd1 et tmjd2

	open(10,file=fileinput, status='old')
	
	do i =1,14
	read(10,81) kom
	enddo
	
81 	format(a120)

	j = 0
	do i=1,100000

	read(10,*,end=10) an,mois,jour,mjd,x1,y1,UT1_UTC1,LOD1,dX1,dY1,sx1,sy1,sUT1_UTC1,sLOD1,sdX1,sdY1
 
	jjm1 = mjd - tmjd0  ! date UTC
	
	if ((jjm1.ge.tmjd1).and.(jjm1.le.tmjd2)) then
	j = j+1
	jjm(j) = jjm1
	x(j) = x1
	y(j) = y1
	ut1_utc(j) = ut1_utc1  - tartuf(jjm1+tmjd0) ! calcul de UT1-TAI pour faire l'interpolation 	
	lod(j) = lod1
	dX(j) = dX1
	dY(j) = dY1
	endif

	enddo
	
10 	continue
	close(10)
	
	endif
		
c interpolation des EOP a la date tdj par splines cubiques

	if (j.ge.3) then

	x1       = spline(jjm,x,j,tmjd)
	y1       = spline(jjm,y,j,tmjd)
	ut1_utc1 = spline(jjm,ut1_utc,j,tmjd) + tartuf(tmjd_utc)  ! on repasse a UT1-UTC apres interpolation
	lod1     = spline(jjm,lod,j,tmjd)
	dX1      = spline(jjm,dX,j,tmjd)
	dY1      = spline(jjm,dY,j,tmjd)

c	x1       = interplin(jjm(1),jjm(j),x(1)      ,x(j),tmjd)
c	y1       = interplin(jjm(1),jjm(j),y(1)      ,y(j),tmjd)
c	ut1_utc1 = interplin(jjm(1),jjm(j),ut1_utc(1),ut1_utc(j),tmjd)
c	lod1     = interplin(jjm(1),jjm(j),lod(1),lod(j),tmjd)
c	dX1      = interplin(jjm(1),jjm(j),dX(1),dX(j),tmjd)
c	dY1      = interplin(jjm(1),jjm(j),dY(1),dY(j),tmjd)
 
	if (tides.eq.1) then ! corrections marees oceaniques
        call PMUT1_OCEANS (tmjd+tmjd0,cor_x,cor_y,cor_ut1,cor_lod)
	x1	= x1       + cor_x 	 
	y1	= y1	   + cor_y 
	ut1_utc1= ut1_utc1 + cor_ut1
	endif
	
c	write(6,*)  'spline',x1,y1,ut1_utc1,lod1,dX1,dY1
	else 
	x1       =  0.d0
	y1       =  0.d0 
	ut1_utc1 =  0.d0
	lod1     =  0.d0
	dX1      =  0.d0
	dY1      =  0.d0
c	write(6,*)  x1,y1,ut1_utc1,lod1,dX1,dY1
	endif

	if (pm.eq.0) then
	x1       = 0.d0
	y1       = 0.d0 
	endif
	if (ut.eq.0) then
	ut1_utc1 = 0.d0
	lod1     = 0.d0 
	endif
	if (nut.eq.0) then
	dX1      = 0.d0
	dY1      = 0.d0
	endif 
		 	
c mise en forme
        
	tmjd = tmjd + tmjd0
    	
	tumjd= dint(tmjd)  - 51544.5d0  ! partie entiere de la date UT1 / J2000 12h TAI
	
	tmjd1 = tmjd - 51544.5                     ! Jour julien - jour julien J2000.0 
	tmjd2 = tmjd1 + tartuf(tmjd_utc)/86400.d0  ! on forme tdj en TAI en ajoutant les sauts de secondes	
	tmjd1 = tmjd2 + 32.184d0/86400.d0         ! on forme tdj en TT	
	tps = tmjd1/36525.d0                      ! Date TT rapportee a J2000.0 
 		
           	
c UT1 - UT1 (J2000.0)  
          
	tu = (tmjd - 51544.5d0) - tumjd + (ut1_utc1 - tartuf(tmjd_utc))/ 86400.d0  ! fraction de jour UT1 depuis 0h TAI               
 	
	xp  = x1
	yp  = y1	

	dX0 = dX1
	dY0 = dY1
              
c	write(6,*) 'tps=',tps,tumjd,tu,x1,y1,dX1,dY1,ut1_utc1,nm,nh,s	                       
	                       
	iclef=1
c	write(6,*) ' tps=',tps, '    tumjd=',tumjd,'   tu=',tu, '  ut1_utc1 =',ut1_utc1 
	call nro_transf(tps, tumjd, tu, xp, yp, dX0, dY0, iclef, xq,
     .                      dxqdxp, dxqdyp, dxqdut1, dxqdx, dxqdy) 	

c        write(6,1040) tdj
c        write(6,*) tps
c        write(6,1042) tumjd
c	 write(6,*) tu
c	 write(6,*) 's,ns=',s,ns

c	kom= xq(1,1)\\xq(1,2)\\xq(1,3)\\xq(2,1)\\xq(2,2)\\xq(2,3)
	
c	do i=1,3
c	write(6,"(3f18.12)") (xq(i,j),j=1,3)
c	enddo

	write(6,"(9f22.12)")
     &                      xq(1,1),xq(1,2),xq(1,3),
     &                      xq(2,1),xq(2,2),xq(2,3),
     &                      xq(3,1),xq(3,2),xq(3,3)

	end

	include "trans.f" 
    	
