c Calcul des composantes du vecteur instantane de rotation dans le repere  
c terrestre et le repere local de l'observateur Oxyz :
c  - Oz : verticale du lieu, 
c  - Ox dans le plan horizontal, aligne selon le meridien local vers le SUD, 
c  - Oy axe horizontal directement perpendiculaire a Ox    
c
c Arguments en ligne d'entree : taper "rotation" 

	implicit none
	
	integer n, ndate,i,j,k,ii,jj,m,l
	parameter (n=20000)
	real*8 jjm1,x1,y1,UT1_UTC1,LOD1,dX1,dY1,tdj1,tdj2,
     &  tdjmax,tdjmin,lon,lat,h,tmjd0,tmjdd,mjd,
     &  sx1,sy1,sUT1_UTC1,sLOD1,sdX1,sdY1
        real*8 tdjdeb,tdjfin,tdj,dj,tdj_bis
	real*8 tmjd_deb,tmjd_fin
	real*8 jjm(n),x(n),y(n),ut1_utc(n),lod(n),dX(n),dY(n)
	real*8 jjm_bis(20),x_bis(20),y_bis(20),ut1_utc_bis(20),
     &         lod_bis(20),dX_bis(20),dY_bis(20)
	real*8 spline,tervon,interplin
	real*8 s,s1,s2,dt,omega,
     &         wx,wy,wz,wx1,wy1,wz1,cor_x,cor_y,cor_ut1,cor_lod
        real*8 tpi,degrad,secrad,masrad
	external dj 
	real*8 sintercal,ds,tmjd_utc,tartuf
	
	real*8 tps    , tumjd    , tu    , tmjd,
     &         tps_inf, tumjd_inf, tu_inf, tmjd_inf, tdj_inf,    
     &         tps_sup, tumjd_sup, tu_sup, tmjd_sup, tdj_sup,
     &         E, D, stheta, snro,XCIP, YCIP ,    
     &         xp, yp, dX0, dY0, xq(3,3), dxq(3,3),
     &         dxqdxp(3,3), dxqdyp(3,3) , dxqdut1(3,3),
     &         dxqdx(3,3) , dxqdy(3,3)  ,
     &         xq_inf(3,3), xq_sup(3,3) , xq_inv(3,3),xr(3,3),xpn(3,3)
	
	integer iclef,option,pm,ut,dP,hf
	integer an1,mois1,jour1,nh1,nm1,ns1,
     &          an2,mois2,jour2,nh2,nm2,ns2,
     &          an , mois,jour ,nh ,nm ,ns, ncom 
        character*100 
     &  clat,clon,ch,coption,cpm,cut,cdP,chf,
     &  cjour1,cmois1,can1,ch1,cm1,cs1,
     &  cjour2,cmois2,can2,ch2,cm2,cs2,
     &  fileinput,fileinputp,kom,kom1
	
	COMMON/conversion/tpi,secrad,degrad

        ! Conversion factor
        tpi= 6.28318530717958648d0
        secrad= tpi / (360.d0*3600.d0)
        degrad= tpi / 360.d0
		
	masrad=4.848d-9
	fileinput='/home/ftp/iers/eop/eopc04_05/eopc04_IAU2000.62-now' !  serie 9000 depuis 1962
        call getarg(1,can1)	    
	
	if (Len_trim(can1).eq.0) then
	write(6,*) '<<   Transf. de coordonnees du repere terrestre international Oxyz >>'
	write(6,*) '<<   au repere celeste international OXYZ et parametres derives    >>'
	write(6,*) 'Usage : rotation                     '
        write(6,*) '     [an1][mois1][jour1][h1][mn1][s1]'
        write(6,*) '     [an2][mois2][jour2][h2][mn2][s2]'
        write(6,*) '     [pas][option][pm][ut1][dP][hf][lat][lon]'
	write(6,*) ' [an1] ...... 1ere date : annee  '
        write(6,*) ' [mois1] .... 1ere date : mois   '
        write(6,*) ' [jour1]  ... 1ere date : jour   '       
        write(6,*) ' [h1]  ...... 1ere date : heure  '	  
        write(6,*) ' [mn1] ...... 1ere date : minute '	  
        write(6,*) ' [s1]  ...... 1ere date : seconde'	  
        write(6,*) ' [an2] ...... date fin  : annee  '
        write(6,*) ' [mois2] .... date fin  : mois   '
        write(6,*) ' [jour2]  ... date fin  : jour   '	 
        write(6,*) ' [h2]  ...... date fin  : heure  '    
        write(6,*) ' [mn2] ...... date fin  : minute '    
        write(6,*) ' [s2]  ...... date fin  : seconde'    
        write(6,*) ' [pas]  ..... pas temporel       '       
	write(6,*) ' [option] ... 1 : matrice de rotation'
	write(6,*) '          ... 2 : axe Oz projete sur le plan eq. J2000.0 OXY (")'
	write(6,*) '          ... 3 : vecteur inst. de rotation projete sur Oxy (nanorad/s)  '
	write(6,*) '          ... 4 : vecteur inst. de rotation projete sur OXY (nanorad/s)  '
	write(6,*) '          ... 5 : matrice de precession-nutation              '
	write(6,*) ' [pm] ....... mouvement du pole   :  oui :1 /non :0'
	write(6,*) ' [ut] ....... UT1 - UTC           :  oui :1 /non :0'
	write(6,*) ' [dP] ....... ecarts pole celeste :  oui :1 /non :0'
	write(6,*) ' [hf] ....... corr. de marees diurnes et semi-diurnes :  oui :1 /non :0'
	write(6,*) ' [lon] ...... longitude          '
	write(6,*) ' [lat] ...... latitude           '
       
	stop     
	endif

        call getarg(2,cmois1)
        call getarg(3,cjour1)   
        call getarg(4,ch1)
        call getarg(5,cm1)
        call getarg(6,cs1)        
        call getarg(7,can2)
        call getarg(8,cmois2)
        call getarg(9,cjour2)
        call getarg(10,ch2)
        call getarg(11,cm2)
        call getarg(12,cs2)
        call getarg(13,ch)
        call getarg(14,coption)
        call getarg(15,cpm)
        call getarg(16,cut)        
	call getarg(17,cdP)
	call getarg(18,chf)
	call getarg(19,clon)
	call getarg(20,clat)

	read(cjour1,*)  jour1
	read(cmois1,*)  mois1
	read(can1,*)    an1
	read(ch1,*)     nh1
	read(cm1,*)     nm1
	read(cs1,*)     s1 		
	read(cjour2,*)  jour2
	read(cmois2,*)  mois2
	read(can2,*)    an2
	read(ch2,*)     nh2
	read(cm2,*)     nm2
	read(cs2,*)     s2
	read(ch,*) h
	read(coption,*) option
	read(cpm,*) pm
	read(cut,*) ut
	read(cdP,*) dP
	read(chf,*) hf
	read(clat,*) lat
	read(clon,*) lon

	lon = lon * degrad
	lat = lat * degrad
	
c conversion en jour julien 
	
	ns1 = dint(s1)
	ns2 = dint(s2)
     
	tmjd_deb = DJ(an1,mois1,jour1,0,0,0.d0) - 2400000.5d0   ! date julienne modifiee UTC a 0hUTC 
	 
        tmjd0 = dint(tmjd_deb)  ! date mjd de reference (permet de minimiser les erreurs sur le calcul des temps par soustraction de cette date)      

	tmjd_deb = tmjd_deb - tmjd0 + dble(nh1)/24.D0 + dble(nm1)/60.D0/24.D0+ s1/86400.D0 
	tmjd_deb = tmjd_deb + tartuf(tmjd_deb + tmjd0)/86400.d0  ! DATE TAI

	tmjd_fin = DJ(an2,mois2,jour2,0,0,0.d0)- 2400000.5	
	tmjd_fin = tmjd_fin - tmjd0 + dble(nh2)/24.D0 + dble(nm2)/60.D0/24.D0+ s2/86400.D0 
	tmjd_fin = tmjd_fin + tartuf(tmjd_fin + tmjd0)/86400.d0  ! DATE TAI
	i = 0     		
	ncom = 13 + abs(an1 -1962) * 365  ! le debut de la serie 9000 est 1962
      
! -------------------------------------
! lecture et selection des EOP observes
! -------------------------------------

	open(10,file=fileinput, status='old')

	do j =1,ncom
	read(10,81) kom
	enddo	
81 	format(a100)
	
	do j=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.tmjd_deb-10.d0).and.(jjm1.le.tmjd_fin+10.d0)) then
		i = i + 1
		jjm(i) = jjm1
		x(i) = x1
		y(i) = y1
		ut1_utc(i) = ut1_utc1 - tartuf(jjm1+tmjd0) ! calcul de UT1-TAI pour faire  l'interpolation sur une echelle de temps
c		write(55,*) jjm1+tmjd0,ut1_utc(i)
		lod(i) = lod1
		dX(i) = dX1
		dY(i) = dY1
		endif
	
	enddo
		
10 	continue
	close(10)
	
	k = i  ! nombre de dates selectionnees		

 	! boucle temporelle - interpolation des EOP avec le pas h  
	ndate=dint((tmjd_fin-tmjd_deb)*86400.d0/h)+2
		
	do i=1,ndate	
	
		tmjd = tmjd_deb + dble(i-1) * h / 86400.d0   ! TAI		
		
		tmjd_utc = tmjd - tartuf(tmjd + tmjd0)/86400.d0 
		do j=1,100
		tmjd_utc = tmjd - tartuf(tmjd_utc + tmjd0)/86400.d0
		enddo
		
		if (k.ne.0) then ! EOP dans la periode selectionnee
99		format(f18.8)
		m = 0
	        	do l = 1,k	 
			if ( (jjm(l).ge.tmjd-10.0).and.(jjm(l).le.tmjd+10.0) ) then
			m = m+1	
			jjm_bis(m)     = jjm(l) 
			x_bis(m)       = x(l)
			y_bis(m)       = y(l) 
			ut1_utc_bis(m) = ut1_utc(l)
			lod_bis(m)     = lod(l) 
			dX_bis(m)      = dX(l) 
			dY_bis(m)      = dY(l) 
			endif	 	  
			enddo	
	
	     		if (m.lt.3) then 
			x1       = 0.d0
			y1       = 0.d0 
			ut1_utc1 = 0.d0
			lod1     = 0.d0 
			dX1      = 0.d0
			dY1      = 0.d0		
			else
			! --------------------  
			! interpolation a tmjd 
			! -------------------- 
	        	x1	= spline(jjm_bis,x_bis,m,tmjd)
	        	y1	= spline(jjm_bis,y_bis,m,tmjd)
	       	 	ut1_utc1= spline(jjm_bis,ut1_utc_bis,m,tmjd) 		
			ut1_utc1 = ut1_utc1  + tartuf(tmjd_utc + tmjd0)  ! on reforme UT1 - UTC	
c			write(55,1040) tmjd_utc + tmjd0,ut1_utc1,tartuf(tmjd_utc + tmjd0)
c1040			format(3f18.10)
			lod1	= spline(jjm_bis,lod_bis,m,tmjd)
	        	dX1	= spline(jjm_bis,dX_bis,m,tmjd)
	        	dY1	= spline(jjm_bis,dY_bis,m,tmjd)	
                			       
			endif
		
c 	       	 	x1	= interplin(jjm_bis(1),jjm_bis(m),x_bis(1),x_bis(m),tmjd)
c	        	y1	= interplin(jjm_bis(1),jjm_bis(m),y_bis(1),y_bis(m),tmjd)
c	        	ut1_utc1= interplin(jjm_bis(1),jjm_bis(m),ut1_utc_bis(1),ut1_utc_bis(m),tmjd)
c	        	lod1	= interplin(jjm_bis(1),jjm_bis(m),lod_bis(1),lod_bis(m),tmjd)
c	        	dX1	= interplin(jjm_bis(1),jjm_bis(m),dX_bis(1),dX_bis(m),tmjd)
c	        	dY1	= interplin(jjm_bis(1),jjm_bis(m),dY_bis(1),dY_bis(m),tmjd)
c			write(6,*) 'jjm_bis',jjm_bis(1),jjm_bis(m),ut1_utc_bis(1),ut1_utc_bis(m),ut1_utc1,'tdj',tdj
        	else ! aucun EOP dans la periode selectionnee
	 	x1       = 0.d0
		y1       = 0.d0 
		ut1_utc1 = 0.d0
		lod1     = 0.d0 
		dX1      = 0.d0
		dY1      = 0.d0	
         	endif 						      
       
		if (hf.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
			lod1	= lod1	   + cor_lod 
		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 (dP.eq.0) then
			dX1      = 0.d0
			dY1      = 0.d0
		endif 
		
c vitesse de rotation axiale 
		omega = 7.2921151467064d-5*(1.d0-lod1/86400.d0)
	      
	        tmjd     = tmjd + tmjd0 
		tmjd_utc = tmjd_utc + tmjd0 
		
		call MJDDATE(tmjd,jour,mois,an,nh,nm,s) ! date civile UTC a partir de tmjd TAI
								   
		tumjd = dint(tmjd) - 51544.5d0 ! partie entiere de la date UT1 / J2000 12h TAI
							
		tdj = tmjd + 2400000.5d0        ! date julienne TAI  
				  
		tdj2 = tdj - 2451545.d0          ! date TAI / J2000.0 en jour julien  
		tdj2 = tdj2 + 32.184d0/86400.d0  ! date TT  / J2000.0 en jour julien	
		tps  = tdj2 / 36525.d0           ! date TT  / J2000.0 en sc julien.           

	        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(56,1030) s,tmjd,tps,tumjd+51544.5d0,tu*86400.d0,tmjd_utc	                       
c1030       format(6f18.10)
	   iclef=1
	
c ###################################
c matrice d'orientation et sa derivee
c ###################################

		if (option.le.4) then 
	
		call nro_transf(tps, tumjd, tu, xp, yp, dX0, dY0, iclef, xq,
     .                  dxqdxp, dxqdyp, dxqdut1, dxqdx, dxqdy) 	
	
			if (option.eq.1) then
	
c matrice d'orientation
c      write(6,"(f22.13,i6,i3,i3,i3,i3,f18.14,/,3f22.12,/,3f22.12,/,3f22.12)")
c     &                tmjd,an, mois, jour, nh, nm, s, 
c     &  		     xq(1,1),xq(1,2),xq(1,3),
c     &  		     xq(2,1),xq(2,2),xq(2,3),
c     &  		     xq(3,1),xq(3,2),xq(3,3)	
		          
												
		       write(6,"(f22.13,i6,i3,i3,i3,i3,f18.12,3f22.12,3f22.12,3f22.12)")
     &                       tmjd,an, mois, jour, nh, nm,s, 
     &  		     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)		

	        	endif

			if (option.eq.2) then
			x1 = xq(1,3) / masrad / 1000.d0  ! conversion en seconde
			y1 = xq(2,3) / masrad / 1000.d0	 ! conversion en seconde
			write(6,"(f15.6,i6,i3,i3,i3,i3,f7.2,2f15.7)") tmjd,an, mois, jour, nh, nm, s, x1,y1 
			endif
	
			if ((option.eq.3).or.(option.eq.4)) then

				dt = 10.d0  ! pas de la derivation temporelle en seconde 

	        		tu_inf = tu - dt / 86400.d0                             

				tdj_inf  = tdj - dt / 86400.d0   
	        		tdj2     = tdj_inf - 2451545.d0 + tartuf(tmjd) / 86400.d0    	
				tdj2     = tdj2 + 32.184d0 / 86400.d0                         	
				tps_inf  = tdj2 / 36525.d0                 

	        		tu_sup = tu + dt/86400.d0            

				tdj_sup  = tdj + dt / 86400.d0      	
				tdj2     = tdj_sup - 2451545.d0 + tartuf(tmjd) / 86400.d0    	
				tdj2     = tdj2 + 32.184d0 / 86400.d0                           	
				tps_sup  = tdj2 / 36525.d0                  

				call nro_transf(tps_inf, tumjd, tu_inf, xp, yp, dX0, dY0, iclef, xq_inf,dxqdxp, dxqdyp, dxqdut1, dxqdx, dxqdy) 	
				call nro_transf(tps_sup, tumjd, tu_sup, xp, yp, dX0, dY0, iclef, xq_sup,dxqdxp, dxqdyp, dxqdut1, dxqdx, dxqdy) 	
	
	
		        	do ii=1,3                           
				do jj=1,3
				dxq(ii,jj)=(xq_sup(ii,jj)-xq_inf(ii,jj))/(2.d0*dt)
				enddo
				enddo

c Soit la matrice de transformation repere terrestre -> repere celeste :
c [X]= Q[x] 
c La formule de base est :
c          | 0  -wz  wy |
c	   | wz  0   -wx| = Q^{-1} d/dt(Q)  pour les coordonnees de w dans
c	   |-wy  wx   0 |                   le repere terrestre.
c
	
				call transmat(xq, xq_inv, 3, 3)
	
				if (option.eq.3) then ! vecteur instantane de rotation dans le repere terrestre 
c inversion de la matrice xq
	 	
	        		call mulmat(xq_inv , dxq, xr, 3, 3, 3)
   
c matrice d'orientation
c	write(6,"(3f22.12,/,3f22.12,/,3f22.12,/)") 
c     &                       xr(1,1)*1d9,xr(1,2)*1d9,xr(1,3)*1d9,
c     &                       xr(2,1)*1d9,xr(2,2)*1d9,xr(2,3)*1d9,
c     &                       xr(3,1)*1d9,xr(3,2)*1d9,xr(3,3)*1d9

				wx = xr(3,2)*1d9
				wy = xr(1,3)*1d9
				wz = omega*1d9

				wx1 = dsin(lat)*dcos(lon) * wx + dsin(lat)*dsin(lon) * wy-dcos(lat)*wz
				wy1 = -dsin(lon) * wx + dcos(lon) * wy 
				wz1 = dcos(lat)*dcos(lon) * wx + dcos(lat)*dsin(lon) * wy+dsin(lat)*wz
	
				x1 =  x1/3600.d0 * degrad * 7.292115d-5 * 1d9
				y1 = -y1/3600.d0 * degrad * 7.292115d-5 * 1d9
		
				write(6,"(f15.6,i6,i3,i3,i3,i3,f7.2,8f15.7)")  tmjd,an, mois, jour, nh, nm, s, wx,wy,wz,wx1,wy1,wz1,x1,y1 
			
				endif
	
	 
				if (option.eq.4) then ! vecteur instantane de rotation dans OXYZ
c inversion de la matrice xq
		
	        		call mulmat(dxq , xq_inv, xr, 3, 3, 3)
   
c matrice d'orientation
c	write(6,"(3f22.12,/,3f22.12,/,3f22.12,/)") 
c     &                       xr(1,1)*1d9,xr(1,2)*1d9,xr(1,3)*1d9,
c     &                       xr(2,1)*1d9,xr(2,2)*1d9,xr(2,3)*1d9,
c     &                       xr(3,1)*1d9,xr(3,2)*1d9,xr(3,3)*1d9

	wx = xr(3,2)*1d9
	wy = xr(1,3)*1d9
	wz = omega*1d9	
	
c composantes du vecteur instantane de rotation dans le repere terrestre 
				write(6,"(f15.6,i6,i3,i3,i3,i3,f7.2,3f15.7)")  tmjd,an, mois, jour, nh, nm, s,wx,wy,wz

				endif
	
			endif  ! fin de la condition : option 3 ou 4

		endif  ! fin de la condition : option <=4 

c ---------------------------------------	
c option=5 matrice de precession-nutation
c ---------------------------------------
		if (option.eq.5) then	
	        call precnut (dX1, dY1, iclef , tps, E, D, SNRO, stheta, XCIP, YCIP, xpn)	
c matrice de precession-nutation
	        write(6,"(f22.13,i6,i3,i3,i3,i3,f18.12,3f22.12,3f22.12,3f22.12)")  
     &                tmjd,an, mois, jour, nh, nm, s, 
     &  		     xpn(1,1),xpn(1,2),xpn(1,3),
     &  		     xpn(2,1),xpn(2,2),xpn(2,3),
     &  		     xpn(3,1),xpn(3,2),xpn(3,3)		
		endif

		
	enddo

	end

	include "trans.f" 

