! BEGIN ROCOTLIB 3.2                                                    
!                                                                       
! XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! XX                                                                  XX
! XX   Robert's Coordinates Transformation Library  ROCOTLIB          XX
! XX                                                                  XX
! XX            Revised  version 2.2, December 2011                   XX
! XX                    ___________________                           XX
! XX                                                                  XX
! XX             Version 1.0 initially supported by                   XX
! XX                                                                  XX
! XX                  EUROPEAN  SPACE  AGENCY                         XX
! XX                                                                  XX
! XX               Study of the Cluster Mission                       XX
! XX                 Planning Related Aspects                         XX
! XX          within the Numerical Simulations Network                XX
! XX                                                                  XX
! XX            Patrick ROBERT, CRPE, November 1992                   XX
! XX                    ___________________                           XX
! XX                                                                  XX
! XX Further versions developped by the Author:                       XX
! XX                                                                  XX
! XX    version 1.0, November 1992                                    xx
! XX    version 1.1, July     1993                                    xx
! XX    version 1.2, January  1995                                    xx
! XX    version 1.3, July     2000 (Julian day 2000 et sun dir)       xx
! XX    version 1.4, June     2001 (for automatic documentation)      xx
! XX    version 1.5, December 2001 (adding cp_sunrise_sunset)         xx
! XX    version 1.6, Juin     2002 (upgrade IGRF -> 2005)             xx
! XX    version 1.7, December 2002 (Version for CDPP)                 xx
! XX    version 1.8, November 2003 (adding t_sr2_to_sr)               xx
! XX    version 1.9, March    2004 (compatibility with IDL)           xx
! XX    version 2.0, November 2006 (Update  IGRF -> 2010)             xx
! XX    version 2.1, November 2010 (Update  IGRF -> 2015)             xx
! XX    version 2.2, December 2011 (cp_sunrise_sunset on polar zone   xx
! XX    version 3.0, February 2016 (IGRF -> table of geomagnetic pole XX
! XX                                -> 2020 + some other coord. sys.) XX
! XX    version 3.1, January  2019 (add name compatibility with V2.2) XX
! XX    version 3.2, February 2020 (add trans. with Euler angles)     XX
! XX                                                                  XX
! XX Copyright 1992, Patrick ROBERT, CNRS-ESA, All Rights reserved    XX
! XX                    ___________________                           XX
! XX                                                                  XX
! XX         For details, see the orginal document untitled:          XX
! XX CLUSTER Software Tools, part I: Coordinate Transformation LibraryXX
! XX             Document de travail DT/CRPE/1231                     XX
! XX             Patrick Robert, CRPE/TID, Juillet 1993               XX
! XX                                                                  XX
! XX         Available in CDPP:                                       XX
! XX             ROCOTLIB: a coordinate Transformation Library        XX
! XX             for Solar-Terrestrial studies                        XX
! XX             Patrick ROBERT, version 1.7 - January 2003,          XX
! XX             Rapport Interne no RI-CETP/02/2003                   XX
! XX                    ___________________                           XX
! XX                                                                  XX
! XX         The previous version 3.0 contains new transformations     X
! XX         (TPN, MVA), and the cp_sunrise_sunset module allowing    XX
! XX         computation of sunrise and sunset anywhere on Earth,     XX
! XX         including polar zones.                                   XX
! XX                                                                  XX
! XX         The determination of the dipole axis direction           XX
! XX         has been replaced by a table available from years        XX
! XX         1900 to 2020, in order to have an easy way to code       XX
! XX         maintenance in the coming years.                         XX
! XX                                                                  XX
! XX         Operations on matrix (somme, product, inversion,         XX
! XX         computation of eigen vectors, diagonalization...)        XX
! XX         required for Minimum Variance Analysis coordinates       XX
! XX         have been added with prefix "mat_".                      XX
! XX                                                                  XX
! XX         The previous and temporary V3.1 is the same as 3.0       XX
! XX         but is compatible with V2.2 version concerning           XX
! XX         subroutine name (ex: tmaggsm -> t_mag_to_gsm)            XX
! XX                                                                  XX
! XX         The present V3.2 version contains transformations        XX
! XX         with Euler angles and their interpolation.               XX
! XX                                                                  XX
! XX         Fortran 90, C, and IDL codes are also available.         XX
! XX                                                                  XX
! XX                                                                  XX
! XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
!                                                                       
!     ***************************************************************0**
!                                                                       
!     Reminder for the developper:                                      
!     Any changes leadind to a new version implies following actions:   
!        - Update the file header above                                 
!        - Update the subroutine g_rocot_version_number                 
!        - Update the subroutine print_rocot_info                       
!        - Update the readme.txt file in the full package               
!                                                                       
!     ***************************************************************0**
!                                                                       
!                                                                       
      subroutine cp_angle_and_ratio(ux,uy,uz,vx,vy,vz,angle,ratio) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_angle_and_ratio beetween U and V vectors         
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : ux,uy,uz                                                 
!              vx,vy,vz                                                 
!                                                                       
! *   Output : sp=U.V                                                   
!              angle=angle beetween U and V (radians)                   
!              ratio= mod(U)/mod(V)                                     
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      double precision u1,u2,u3,v1,v2,v3,dp,ru,rv,cot 
!                                                                       
      u1=dble(ux) 
      u2=dble(uy) 
      u3=dble(uz) 
      v1=dble(vx) 
      v2=dble(vy) 
      v3=dble(vz) 
!                                                                       
      dp= u1*v1 + u2*v2 + u3*v3 
      ru= dsqrt(u1*u1 + u2*u2 + u3*u3) 
      rv= dsqrt(v1*v1 + v2*v2 + v3*v3) 
      cot=dp/(ru*rv) 
      cot=cot -sign(1.d-13,cot) 
      ratio=real(ru/rv) 
                                                                        
      if(cot.ge.1.d0)  then 
                 print*, ' *** Rocotlib/cp_angle_and_ratio: cos > 1 !!!' 
                 print*, '                       angle set to 0.' 
                 angle=0. 
                 return 
                 endif 
                                                                        
      if(cot.lt.-1.d0) then 
                 print*, ' *** Rocotlib/cp_angle_and_ratio: cos < 1 !!!' 
                 print*, '                       angle set to 0.' 
                 angle=0. 
                 return 
                 endif 
!                                                                       
      angle=real(dacos(cot)) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_Euler_interp(a1,b1,c1,a2,b2,c2,ti,dt,ai,bi,ci) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_Euler_angles_interpolation                       
! *   Author : P. Robert, SDev, 2020                                    
!                                                                       
! *   Input  :  a1,b1,c1 Euler angles at time t1                        
!               a2,b2,c2 Euler angles at time t2                        
!               ti time for interpolation t1 < ti < t2                  
!               dt = t2 - t1                                            
!                                                                       
! *   Output :  ai,bi,ci Euler angles interpolated at time ti           
! ----------------------------------------------------------------------
!                                                                       
                                                                        
      Diff_a = a2 -a1 
      Diff_b = b2 -b1 
      Diff_c = c2 -c1 
                                                                        
!     WARNING: a does not vary linearly, it turns fast                  
!     while b and c vary slowly,                                        
!     so we can do the linear interpolation on 2 and 3                  
!     but not on a                                                      
                                                                        
      Diff_a= modulo(Diff_a,360.) 
                                                                        
!     for 2 and 3 management of zero crossings                          
                                                                        
      if(abs(Diff_b) > 180.) then 
              if(Diff_b < 0) then 
                                Diff_b= Diff_b +360. 
                                else 
                                Diff_b= Diff_b -360. 
              endif 
      endif 
                                                                        
      if(abs(Diff_c) > 180.) then 
              if(Diff_c < 0) then 
                                Diff_c= Diff_c +360. 
                                else 
                                Diff_c= Diff_c -360. 
              endif 
      endif 
                                                                        
!     interpolate Euler angles                                          
                                                                        
      ai = a1 +Diff_a*ti/dt 
      bi = b1 +Diff_b*ti/dt 
      ci = c1 +Diff_c*ti/dt 
                                                                        
      ai = modulo(a2,360.) 
      bi = modulo(b2,360.) 
      ci = modulo(c2,360.) 
                                                                        
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_geo_dipole_dir(iyear,idoy,d1,d2,d3) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_dipole_direction in GEO system                   
! *   Author : P. Robert, LPP , 2016                                    
!                                                                       
! *   Input  :  iyear (1900 - 2020), idoy= day of year (1/1=1)          
! *   Output :  d1,d2,d3  cartesian dipole components in GEO system     
! ----------------------------------------------------------------------
!                                                                       
      real year(25),rlat(25),rlon(25) 
                                                                        
      data iy,id,ipr /-1,-1,-1/ 
      save iy,id,ipr 
      save d1p,d2p,d3p 
!                                                                       
! *** table of geomagnetic Nort pole                                    
!     International Geomagnetic Reference Field: the 12th generation,   
!     Thébault et al. Earth, Planets and Space (2015) 67:79 ,           
!     DOI 10.1186/s40623-015-0228-9                                     
!                                                                       
!     geomagnetice pole (N & S) are symmetric, so there are used to     
!     define dipole axis. Do not mix up within the the North and south  
!     magnetic pole (not symetric).                                     
                                                                        
!                                                                       
      data year /1900.0, 1905.0, 1910.0, 1915.0, 1920.0, 1925.0, 1930.0,&
     &           1935.0, 1940.0, 1945.0, 1950.0, 1955.0, 1960.0, 1965.0,&
     &           1970.0, 1975.0, 1980.0, 1985.0, 1990.0, 1995.0, 2000.0,&
     &           2005.0, 2010.0, 2015.0, 2020.0/                        
!                                                                       
      data rlat / 78.68, 78.68, 78.66, 78.64, 78.63, 78.62, 78.60, 78.57&
     &          , 78.55, 78.55, 78.55, 78.54, 78.58, 78.60, 78.66, 78.76&
     &          , 78.88, 79.04, 79.21, 79.39, 79.61, 79.82, 80.09, 80.37&
     &          , 80.65/                                                
!                                                                       
      data rlon /-68.79,-68.75,-68.72,-68.57,-68.38,-68.27,-68.26,-68.36&
     &          ,-68.51,-68.53,-68.85,-69.16,-69.47,-69.85,-70.18,-70.47&
     &          ,-70.76,-70.90,-71.13,-71.42,-71.57,-71.81,-72.21,-72.63&
     &          ,-73.17/                                                
!                                                                       
!                        ------------------------                       
!                                                                       
! *** Computation are not done if date is the same as previous call     
!                                                                       
      if(iyear.eq.iy.and.idoy.eq.id) then 
        d1=d1p 
        d2=d2p 
        d3=d3p 
        return 
      endif 
!                                                                       
      iy=iyear 
      id=idoy 
!                                                                       
! *** Check date interval of validity                                   
!                                                                       
!     we are restricted by the interval 1965-2010, for which the igrf   
!     coefficients are known;                                           
!     if iyear is outside this interval, then the subroutine uses the   
!     nearest limiting value and prints a warning:                      
!                                                                       
      if(iy.lt.1900) then 
                     iy=1900 
                     if(ipr.ne.1) print 10, iyear, iy 
                     ipr=1 
                     alat=rlat(1) 
                     alon=rlon(1) 
                     endif 
                                                                        
      if(iy.gt.2020) then 
                     iy=2020 
                     if(ipr.ne.1) print 10, iyear, iy 
                     ipr=1 
                     alat=rlat(25) 
                     alon=rlon(25) 
                     endif 
!                                                                       
! *** decimal year                                                      
!                                                                       
      dyear=float(iy) +float(id)/365.25 
!                                                                       
! *** find two known intervals                                          
                                                                        
      do i=2,25 
         if(dyear.gt.year(i-1) .and. dyear.lt.year(i)) then 
             alat= rlat(i-1) + (rlat(i)-rlat(i-1))*(dyear-year(i-1))/5. 
             alon= rlon(i-1) + (rlon(i)-rlon(i-1))*(dyear-year(i-1))/5. 
             go to 20 
         endif 
      enddo 
                                                                        
   20 continue 
                                                                        
!                                                                       
! *** direction of dipole axis in GEO system:                           
!                                                                       
      clat=(90. -alat)*3.141593/180. 
      clon=alon*3.141593/180. 
                                                                        
      call t_sph_to_car(1.,clat,clon,d1,d2,d3) 
!                                                                       
      d1p=d1 
      d2p=d2 
      d3p=d3 
!                                                                       
   10 format(' * ROCOTLIB/cp_geo_dipole_dir: Warning! year=',i4.4,      &
     &       '   dipole direction can be computed between 1900-2020.',  &
     &       '   It will be computed for year ',i4.4)                   
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_gei_sun_dir(iyear,idoy,ih,im,is,                    &
     &                          gst,slong,sra,sdec,obliq)               
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_sun_direction in GEI system                      
! *   Author : CT.Russel, CE-D, 1971, rev. P.R., 1992,2001,2002         
!                                                                       
! *   Comment: calculates four quantities in gei system necessary for   
!              coordinate transformations dependent on sun position     
!              (and, hence, on universal time and season)               
!              Initial code from C.T. Russel, cosmic electro-dynamics,  
!              v.2, 184-196, 1971.                                      
!              accuracy: 0.006 degree.                                  
!              Adaptation P.Robert, November 1992.                      
!              Revised and F90 compatibility, P. Robert June 2001.      
!              Optimisation of DBLE computations and comments,          
!              P. Robert, December 2002                                 
!                                                                       
! *   Input  : iyear : year (1901-2099)                                 
!              idoy : day of the year (1 for january 1)                 
!              ih,im,is : hours, minutes, seconds U.T.                  
!                                                                       
! *   Output : gst      greenwich mean sideral time (radians)           
!              slong    longitude along ecliptic (radians)              
!              sra      right ascension (radians)                       
!              sdec     declination of the sun (radians)                
!              obliq    inclination of Earth's axis (radians)           
! ----------------------------------------------------------------------
!                                                                       
      double precision dj,fday 
!                                                                       
!                                                                       
      if(iyear.lt.1901.or.iyear.gt.2099) then 
      print*,'*** Rocotlib/cp_gei_sun_dir: year = ',iyear 
      print*,'*** Rocotlib/cp_gei_sun_dir: year must be >1900 and <2099' 
      stop   '*** Rocotlib/cp_gei_sun_dir: year must be >1901 and <3000' 
                                         endif 
!                                                                       
      pi=acos(-1.) 
      pisd= pi/180. 
!                                                                       
! *** Julian day and greenwich mean sideral time                        
!                                                                       
      fday=dble(ih*3600+im*60+is)/86400.d0 
      ileap=(iyear-1901)/4 
!     Note: year 2000 is a leap year, line above is right               
      dj=dble(365*(iyear-1900) +ileap +idoy) -0.5d0 +fday 
      gst=real(dmod(279.690983d0 +0.9856473354d0*dj+360.d0*fday +180.d0,&
     &               360.d0) )*pisd                                     
!                                                                       
! *** longitude along ecliptic                                          
!                                                                       
      vl= real( dmod(279.696678d0 +0.9856473354d0*dj,360.d0) ) 
      t=real(dj/36525.d0) 
      g=real( dmod(358.475845d0 +0.985600267d0*dj, 360.d0) )*pisd 
      slong=(vl+(1.91946 -0.004789*t)*sin(g) +0.020094*sin(2.*g))*pisd 
!                                                                       
! *** inclination of Earth's axis                                       
!                                                                       
      obliq=(23.45229 -0.0130125*t)*pisd 
      sob=sin(obliq) 
      cob=cos(obliq) 
!                                                                       
!     precession of declination (about 0.0056 deg., ajou P. Robert)     
!                                                                       
      pre= (0.0055686 - 0.025e-4*t)*pisd 
!                                                                       
! *** declination of the sun                                            
!                                                                       
      slp=slong -pre 
      sind=sob*sin(slp) 
      cosd=sqrt(1. -sind**2) 
      sc=sind/cosd 
      sdec=atan(sc) 
!                                                                       
! *** right ascension of the sun                                        
!                                                                       
      sra=pi -atan2((cob/sob)*sc, -cos(slp)/cosd) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_sunrise_sunset(iyear,imon,iday,rlat,rlon,           &
     &                   tmer,tris,tset,durd,                           &
     &                   elemer,azimer,eleris,aziris,eleset,aziset,icor)
                                                                        
      common /sunVDH/ suntim(1441),sunele(1441),sunazi(1441) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_sunset_time and others                           
! *   Author : P. Robert, CRPE, 2001 Revised Dec. 2011                  
! *   Comment: Spherical Earth assumed                                  
!                                                                       
! *   Input  : iyear,imon,iday, rlat, rlong (rad)                       
!              icor: type de correction                                 
!                  0 pas de correction                                  
!                  1: correction de la refraction atmospherique seule   
!                     (calcul comparable a l'IMC, ancien BDL)           
!                  2: correction de la refraction et du bord superieur  
!                     du soleil (calcul de la SAF).                     
!                  3: correction de la refraction avec l'elevation      
!                     pour le passage au meridien (utile pour les zones 
!                     polaires, ou le soleil peut raser l'horizon       
!                     a midi)                                           
!                                                                       
!                  voir explications plus precises ci-dessous           
!                                                                       
! *   Output : tmer,tris,tset,durd on char*8 format, as '23:42:37'      
!              tmer: Sun Meridian time (local noon)                     
!              tris: Sunrise time                                       
!              tset: Sunset  time                                       
!              durd: Duration of the day                                
!                                                                       
!              elemer; elevation, en degres, du soleil au meridin       
!              azimer: azimuth, en degres, a partir du nord vers l'est  
!                      du soleil au meridien                            
!              eleris,aziris: meme chose pour le lever de soleil        
!              eleset,aziset: meme chose pour le coucher                
!                                                                       
!              /sunVDH/ trajectoire du Soleil pendant la journée        
!              Accuracy: 10 sec.                                        
!                                                                       
!    Les valeures pour Paris (latitude moyenne) et Tromso (zone polaire)
!    ont ete comparees aux ephemerides de l'institut de mecanique       
!    celeste (anciennement bureau des longitudes) avec icor=1           
!    Les resultats sont comparables a la minute près pour les levers    
!    et couchers, et a quelques secondes d'arc pour les elevations et   
!    azimuts, y compris dans les zones polaires et pour les transitions 
!    avec les nuits polaires ou jours polaires (pas de lever ni coucher)
!                                                                       
! ----------------------------------------------------------------------
!                                                                       
! *** calcul des heures de lever et coucher de soleil:                  
!     correction du diametre apparent du soleil et de la refraction     
!     a l'horizon. Les valeurs prises sont extraites de l'ouvrage       
!     "introduction aux ephemerides astronomiques", publie par le       
!     Bureau Des Longitudes.                                            
!     Diametre apparent du soleil: 32'                                  
!     Refraction a l'horizon     : 34' (36.6 pour le BDL)               
!     correction= 32'/2 + 34' = 50'                                     
!                                                                       
! *** Attention :                                                       
!     Les Ephemrides Astronomiques du Bureau des Longitudes ne tiennent 
!     compte que de la refraction a laquelle elles donnent la valeur    
!     r = 36,6'.                                                        
!     Les Ephemrides Astronomiques de la S.A.F. considere le lever comme
!     l'apparition de son bord superieur.                               
!     Elles utilisent la valeur ht = -50'.                              
!                                                                       
! *   variation avec la temperature et la pression:                     
!     on peut utiliser la formule de Bennet, pour corriger la refraction
!     en la multipliant par le coefficient:                             
!     (P/1010)(283/(273+T))                                             
!     avec P en millibar et T en Celsius.                               
!     Note: Le BDL prend 1 pour cette correction, donc                  
!     suppose une presion P=1010 et une temperature T=10 deg.           
!                                                                       
! *   variation de R avec la hauteur apparente:                         
!     Une formule assez simple est donne par Bennet :                   
!     R = 1/Tan[h+7,31/(h+4,4)]                                         
!     avec h en degre et R en minute d'arc                              
!     soit pour h=0 on trouve R= 34.5                                   
!          pour h=45          R=  0.99'                                 
!     (34' et 1' etant les valeurs couramment admises)                  
!                                                                       
! *   Crepuscules : par definition la fin (le soir) ou le debut         
!     (le matin) des crepuscules civil, nautique et astronomique        
!     se produit quand le centre du Soleil est abaisse de 6, 12 et 18   
!     degres sous l'horizon.                                            
! ----------------------------------------------------------------------
!                                                                       
      character*(*) tmer,tris,tset,durd 
      dimension mimer(2), vemer(2) 
!                                                                       
      lmer=len(tmer) 
      lris=len(tris) 
      lset=len(tset) 
      ldur=len(durd) 
!                                                                       
      if(lmer.lt.8) then 
           print*, '*** Rocotlib/cp_sunrise_sunset: tmer must be char*8' 
           print*, '    len=',lmer 
           stop 'cp_sunrise_sunset: tmer must be char*8 variable' 
                    endif 
!                                                                       
      if(lris.lt.8) then 
           print*, '*** Rocotlib/cp_sunrise_sunset: tris must be char*8' 
           print*, '    len=',lris 
           stop 'cp_sunrise_sunset: tris must be char*8 variable' 
                    endif 
!                                                                       
      if(lset.lt.8) then 
           print*, '*** Rocotlib/cp_sunrise_sunset: tset must be char*8' 
           print*, '    len=',lset 
           stop 'cp_sunrise_sunset: tset must be char*8 variable' 
                    endif 
!                                                                       
      if(ldur.lt.8) then 
           print*, '*** Rocotlib/cp_sunrise_sunset: durd must be char*8' 
           print*, '    len=',ldur 
           stop 'cp_sunrise_sunset: durd must be char*8 variable' 
                    endif 
!                                                                       
      milris= -1 
      milset= -1 
      milmer= -1 
      milmin= -1 
                                                                        
      imer=0 
      vemer(1)= -2. 
      vemer(2)= -2. 
!                                                                       
      sdp=0. 
      svp=0. 
                                                                        
      deg=180./acos(-1.) 
!                                                                       
! *   tolerance de 1 degre sur le meridien quand le soleil est          
!     au ras de l'horizon et donc sujet a la refraction                 
!                                                                       
      cormax=1. 
      svmin= -sin(cormax/deg) 
                                                                        
! *** initialisation de sdp=sd et svp=sv pour le jour d'avant           
                                                                        
      call cv_date_to_jul2000(iyear,imon,iday,jd00) 
      jd00=jd00-1 
      call cv_jul2000_to_date(jd00,iyearp,imonp,idayp) 
      call cp_time_param(iyearp,imonp,idayp,23,59,30) 
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
      call t_geo_to_vdh(sxgeo,sygeo,szgeo,rlat,rlon,sv,sd,sh) 
                                                                        
! *   correction de la refraction (voir + loin)                         
                                                                        
      call t_car_to_sph(sd,sh,sv,r,tetnc,phi) 
! *   pas de correction de refraction par defaut                        
      tet=tetnc 
! *   calcul de la refraction a l'horizon                               
      if(icor.gt.0) then 
                    if(icor.eq.1) then 
                                  cor= 36.6/60. 
                                  else 
                                  cor= 50./60. 
                    endif 
                    tet=tetnc -cor/deg 
                    call t_sph_to_car(r,tet,phi,sd,sh,sv) 
      endif 
                                                                        
      sdp=sd 
      svp=sv 
      iyearp=iyear 
      imonp =imon 
      idayp =iday 
      ismax=50 
!                                                                       
!     ----------------------------------------------                    
!     boucle sur le temps de la journee                                 
!     ----------------------------------------------                    
                                                                        
      do 10 imin=0,1440 
                                                                        
        if (imin.eq.1440) then 
          call cv_date_to_jul2000(iyear,imon,iday,jd00) 
          jd00=jd00+1 
          call cv_jul2000_to_date(jd00,iyear,imon,iday) 
          ismax=30 
      endif 
                                                                        
      do  8 is  =0,ismax,10 
                                                                        
      if(imin.lt.1440) then 
                       isec=imin*60+is 
                       else 
                       isec=is 
      endif 
!                                                                       
! *** calcul de la direction du Soleil dans le VDH pour chaque pas de   
!     temps de la journee en cours ( precision 10 sec)                  
                                                                        
      millsec=isec*1000 
                                                                        
      call cv_msotd_to_hmsms(millsec,ih,im,is,ims) 
      call cp_time_param(iyear,imon,iday,ih,im,is) 
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
      call t_geo_to_vdh(sxgeo,sygeo,szgeo,rlat,rlon,sv,sd,sh) 
!                                                                       
! *** calcul des heures de lever et coucher de soleil:                  
!     correction du diametre apparent du soleil et de la refraction     
!     a l'horizon. Les valeurs prises sont extraites de l'ouvrage       
!     "introduction aux ephemerides astronomiques", publie par le       
!     Bureau Des Longitudes.                                            
!     Diametre apparent du soleil: 32'                                  
!     Refraction a l'horizon     : 34'                                  
!     correction= 32'/2 + 34' = 50'                                     
!                                                                       
! *** Attention :                                                       
!     Les Ephemrides Astronomiques du Bureau des Longitudes ne tiennent 
!     compte que de la refraction a laquelle elles donnent la valeur    
!     r = 36,6'.                                                        
!     Les Ephemrides Astronomiques de la S.A.F. considere le lever comme
!     l'apparition de son bord superieur.                               
!     Elles utilisent la valeur ht = -50'.                              
!                                                                       
!                                                                       
! *   dans le reperes DHV, teta est l'angle entre la verticale sortante 
!     et la direction du Soleil                                         
                                                                        
      call t_car_to_sph(sd,sh,sv,r,tetnc,phi) 
                                                                        
! *   pas de correction de refraction par defaut                        
                                                                        
      tet=tetnc 
                                                                        
! *   calcul de la refraction a l'horizon                               
                                                                        
      if(icor.gt.0) then 
                    if(icor.eq.1) then 
                                  cor= 36.6/60. 
                                  else 
                                  cor= 50./60. 
                    endif 
                                                                        
                    tet=tetnc -cor/deg 
                    call t_sph_to_car(r,tet,phi,sd,sh,sv) 
                                                                        
! *   calcul de la refraction a teta pour le calcul de la trajectoire   
!     Quand le Soleil est vertical (tet=0), la correction est minime;   
!     Quand le soleil est proche de l'horizon (tet=90) cor=34'          
                                                                        
                    if(icor.eq.3) then 
                                  h  = 90. -tet*deg 
                                  cor= 1./tan( (h +7.31/(h+4.4))/deg ) 
                                  cor=cor/60. 
                    endif 
                                                                        
                    tet=tetnc -cor/deg 
      endif 
                                                                        
! *** chargement du common /sunVDH/                                     
!     resolution en temps: 1mn                                          
!     l'azimut est compté positivement depuis le Nors vers l'Est        
                                                                        
      if(is.eq.0) then 
          suntim(imin+1)= float(imin) 
          sunele(imin+1)= 90. -tet*deg 
          sunazi(imin+1)= 90. -phi*deg 
      endif 
!                                                                       
! *** calcul du midi local, quand D (vers l'est) change de signe        
!     attention: dans le cas du soleil de minuit, il y a deux meridien  
!     ou le soleil est visible.                                         
!     On prend celui ou le soleil est au plus haut                      
!     (a midi, et non pas a minuit)                                     
!                                                                       
      if(sv.gt.svmin.and.sd*sdp.lt.0.) then 
           imer=imer+1 
           if(imer.gt.2)then 
                      print*,'**Rocotlib more than 2 meridian crossing!' 
                      print*,'           last is taken into account' 
                        go to 21 
           endif 
                                                                        
           mimer(imer)=isec*1000 
           vemer(imer)=sv 
      endif 
                                                                        
   21 continue 
!                                                                       
! *** calcul du temps ou l'elevation est minimum (soleil de minuit)     
!                                                                       
      if(sv.ge.0. .and. svp.ge.0. .and. sv.lt.svp) milmin= isec*1000 
!                                                                       
! *** on a une transition (lever ou couche) quand V (vertical sortante) 
!     change de signe (unite: millisec du jour)                         
!                                                                       
! *   cas ou sv est nul (le soleil est a l'horizon)                     
!                                                                       
      epsi=1.e-6 
                                                                        
!cc   if(sv.eq.0.) then                                                 
      if(abs(sv).lt.epsi) then 
                   sv=0. 
                   if(svp.lt.0.) then 
!                    juste avant, le soleil etait sous l'horizon        
!                    la transition est donc un lever de soleil          
                                 milris= isec*1000 
                                 else 
!                    juste avant, le soleil etait au dessus de l'horizon
!                    la transition est donc un coucher de soleil        
                                 milset= isec*1000 
                    endif 
!cc   print*, '          cas 1 sv,svp=',sv,svp,' ris,set=',milris,milset
      go to 20 
      endif 
                                                                        
! *   cas ou svp est nul (le soleil etait a l'horizon)                  
!                                                                       
      if(abs(svp).lt.epsi) then 
                   if(sv.gt.0.) then 
!                    le soleil est au dessus de l'horizon               
!                    la transition est donc un lever de soleil          
                                 milris= isec*1000 
                                 else 
!                    le soleil est au dessous de l'horizon              
!                    la transition est donc un coucher de soleil        
                                 milset= isec*1000 
                    endif 
!cc   print*, '          cas 2 sv,svp=',sv,svp,' ris,set=',milris,milset
      go to 20 
      endif 
!                                                                       
! *   sv et svp sont non nuls                                           
!                                                                       
      if(sv*svp.lt.0.) then 
!              si sv et svp sont de signe contraire, on a une transition
!              (lever ou coucher)                                       
                       if(sv.gt.0.) then 
!                    le soleil passe au passe au dessus de l'horizon    
!                    la transition est donc un lever de soleil          
                                    milris= isec*1000 
                                    else 
!                    le soleil passe au passe au dessous de l'horizon   
!                    la transition est donc un coucher de soleil        
                                    milset= isec*1000 
                       endif 
!cc   print*, '         cas 3 sv,svp=',sv,svp,' ris,set=',milris,milset 
      endif 
                                                                        
                                                                        
   20 continue 
      sdp=sd 
      svp=sv 
                                                                        
    8 continue 
   10 continue 
!                                                                       
!     ----------------------------------------------                    
!     la boucle sur le temps de la journee est finie                    
!     ----------------------------------------------                    
                                                                        
! *   on restaure le jour en cours                                      
                                                                        
      iyear=iyearp 
      imon =imonp 
      iday =idayp 
                                                                        
! *** encodage des resultats sous la forme HH:MM:SS                     
!     et calcul de l'elevation et de l'azimuth pour chacun des 3 cas    
!                                                                       
      tmer='night   ' 
      tris='no SR   ' 
      tset='no SS   ' 
      durd='00:00:00' 
                                                                        
      elemer=999. 
      eleris=999. 
      eleset=999. 
                                                                        
      azimer=999. 
      aziris=999. 
      aziset=999. 
!                                                                       
! *   choix du midi local si soleil de minuit                           
!                                                                       
      if(imer.gt.0) then 
                    if(vemer(1).gt.vemer(2)) then 
                                             milmer=mimer(1) 
                                             else 
                                             milmer=mimer(2) 
                    endif 
      endif 
!                                                                       
! *   calcul des positions au temps du meridien                         
!                                                                       
      if(milmer.gt.-1) then 
            call cv_msotd_to_hmsms(milmer,ih,im,is,ims) 
            write(tmer,100) ih,im,is 
                                                                        
            call cp_time_param(iyear,imon,iday,ih,im,is) 
            call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
            call t_geo_to_vdh(sxgeo,sygeo,szgeo,rlat,rlon,sv,sd,sh) 
            call t_car_to_sph(sd,sh,sv,r,tet,phi) 
!                                                                       
            if(icor.eq.3) then 
!                         calcul de la refraction a teta                
                          h  = 90. -tet*deg 
                          cor= 1./tan( (h +7.31/(h+4.4))/deg) 
                          cor=cor/60. 
                          tet=tet -cor/deg 
            endif 
                                                                        
            elemer= 90. -tet*deg 
            azimer= 90. -phi*deg 
      endif 
!                                                                       
! *   duree du jour calculee seuleument si le lever et le coucher       
!     sont definis                                                      
!                                                                       
      if(milris .gt. -1 .and. milset .gt. -1) then 
                 idiff= milset-milris 
                 if(idiff.lt.0) idiff=idiff +86400000 
                 call cv_msotd_to_hmsms(idiff,ih,im,is,ims) 
                 write(durd,200) ih,im,is 
                                            else 
                 if(elemer.gt.0. .and. elemer.lt.180.) durd='24:00:00' 
      endif 
!                                                                       
  100 format(i2.2,':',i2.2,':',i2.2) 
  200 format(i2  ,':',i2.2,':',i2.2) 
!                                                                       
! *** cas du soleil de minuit: lever = coucher= temp elevation min      
!                                                                       
      if(durd.eq.'24:00:00' .and. milmin.gt.-1)  then 
                      milris=milmin 
                      milset=milmin 
      endif 
!                                                                       
      if(milris.gt.-1) then 
            call cv_msotd_to_hmsms(milris,ih,im,is,ims) 
            write(tris,100) ih,im,is 
                                                                        
            call cp_time_param(iyear,imon,iday,ih,im,is) 
            call g_gei_geo_sun_dir(sxgei,sygei,szgei, sxgeo,sygeo,szgeo) 
            call t_geo_to_vdh(sxgeo,sygeo,szgeo,rlat,rlon,sv,sd,sh) 
            call t_car_to_sph(sd,sh,sv,r,tet,phi) 
                                                                        
            eleris= 90. -tet*deg 
            aziris= 90. -phi*deg 
      endif 
!                                                                       
      if(milset.gt.-1) then 
            call cv_msotd_to_hmsms(milset,ih,im,is,ims) 
            write(tset,100) ih,im,is 
                                                                        
            call cp_time_param(iyear,imon,iday,ih,im,is) 
            call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
            call t_geo_to_vdh(sxgeo,sygeo,szgeo,rlat,rlon,sv,sd,sh) 
            call t_car_to_sph(sd,sh,sv,r,tet,phi) 
                                                                        
            eleset= 90. -tet*deg 
            aziset= 90. -phi*deg 
      endif 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_time_param(iyear,imonth,iday,ih,im,is) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_time_parameters and time-dependent matrix        
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Comment: Prepare all time varying quantities for computations of  
!              coordinate transforms of the library.                    
!                                                                       
! *   Input  : iyear,imonth,iday (1901<year<2099)                       
!              ih,im,is : hours, minutes, seconds U.T.                  
! *   Output : in common statements                                     
! ----------------------------------------------------------------------
!                                                                       
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq 
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3,pr1,pr2,pr3 
!                                                                       
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3 
      common /timp04/ grgs1,grgs2,grgs3,  grge1,grge2,grge3 
      common /timp05/ gdps1,gdps2,gdps3 
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3 
!                                                                       
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3 
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3 
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3 
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3 
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3 
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3 
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3 
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3 
!                                                                       
      data idoy /-1/ 
!                                                                       
! *** day of year                                                       
!                                                                       
      call cv_date_to_doty(iyear,imonth,iday,idoy) 
!                                                                       
! *** Dipole direction in GEO system                                    
!                                                                       
      call cp_geo_dipole_dir(iyear,idoy,gd1,gd2,gd3) 
!                                                                       
      qd12=sqrt(gd1**2 + gd2**2) 
!                                                                       
! *** Sun direction in GEI (from right ascension and declination),      
!     greenwich mean sideral time and longitude along ecliptic          
!                                                                       
      call cp_gei_sun_dir(iyear,idoy,ih,im,is,                          &
     &                    gst,slong,srasn,sdecl,obliq)                  
!                                                                       
      gs1=cos(srasn)*cos(sdecl) 
      gs2=sin(srasn)*cos(sdecl) 
      gs3=sin(sdecl) 
!                                                                       
! *** sin and cos of GMST                                               
!                                                                       
      sgst=sin(gst) 
      cgst=cos(gst) 
!                                                                       
! *** ecliptic pole in GEI system                                       
!                                                                       
      ge1=  0. 
      ge2= -sin(obliq) 
      ge3=  cos(obliq) 
!                                                                       
! *** direction of the rotation axis of the sun in GEI system           
!     (from C.T. Russell, Cosmic Electro-dynamics, V.2, 184-196, 1971.) 
!                                                                       
      rad=  57.29578 
!                                                                       
      ras= -74.0/rad 
      dec=  63.8/rad 
!                                                                       
      gr1=cos(ras)*cos(dec) 
      gr2=sin(ras)*cos(dec) 
      gr3=sin(dec) 
!                                                                       
! *** dipole direction in GEI system                                    
!                                                                       
      gm1= gd1*cgst - gd2*sgst 
      gm2= gd1*sgst + gd2*cgst 
      gm3= gd3 
!                                                                       
      qm12=sqrt(gm1**2 + gm2**2) 
!                                                                       
! *** direction of the sun in GEO system                                
!                                                                       
      ps1=  gs1*cgst + gs2*sgst 
      ps2= -gs1*sgst + gs2*cgst 
      ps3=  gs3 
!                                                                       
! *** direction of the ecliptic in GEO system                           
!                                                                       
      pe1=  ge1*cgst + ge2*sgst 
      pe2= -ge1*sgst + ge2*cgst 
      pe3=  ge3 
!                                                                       
! *** direction of the rotation axis of the sun in GEO system           
!                                                                       
      pr1=  gr1*cgst + gr2*sgst 
      pr2= -gr1*sgst + gr2*cgst 
      pr3=  gr3 
!                                                                       
! *** cross product MxS in GEI system                                   
!                                                                       
      gmgs1= gm2*gs3 - gm3*gs2 
      gmgs2= gm3*gs1 - gm1*gs3 
      gmgs3= gm1*gs2 - gm2*gs1 
!                                                                       
      rgmgs=sqrt(gmgs1**2 + gmgs2**2 + gmgs3**2) 
!                                                                       
! *** cross product ExS in GEI system                                   
!                                                                       
      gegs1= ge2*gs3 - ge3*gs2 
      gegs2= ge3*gs1 - ge1*gs3 
      gegs3= ge1*gs2 - ge2*gs1 
!                                                                       
! *** cross product RxS in GEI system                                   
!                                                                       
      grgs1= gr2*gs3 - gr3*gs2 
      grgs2= gr3*gs1 - gr1*gs3 
      grgs3= gr1*gs2 - gr2*gs1 
!                                                                       
      rgrgs=sqrt(grgs1**2 + grgs2**2 + grgs3**2) 
!                                                                       
! *** cross product RxE in GEI system                                   
!                                                                       
      grge1= gr2*ge3 - gr3*ge2 
      grge2= gr3*ge1 - gr1*ge3 
      grge3= gr1*ge2 - gr2*ge1 
!                                                                       
! *** cross product DxS in GEO system                                   
!                                                                       
      gdps1= gd2*ps3 - gd3*ps2 
      gdps2= gd3*ps1 - gd1*ps3 
      gdps3= gd1*ps2 - gd2*ps1 
!                                                                       
      rgdps=sqrt(gdps1**2 + gdps2**2 + gdps3**2) 
!                                                                       
! *** cross product ExS in GEO system                                   
!                                                                       
      peps1= pe2*ps3 - pe3*ps2 
      peps2= pe3*ps1 - pe1*ps3 
      peps3= pe1*ps2 - pe2*ps1 
!                                                                       
! *** cross product RxS in GEO system                                   
!                                                                       
      prps1= pr2*ps3 - pr3*ps2 
      prps2= pr3*ps1 - pr1*ps3 
      prps3= pr1*ps2 - pr2*ps1 
!                                                                       
      rprps=sqrt(prps1**2 + prps2**2 + prps3**2) 
!                                                                       
! *** computation of gei to mag vectors                                 
!                                                                       
      xeima1=  gm1*gm3/qm12 
      xeima2=  gm2*gm3/qm12 
      xeima3= -qm12 
!                                                                       
      yeima1= -gm2/qm12 
      yeima2=  gm1/qm12 
      yeima3=  0. 
!                                                                       
! *** computation of gei to sm vectors                                  
!                                                                       
      yeism1=gmgs1/rgmgs 
      yeism2=gmgs2/rgmgs 
      yeism3=gmgs3/rgmgs 
!                                                                       
      xeism1= yeism2*gm3 - yeism3*gm2 
      xeism2= yeism3*gm1 - yeism1*gm3 
      xeism3= yeism1*gm2 - yeism2*gm1 
!                                                                       
! *** computation of gei to gsm vectors                                 
!                                                                       
      yeigm1= gmgs1/rgmgs 
      yeigm2= gmgs2/rgmgs 
      yeigm3= gmgs3/rgmgs 
!                                                                       
      zeigm1= gs2*yeigm3 - gs3*yeigm2 
      zeigm2= gs3*yeigm1 - gs1*yeigm3 
      zeigm3= gs1*yeigm2 - gs2*yeigm1 
!                                                                       
! *** computation of gei to gseq vectors                                
!                                                                       
      yeigq1= grgs1/rgrgs 
      yeigq2= grgs2/rgrgs 
      yeigq3= grgs3/rgrgs 
!                                                                       
      zeigq1= gs2*yeigq3 - gs3*yeigq2 
      zeigq2= gs3*yeigq1 - gs1*yeigq3 
      zeigq3= gs1*yeigq2 - gs2*yeigq1 
!                                                                       
! *** computation of tetq angle                                         
!                                                                       
      stetq= (grge1*gs1 + grge2*gs2 + grge3*gs3)/rgrgs 
      ctetq= sqrt(1.-stetq**2) 
!                                                                       
! *** computation of mu angle                                           
!                                                                       
      smu= ps1*gd1 + ps2*gd2 + ps3*gd3 
      cmu= sqrt(1.-smu*smu) 
!                                                                       
! *** computation of dzeta angle                                        
!                                                                       
      cdze= (ge1*gm1   + ge2*gm2   + ge3*gm3)/rgmgs 
      sdze= (ge1*gmgs1 + ge2*gmgs2 + ge3*gmgs3)/rgmgs 
      epsi=1.e-5 
      if(abs(sdze**2 +cdze**2 -1.).gt.epsi) then
                                           print*, 'sdze,cdze=',sdze,cdze, sdze**2 +cdze**2
                                           stop '*** Rocotlib error 3' 
      endif
!                                                                       
! *** computation of phi angle                                          
!                                                                       
      q=qd12*rgdps 
!                                                                       
      cphi=(gd1*gd3*ps1 + gd2*gd3*ps2 - (gd1**2+gd2**2)*ps3)/q 
      sphi=(gd2*ps1 - gd1*ps2)/q 
      if(abs(sphi**2 +cphi**2-1.).gt.epsi) then
                                           print*,'sphi,cphi=',sphi,cphi
                                           stop '*** Rocotlib error 4' 
      endif
!                                                                       
! *** computation of geo to mag vectors                                 
!                                                                       
      yeoma1= -gd2/qd12 
      yeoma2=  gd1/qd12 
      yeoma3=  0. 
!                                                                       
      xeoma1=  yeoma2*gd3 
      xeoma2= -yeoma1*gd3 
      xeoma3=  yeoma1*gd2 - yeoma2*gd1 
!                                                                       
! *** computation of geo to sm vectors                                  
!                                                                       
      yeosm1= gdps1/rgdps 
      yeosm2= gdps2/rgdps 
      yeosm3= gdps3/rgdps 
!                                                                       
      xeosm1= yeosm2*gd3 - yeosm3*gd2 
      xeosm2= yeosm3*gd1 - yeosm1*gd3 
      xeosm3= yeosm1*gd2 - yeosm2*gd1 
!                                                                       
! *** computation of geo to gsm vectors                                 
!                                                                       
      yeogm1=yeosm1 
      yeogm2=yeosm2 
      yeogm3=yeosm3 
!                                                                       
      zeogm1= ps2*yeogm3 - ps3*yeogm2 
      zeogm2= ps3*yeogm1 - ps1*yeogm3 
      zeogm3= ps1*yeogm2 - ps2*yeogm1 
!                                                                       
! *** computation of geo to gsq vectors                                 
!                                                                       
      yeogq1= prps1/rprps 
      yeogq2= prps2/rprps 
      yeogq3= prps3/rprps 
!                                                                       
      zeogq1= ps2*yeogq3 - ps3*yeogq2 
      zeogq2= ps3*yeogq1 - ps1*yeogq3 
      zeogq3= ps1*yeogq2 - ps2*yeogq1 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_time_param2(jd1950,houday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_time_parameters and time-dependent matrix        
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Comment: Prepare all time varying quantities for computations of  
!              coordinate transforms of the library.                    
!              Same as cp_time_param, only input arguments are changed. 
!                                                                       
! *   Input  : jd1950: Julian day 1950 (0= 1/1/1950)                    
!              houday: decimal hour of the day (U.T.)                   
! *   Output : in common statements                                     
! ----------------------------------------------------------------------
!                                                                       
      data iyear,imonth,iday /-1,-1,-1/ 
      data ih,im,is /-1,-1,-1/ 
!                                                                       
      call cv_jul1950_to_date(jd1950,iyear,imonth,iday) 
      call cv_dech_to_hms(houday,ih,im,is) 
!                                                                       
      call cp_time_param(iyear,imonth,iday,ih,im,is) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_time_param3(jd2000,houday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_time_parameters and time-dependent matrix        
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Comment: Prepare all time varying quantities for computations of  
!              coordinate transforms of the library.                    
!              Same as cp_time_param, only input arguments are changed. 
!                                                                       
! *   Input  : jd2000: Julian day 2000 (0= 1/1/2000)                    
!              houday: decimal hour of the day (U.T.)                   
! *   Output : in common statements                                     
! ----------------------------------------------------------------------
!                                                                       
      data iyear,imonth,iday /-1,-1,-1/ 
      data ih,im,is /-1,-1,-1/ 
!                                                                       
      call cv_jul1950_to_date(jd2000,iyear,imonth,iday) 
      call cv_dech_to_hms(houday,ih,im,is) 
!                                                                       
      call cp_time_param(iyear,imonth,iday,ih,im,is) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_tpn_param(xo,yo,zo,xs, Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz) 
!                                                                       
      real Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz 
                                                                        
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_TPN_system                                       
! *   Author : P. Robert, CETP, 2004                                    
! *   Comment: Compute TPN vector in GSE system or any system having X  
!              axis towards the SUN.                                    
!                                                                       
!         N: Output normal to the paraboloid                            
!         T: tengente to the paraboloid, towards the summit             
!         P: tengente to the paraboloid, P=N X T                        
!                                                                       
!     The paraboloid is defined by its summit Xs and the local point    
!     Xo, Yo, Zo                                                        
!                                                                       
!     Note: the paraboloid is close to the magnetopause if the summit   
!           is defined as the subsolar point (by T87,T89 model or other)
!           and if the local point Xo,Yo,Zo correspond to magnetopause  
!           crossing.                                                   
!                                                                       
! *   Input : xo,yo,zo,xs                                               
! *   Output: Nx,Ny,Nz,Tx,Ty,Tz,Px,Py,Pz                                
!     ---------------------------------------------------------------+--
!                                                                       
!                                                                       
      r= sqrt(yo**2 +zo**2) 
      dx=xs-xo 
!                                                                       
      Nx= r 
      Ny= 2.*dx*yo/r 
      Nz= 2.*dx*zo/r 
!                                                                       
      Tx= 2.*dx 
      Ty= -yo 
      Tz= -zo 
!                                                                       
      Px= Ny*Tz -Nz*Ty 
      Py= Nz*Tx -Nx*Tz 
      Pz= Nx*Ty -Ny*Tx 
!                                                                       
!     normalisation                                                     
!                                                                       
      rnoN= sqrt(Nx**2 +Ny**2 +Nz**2) 
      rnoT= sqrt(Tx**2 +Ty**2 +Tz**2) 
      rnoP= sqrt(Px**2 +Py**2 +Pz**2) 
!                                                                       
      Nx=Nx/rnoN 
      Ny=Ny/rnoN 
      Nz=Nz/rnoN 
!                                                                       
      Tx=Tx/rnoT 
      Ty=Ty/rnoT 
      Tz=Tz/rnoT 
!                                                                       
      Px=Px/rnoP 
      Py=Py/rnoP 
      Pz=Pz/rnoP 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_nbday_in_month(iyear,imonth,nbday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : compute_number_of_day_of_the_month                       
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : iyear,imonth (1-12)                                      
! *   Output : nbday                                                    
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      dimension nday(12) 
      data nday/31,28,31,30,31,30,31,31,30,31,30,31/ 
      data ily /-1/ 
!                                                                       
      if(imonth.lt.1.or.imonth.gt.12) then 
       print*,'*** Rocotlib/cp_nbday_in_month: month= ',imonth 
       print*,'*** Rocotlib/cp_nbday_in_month: month must be >0 and <13' 
       stop   '*** Rocotlib/cp_nbday_in_month: month must be >0 and <13' 
                                      endif 
!                                                                       
      call cp_leap_year(iyear,ily) 
!                                                                       
      if(ily.eq.1) then 
                   nday(2)=29 
                   else 
                   nday(2)=28 
                   endif 
!                                                                       
      nbday=nday(imonth) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_en_day_name(iday,cday,nbcha) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : compute_english_day_name, ex: 'Monday' for iday=1        
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : iday (1-7, otherwise modulo is done)                     
! *   Output : cday,nchar                                               
! ----------------------------------------------------------------------
!                                                                       
      character*(*) cday 
      character*9 days (7) 
      integer nbca(7) 
!                                                                       
      data days /'Monday','Tuesday','Wednesday','Thurday','Friday',     &
     &           'Saturday','Sunday'/                                   
      data nbca /6,7,9,7,6,8,6/ 
!                                                                       
      iweek= iday/7 
      nday= iday -7*iweek 
      if(nday.lt.1) nday=nday+7 
!                                                                       
      cday= days(nday) 
      nbcha=nbca(nday) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_en_month_name(imonth,cmonth,nchar) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : compute_english_month_name                               
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : imonth (1-12)                                            
! *   Output : cmonth,nchar                                             
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      integer nbca(12) 
      character*9 cara(12) 
      character*(*) cmonth 
!                                                                       
      data cara/'January','February','March','April','May','June',      &
     &     'July','August','September','October','November','December'/ 
      data nbca /7,8,5,5,3,4,4,6,9,7,8,8/ 
!                                                                       
!                                                                       
      if(imonth.lt.1.or.imonth.gt.12) then 
        print*,'*** Rocotlib/cp_en_month_name: month= ',imonth 
        print*,'*** Rocotlib/cp_en_month_name: month must be >0 and <13' 
        stop   '*** Rocotlib/cp_en_month_name: month must be >0 and <13' 
                                      endif 
!                                                                       
      cmonth=cara(imonth) 
      nchar =nbca(imonth) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_fr_day_name(iday,cday,nbcha) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : compute_french_day_name, ex: 'Lundi' for iday=1          
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : iday (1-7, otherwise modulo is done)                     
! *   Output : cday,nchar                                               
! ----------------------------------------------------------------------
!                                                                       
      character*(*) cday 
      character*8 days (7) 
      integer nbca(7) 
!                                                                       
      data days /'Lundi','Mardi','Mercredi','Jeudi','Vendredi',         &
     &           'Samedi','Dimanche'/                                   
      data nbca /5,5,8,5,8,6,8/ 
!                                                                       
      iweek= iday/7 
      nday= iday -7*iweek 
      if(nday.lt.1) nday=nday+7 
!                                                                       
      cday= days(nday) 
      nbcha=nbca(nday) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_fr_month_name(imonth,cmonth,nchar) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : compute_french_month_name                                
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : imonth (1-12)                                            
! *   Output : cmonth,nchar                                             
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      integer nbca(12) 
      character*9 cara(12) 
      character*(*) cmonth 
!                                                                       
      data cara /'Janvier','Février','Mars','Avril','Mai','Juin',       &
     &           'Juillet','Août','Septembre','Octobre','Novembre',     &
     &           'Décembre'/                                            
!                                                                       
      data nbca /7,7,4,5,3,4,7,4,9,7,8,8/ 
!                                                                       
!                                                                       
      if(imonth.lt.1.or.imonth.gt.12) then 
        print*,'*** Rocotlib/cp_fr_month_name: month= ',imonth 
        print*,'*** Rocotlib/cp_fr_month_name: month must be >0 and <13' 
        stop   '*** Rocotlib/cp_fr_month_name: month must be >0 and <13' 
                                      endif 
!                                                                       
      cmonth=cara(imonth) 
      nchar =nbca(imonth) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_leap_year(iyear,ileap) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : compute_leap_year with ileap=1 for leap year, 0 if not   
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : iyear (ex: 1980)                                         
! *   Output : ileap (1 or 0 if iyear is or not a leap year)            
! ----------------------------------------------------------------------
!                                                                       
      if(iyear.lt.1900) then 
                print*,'*** Rocotlib/cp_leap_year: iyear= ',iyear 
                print*,'*** Rocotlib/cp_leap_year: iyear must be > 1900' 
                stop   '*** Rocotlib/cp_leap_year: iyear must be > 1900' 
                        endif 
!                                                                       
      ir=iyear-(iyear/4)*4 
      if(ir.eq.0) then 
                  ileap=1 
                  else 
                  ileap=0 
                  endif 
!                                                                       
      is=iyear-(iyear/100)*100 
      if(is.eq.0) then 
                  ir=iyear-(iyear/400)*400 
                  if(ir.eq.0) then 
                              ileap=1 
                              else 
                              ileap=0 
                              endif 
                  else 
                  return 
                  endif 
!                                                                       
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cp_seasons(iyear,id_sso,id_wso,id_seq,id_feq,          &
     &                            ct_sso,ct_wso,ct_seq,ct_feq)          
!                                                                       
      character*5 ct_sso,ct_wso,ct_seq,ct_feq 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : compute_seasons, i.e. solstice & equinox                 
! *   Author : P. Robert, SDev, 2017                                    
!                                                                       
! *   Input  : iyear (ex: 1980)                                         
! *   Output : id_sso,id_wso : June and december day of summer and      
!                              winter solstice.                         
!              id_seq,id_feq : same for march and september spring      
!                              and fall equinoxes.                      
!              ct_sso,ct_wso : June and december time of summer and     
!                              winter solstice.                         
!              ct_seq,ct_feq : same for march and september spring      
!                             and fall equinoxes; Cha*5 (hh:mm)         
! *                                                                     
! *   Comment: calcul des saisons: solstices & equinoxes,               
! *            precision 1 heure environ                                
! ----------------------------------------------------------------------
                                                                        
      id_sso=0 
      id_wso=0 
      id_seq=0 
      id_feq=0 
                                                                        
      ct_sso='?' 
      ct_wso='?' 
      ct_seq='?' 
      ct_feq='?' 
!                                                                       
! *** calcul de l'equinoxe de printemps                                 
!                                                                       
      imon=3 
      tetmax=0. 
                                                                        
      call cp_nbday_in_month(iyear,imon,nbday) 
                                                                        
      do iday=1,nbday 
      do ih=0,23 
      do im=0,50,10 
      is=0 
                                                                        
      call cp_time_param(iyear,imon,iday,ih,im,is) 
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
                                                                        
      teta=acos(szgeo)*180./3.14159 
                                                                        
      if(teta.gt.tetmax .and. szgeo.gt.0.) then 
                         tetmax=teta 
                         id_seq=iday 
                         write(ct_seq,"(i2.2,':',i2.2)") ih,im 
      endif 
                                                                        
      enddo 
      enddo 
      enddo 
!                                                                       
! *** calcul de l'equinoxe d'automne                                    
!                                                                       
      imon=9 
      tetmax=0. 
                                                                        
      call cp_nbday_in_month(iyear,imon,nbday) 
                                                                        
      do iday=1,nbday 
      do ih=0,23 
      do im=0,50,10 
      is=0 
                                                                        
      call cp_time_param(iyear,imon,iday,ih,im,is) 
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
                                                                        
      teta=acos(szgeo)*180./3.14159 
                                                                        
      if(teta.gt.tetmax .and. szgeo.gt.0.) then 
                         tetmax=teta 
                         id_feq=iday 
                         write(ct_feq,"(i2.2,':',i2.2)") ih,im 
      endif 
                                                                        
      enddo 
      enddo 
      enddo 
!                                                                       
! *** calcul du solstice d'ete                                          
!                                                                       
      imon=6 
      zmax=0. 
                                                                        
      call cp_nbday_in_month(iyear,imon,nbday) 
                                                                        
      do iday=1,nbday 
      do ih=0,23 
      do im=0,50,10 
      is=0 
                                                                        
      call cp_time_param(iyear,imon,iday,ih,im,is) 
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
                                                                        
      if(szgeo.gt.zmax) then 
                        zmax=szgeo 
                        id_sso=iday 
                        write(ct_sso,"(i2.2,':',i2.2)") ih,im 
      endif 
                                                                        
      enddo 
      enddo 
      enddo 
!                                                                       
! *** calcul du solstice d'hiver                                        
!                                                                       
      imon=12 
      zmin=1. 
                                                                        
      call cp_nbday_in_month(iyear,imon,nbday) 
                                                                        
      do iday=1,nbday 
      do ih=0,23 
      do im=0,50,10 
      is=0 
                                                                        
      call cp_time_param(iyear,imon,iday,ih,im,is) 
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
                                                                        
      if(szgeo.lt.zmin) then 
                        zmin=szgeo 
                        id_wso=iday 
                        write(ct_wso,"(i2.2,':',i2.2)") ih,im 
      endif 
                                                                        
      enddo 
      enddo 
      enddo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_doty_to_date(idoy,iyear,imonth,iday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_day_of_year for a given year in date             
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : iyear,idoy (idoy=1 for january 1)                        
! *   Output : imonth,iday                                              
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      dimension month(12) 
      data month/31,28,31,30,31,30,31,31,30,31,30,31/ 
      data ily /-1/ 
!                                                                       
      call cp_leap_year(iyear,ily) 
!                                                                       
      if(idoy.gt.365.and.ily.eq.0) then 
           print*,'*** Rocotlib/cv_doty_to_date: idoy= ',idoy 
           print*,'*** Rocotlib/cv_doty_to_date: iyear not a leap year' 
           print*,'                      no more than 365 day' 
           stop   '*** Rocotlib/cv_doty_to_date: iyear not a leap year' 
                                    endif 
!                                                                       
      if(idoy.lt.1) then 
        print*,'*** Rocotlib/cv_doty_to_date: idoy= ',idoy 
        print*,'*** Rocotlib/cv_doty_to_date: idoy must be > 0' 
        stop   '*** Rocotlib/cv_doty_to_date: idoy must be > 0' 
                     endif 
!                                                                       
      if(ily.eq.1) then 
                   month(2)=29 
                   else 
                   month(2)=28 
                   endif 
!                                                                       
      m=0 
!                                                                       
      do 10 im=1,12 
      mp=m 
      m=m+month(im) 
      if(idoy.le.m) go to 20 
   10 continue 
   20 continue 
!                                                                       
      imonth=im 
      iday=idoy-mp 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_jul2000_to_date(jd00,iyear,imonth,iday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_julian_day_2000 in date                          
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: compute date as year, month, day from julian day 2000    
!                                                                       
! *   Input  : jd00 julian day 2000 (0= 1/1/2000)                       
! *   Output : iyear,imonth,iday                                        
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      jd50 = jd00 +18262 
!                                                                       
      call cv_jul1950_to_date(jd50,iyear,imonth,iday) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_jul1950_to_date(jd50,iyear,imonth,iday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_julian_day_1950 in date                          
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: compute date as year, month, day from julian day 1950    
!                                                                       
! *   Input  : jd50  julian day 1950 (0= 1/1/1950)                      
! *   Output : iyear,imonth,iday                                        
! ----------------------------------------------------------------------
!                                                                       
      data ily /-1/ 
!                                                                       
!                                                                       
      jd1= -1 
!                                                                       
      do 10 iy=1950,3000 
      call cp_leap_year(iy,ily) 
      jdp=jd1 
      if(ily.eq.1) then 
                   jd1=jd1+366 
                   else 
                   jd1=jd1+365 
                   endif 
      if(jd1.ge.jd50) then 
                       iyear=iy 
                       go to 20 
                       endif 
!                                                                       
   10 continue 
      print*,'*** Rocotlib/cv_jul1950_to_date: jd50= ',jd50 
      print*,'*** Rocotlib/cv_jul1950_to_date: jd50 leads year > 3000' 
      stop   '*** Rocotlib/cv_jul1950_to_date: jd50 leads year > 3000' 
!                                                                       
   20 continue 
      jd=jd50-jdp 
      call cv_doty_to_date(jd,iy,imonth,iday) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_weekn_to_date(iweek,iyear,imonth,iday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_first_day_of_week_number in date                 
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : iweek,iyear                                              
! *   Output : imonth,iday                                              
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
!  calcul de la date correspondant au premier jour de la semaine        
!                                                                       
! *** premier lundi de l'annee                                          
!                                                                       
      do 10 ipl=1,7 
      call cv_date_to_dotw(iyear,1,ipl,idow) 
      if(idow.eq.1) go to 20 
   10 continue 
      print*, '*** Rocotlib/cv_weekn_to_date: error in idow computation' 
      stop    '*** Rocotlib/cv_weekn_to_date: error in idow computation' 
   20 continue 
!                                                                       
! *** semaine correspondant au premier lundi                            
!                                                                       
      call cv_date_to_weekn(iyear,1,ipl,iweek1) 
!                                                                       
! *** jour de l'annee correspondant au lundi de la semaine              
!                                                                       
      idoy= ipl +(iweek-iweek1)*7 
                                                                        
      if(idoy.lt.1) then 
                    iyear=iyear-1 
                    imonth=12 
                    iday=31+idoy 
                    return 
                    endif 
!                                                                       
! *** date correspondante                                               
!                                                                       
      call cv_doty_to_date(idoy,iyear,imonth,iday) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_date_to_dotw(iyear,imonth,iday,idow) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_date in day_of_the_week                          
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : iyear,imonth,iday                                        
! *   Output : idow (1-7)                                               
! ----------------------------------------------------------------------
!                                                                       
      data jdref,julday /-100000,-100000/ 
!                                                                       
!                                                                       
! *** jour julien du lundi de reference (le 1/1/2001 etait un lundi)    
!                                                                       
      call cv_date_to_jul2000(2001,1,1,jdref) 
!                                                                       
! *** jour julien de la date demandee                                   
!                                                                       
      call cv_date_to_jul2000(iyear,imonth,iday,julday) 
!                                                                       
! *** calcul du jour de la semaine                                      
!                                                                       
      idiff= julday-jdref 
!                                                                       
      iweek= idiff/7 
      idow= idiff -7*iweek +1 
      if(idow.lt.1) idow=idow+7 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_date_to_doty(iyear,imonth,iday,idoy) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_date in day_of_year with idoy=1 for january 1    
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : iyear,imonth,iday   ex: 1990,10,17                       
! *   Output : idoy                ex: 290                              
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      dimension month(12) 
      data month/31,28,31,30,31,30,31,31,30,31,30,31/ 
      data ily /-1/ 
!                                                                       
      call cp_leap_year(iyear,ily) 
!                                                                       
      if(ily.eq.1) then 
                   month(2)=29 
                   else 
                   month(2)=28 
                   endif 
!                                                                       
      if(imonth.lt.1.or.imonth.gt.12) then 
       print*, '*** Rocotlib/cv_date_to_doty: imonth = ',imonth 
       print*, '*** Rocotlib/cv_date_to_doty: imonth must be >0 and <13' 
       stop    '*** Rocotlib/cv_date_to_doty: imonth must be >0 and <13' 
                                      endif 
!                                                                       
      if(iday.gt.month(imonth)) then 
           print*, '*** Rocotlib/cv_date_to_doty: iday= ',iday 
           print*, '*** Rocotlib/cv_date_to_doty: this month has only', &
     &                                    month(imonth),'days'          
              stop '*** Rocotlib: error iday in cv_date_to_doty' 
                                endif 
!                                                                       
      idoy=iday 
      do 10 i=1,imonth-1 
      idoy=idoy+month(i) 
   10 continue 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_hms_to_dech(ih,im,is,houday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_hours_minutes_seconds in decimal hour of the day 
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : ih,im,is                                                 
! *   Output : houday  decimal hour of the day                          
! ----------------------------------------------------------------------
!                                                                       
      houday=float(ih)+float(im)/60.+float(is)/3600. 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_date_to_jul1950(iyear,imonth,iday,jd50) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_date in julian_day_1950 with jd50=0 for jan 1    
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : iyear,imonth,iday   ex: 1990,10,17                       
! *   Output : jd50                                                     
! ----------------------------------------------------------------------
!                                                                       
      data ily,idoy /-1,-1/ 
!                                                                       
!                                                                       
      if(iyear.lt.1950) then 
         print*, '*** Rocotlib/cv_date_to_jul1950: iyear= ',iyear 
         print*, '*** Rocotlib/cv_date_to_jul1950: iyear must be > 1950' 
         stop    '*** Rocotlib/cv_date_to_jul1950: iyear must be > 1950' 
                        endif 
!                                                                       
      call cv_date_to_doty(iyear,imonth,iday,idoy) 
!                                                                       
      jd50=idoy-1 
      do 10 i=1950,iyear-1 
      call cp_leap_year(i,ily) 
      if(ily.eq.1) then 
                   jd50=jd50+366 
                   else 
                   jd50=jd50+365 
                   endif 
   10 continue 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_date_to_jul2000(iyear,imonth,iday,jd00) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_date in julian_day_2000 with jd00=0 for january 1
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : iyear,imonth,iday   ex: 2001,10,17                       
! *   Output : jd00 (may be negative)                                   
! ----------------------------------------------------------------------
!                                                                       
      data jd50 /-1/ 
!                                                                       
!                                                                       
      if(iyear.lt.1950) then 
         print*, '*** Rocotlib/cv_date_to_jul2000: iyear= ',iyear 
         print*, '*** Rocotlib/cv_date_to_jul2000: iyear must be > 1950' 
         stop    '*** Rocotlib/cv_date_to_jul2000: iyear must be > 1950' 
                        endif 
!                                                                       
      call cv_date_to_jul1950(iyear,imonth,iday,jd50) 
!                                                                       
      jd00= jd50 - 18262 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_dhms_to_msotd(ih,im,is,ims,milday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_hours_minutes_seconds_ms in millisec_of_day      
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : ih,im,is,ims                                             
! *   Output : milday  millisecond of the day                           
! ----------------------------------------------------------------------
!                                                                       
      milday= ih*3600000 + im*60000 + is*1000 +ims 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_dech_to_hms(houday,ih,im,is) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_decimal hour of the day in time                  
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : houday decimal hour of the day                           
! *   Output : ih,im,is                                                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      nj=int(abs(houday)/24.) 
      rdech=abs(houday)-float(nj*24) 
!                                                                       
      ih=int(rdech) 
      im=int((rdech-float(ih))*60.) 
      is=int((rdech-float(ih)-float(im)/60.)*3600. +0.001) 
!                                                                       
      if(is.eq.60) then 
                   is=0 
                   im=im +1 
                   endif 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_msotd_to_hmsms(milday,ih,im,is,ims) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_millisec. of the day in time                     
! *   Author : P. Robert, CRPE, 2001                                    
!                                                                       
! *   Input  : milday millisecond of the day                            
! *   Output : ih,im,is,ims                                             
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      nj= milday/86400000 
      mil2= milday-nj*86400000 
!                                                                       
      if(nj.ne.0) then 
          print*, '*** Rocotlib/cv_msotd_to_hmsms: milday=',milday,     &
     &                                           'is > 86400000'        
          print*, '                      assumed:',mil2 
                  endif 
!                                                                       
      ih=  milday/3600000 
      im= (milday-ih*3600000)/60000 
      is= (milday-ih*3600000-im*60000)/1000 
      ims= milday-ih*3600000 -im*60000 -is*1000 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cv_date_to_weekn(iyear,imonth,iday,iweek) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : calendar modules of Rocotlib Software                    
! *   Object : convert_date in week_of_the_year                         
! *   Author : P. Robert, CRPE, 2001, revised SD April 2017             
!                                                                       
! *   Input  : iyear,imonth,iday                                        
! *   Output : iweek                                                    
! ----------------------------------------------------------------------
!                                                                       
!     la semaine 1 de toute année est celle qui contient le 4 janvier   
!     ou celle qui contient le 1er jeudi de janvier(norme ISO-8601).    
                                                                        
! *** calcul du jour de la semaine du 1er janvier de l'annee            
                                                                        
      call cv_date_to_dotw(iyear,1,1,idow) 
                                                                        
! *** calcul du lundi commencant la semaine 2                           
                                                                        
      if(idow.eq.1) idoty2= 8 
      if(idow.eq.2) idoty2= 7 
      if(idow.eq.3) idoty2= 6 
      if(idow.eq.4) idoty2= 5 
      if(idow.eq.5) idoty2=11 
      if(idow.eq.6) idoty2=10 
      if(idow.eq.7) idoty2= 9 
!                                                                       
! *** jour de l'annee de la date demandee                               
!                                                                       
      call cv_date_to_doty(iyear,imonth,iday,idoty) 
!                                                                       
! *** jour de la semaine de la date demandee                            
!                                                                       
      call cv_date_to_dotw(iyear,imonth,iday,idow) 
!                                                                       
! *** jour de l'annee du lundi de la semaine de la date demandee        
!                                                                       
      idotL= idoty -idow +1 
                                                                        
! *** numero de la semaine                                              
                                                                        
      iweek= (idotL -idoty2)/7 +2 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine g_gei_geo_dipole_dir(dxgei,dygei,dzgei,                &
     &                                dxgeo,dygeo,dzgeo)                
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : give modules of Rocotlib Software                        
! *   Object : give_dipole_direction in GEI and GEO system              
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: values are extracted from common                         
!                                                                       
! *   Input  : none                                                     
! *   Output : dxgei,dygei,dzgei  cartesian dipole GEI coordinates      
!              dxgeo,dygeo,dzgeo  cartesian dipole GEO coordinates      
! ----------------------------------------------------------------------
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
!                                                                       
      dxgei=gm1 
      dygei=gm2 
      dzgei=gm3 
!                                                                       
      dxgeo=gd1 
      dygeo=gd2 
      dzgeo=gd3 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine g_gsm_dipole_tilt_angle(diptan) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : give modules of Rocotlib Software                        
! *   Object : give_dipole_tilt_angle in radians                        
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: values are extracted from common                         
!                                                                       
! *   Input  : none                                                     
! *   Output : diptan (radians)                                         
! ----------------------------------------------------------------------
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      diptan=atan2(smu,cmu) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine g_gei_geo_ecliptic_dir(exgei,eygei,ezgei,              &
     &                                  exgeo,eygeo,ezgeo)              
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : give modules of Rocotlib Software                        
! *   Object : give_ecliptic_direction in GEI and GEO system            
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: values are extracted from common                         
!                                                                       
! *   Input  : none                                                     
! *   Output : exgei,eygei,ezgei  cartesian ecliptic GEI coordinates    
!              exgeo,eygeo,ezgeo  cartesian ecliptic GEO coordinates    
! ----------------------------------------------------------------------
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3,pr1,pr2,pr3 
!                                                                       
      exgei=ge1 
      eygei=ge2 
      ezgei=ge3 
!                                                                       
      exgeo=pe1 
      eygeo=pe2 
      ezgeo=pe3 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine g_gei_geo_sun_rot(rxgei,rygei,rzgei,rxgeo,rygeo,rzgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : give modules of Rocotlib Software                        
! *   Object : give_sun_rotation_direction in GEI and GEO system        
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: values are extracted from common                         
!                                                                       
! *   Input  : none                                                     
! *   Output : rxgei,rygei,rzgei cartesian sun rotation GEI coordinates 
!              rxgeo,rygeo,rzgeo cartesian sun rotation GEO coordinates 
! ----------------------------------------------------------------------
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3,pr1,pr2,pr3 
!                                                                       
      rxgei=gr1 
      rygei=gr2 
      rzgei=gr3 
!                                                                       
      rxgeo=pr1 
      rygeo=pr2 
      rzgeo=pr3 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : give modules of Rocotlib Software                        
! *   Object : give_sun_direction in GEI and GEO system                 
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: values are extracted from common                         
!                                                                       
! *   Input  : none                                                     
! *   Output : sxgei,sygei,szgei  cartesian sun GEI coordinates         
!              sxgeo,sygeo,szgeo  cartesian sun GEO coordinates         
! ----------------------------------------------------------------------
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3,pr1,pr2,pr3 
!                                                                       
      sxgei=gs1 
      sygei=gs2 
      szgei=gs3 
!                                                                       
      sxgeo=ps1 
      sygeo=ps2 
      szgeo=ps3 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine g_gei_sun_param(gmst,slon,sras,sdec,obli) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : give modules of Rocotlib Software                        
! *   Object : give_sun_parameter dependant of time in GEI system       
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: values are extracted from common                         
!                                                                       
! *   Input  : none                                                     
! *   Output : gmst   greenwich mean sideral time (radians)             
!              slon   longitude along ecliptic (radians)                
!              sras   right ascension (radians)                         
!              sdec   declination of the sun (radians)                  
! ----------------------------------------------------------------------
!                                                                       
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq 
!                                                                       
      gmst= gst 
      slon= slong 
      sras= srasn 
      sdec= sdecl 
      obli= obliq 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine g_rocot_version_number(vernum,verdat) 
!                                                                       
      character*14 verdat 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : give modules of Rocotlib Software                        
! *   Object : give_version_number and modification date of the library 
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: values are extracted from common                         
!                                                                       
! *   Input  : none                                                     
! *   Output : vernum (ex 1.0) and verdat (ex: 'January 1995')          
! ----------------------------------------------------------------------
!                                                                       
!     For previous versions informations, see print_rocot_info subrout. 
!                                                                       
      vernum=3.2 
      verdat='Fev 2020' 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine print_rocot_info 
!                                                                       
      character*14 verdat 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : print modules of Rocotlib Software                       
! *   Object : print_library_informations                               
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: could be a main program                                  
!                                                                       
! *   Input  : none                                                     
! *   Output : none; print infos on output                              
! ----------------------------------------------------------------------
!                                                                       
      call g_rocot_version_number(vernum,verdat) 
!                                                                       
!                                                                       
      print 100, ' ' 
      print 100, '*****************************************************' 
      print 100, ' ' 
      print 100, '    Coordinates Transformation Library  ROCOTLIB' 
      print 100, ' ' 
      print 200, '        Revised Version ',vernum, ' - ',verdat 
      print 100, ' ' 
      print 100, '               ___________________' 
      print 100, ' ' 
      print 100, '              initially supported by ' 
      print 100, ' ' 
      print 100, '              EUROPEAN  SPACE  AGENCY' 
      print 100, ' ' 
      print 100, '           Study of the Cluster Mission' 
      print 100, '             Planning Related Aspects' 
      print 100, '      within the Numerical Simulations Network' 
      print 100, ' ' 
      print 100, '        Patrick ROBERT, CRPE, November 1992' 
      print 100, '               ___________________' 
      print 100, ' ' 
      print 100, ' version 1.0, November 1992                          ' 
      print 100, ' version 1.1, July     1993                          ' 
      print 100, ' version 1.2, January  1995                          ' 
      print 100, ' version 1.3, July     2000 (Jul. day 2000/sun dir)  ' 
      print 100, ' version 1.4, June     2001 (for automatic docum.)   ' 
      print 100, ' version 1.5, December 2001 (add cp_sunrise_sunset   ' 
      print 100, ' version 1.6, Juin     2002 (upgrade IGRF -> 2005)   ' 
      print 100, ' version 1.7, December 2002 (Version for CDPP)       ' 
      print 100, ' version 1.8, November 2003 (add t_sr2_to_sr)        ' 
      print 100, ' version 1.9, March    2004 (compatibility with IDL) ' 
      print 100, ' version 2.0, November 2006 (Update  IGRF -> 2010)   ' 
      print 100, ' version 2.1, November 2006 (Update  IGRF -> 2015)   ' 
      print 100, ' version 2.2, December 2011 (cp_sunrise_sunset polar)' 
      print 100, ' version 3.0, May      2017 (IGRF->table geomag. pole' 
      print 100, '                     -> 2020 +some other coord. sys.)' 
      print 100, ' version 3.1, January  2019 (compatibility with V2.2)' 
      print 100, ' version 3.2, February 2020 (add trans. Euler angles)' 
      print 100, ' ' 
      print 100, ' ' 
      print 100, '      Copyright 1992, Patrick ROBERT, CNRS-ESA,' 
      print 100, '               All Rights reserved' 
      print 100, '               ___________________' 
      print 100, ' ' 
      print 100, '    For details, see the orginal document untitled:' 
      print 100, ' ' 
      print 100, '              CLUSTER Software Tools' 
      print 100, '      Part I: Coordinate Transformation Library' 
      print 100, '          Document de travail DT/CRPE/1231' 
      print 100, '        Patrick Robert, CRPE/TID, Juillet 1993' 
      print 100, ' ' 
      print 100, '    Available at CDPP:' 
      print 100, ' ' 
      print 100, '    ROCOTLIB: a coordinate Transformation Library' 
      print 100, '              for Solar-Terrestrial studies' 
      print 100, '    Patrick ROBERT, version 1.7 - January 2003,' 
      print 100, '        Rapport Interne no RI-CETP/02/2003' 
      print 100, '               ___________________' 
      print 100, ' ' 
      print 100, ' The present version 3.0 contains new transformations' 
      print 100, ' (TPN, MVA), and the cp_sunrise_sunset sub. allowing ' 
      print 100, ' computation of sunrise and sunset anywhere on Earth,' 
      print 100, ' including polar zones.          ' 
      print 100, ' ' 
      print 100, ' The determination of the dipole axis direction' 
      print 100, ' has been replaced by a table available from years' 
      print 100, ' 1900 to 2020, in order to have an easy way to code ' 
      print 100, ' maintenance in the coming years.' 
      print 100, ' ' 
      print 100, ' Operations on matrix (somme, product, inversion, ' 
      print 100, ' computation of eigen vectors, diagonalization...)' 
      print 100, ' required for Minimum Variance Analysis coordinates' 
      print 100, ' have been added with prefix "mat_".' 
      print 100, ' '
      print 100, ' The previous and temporary V3.1 is the same as 3.0 '
      print 100, ' but is compatible with V2.2 version concerning '
      print 100, ' subroutine name (ex: tmaggsm -> t_mag_to_gsm) '
      print 100, ' '
      print 100, ' The present V3.2 version contains transformations '
      print 100, ' with Euler angles and their interpolation. '
      print 100, ' ' 
      print 100, ' Original code developped in Fortran 77.' 
      print 100, ' Fortran 90, C, and IDL codes are also available.' 
      print 100, ' ' 
      print 100, '*****************************************************' 
      print 100, ' ' 
!                                                                       
  100 format(3a) 
  200 format(a,f4.1,2a) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine r_coordinate_values(x,y,z,cs) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : read modules of Rocotlib Software                        
! *   Object : read coordinate values from input                        
! *   Author : P. Robert, CRPE, 2002                                    
!                                                                       
! *   Comment: read cs and x,y,z cartesian  or spherical coordinates    
! *            print error if cs is not valid, and ask again            
!                                                                       
! *   Input  : cs (c or s) and x,y,z on standard input                  
! *   Output : cs and x,y,z always in cartesian coordinates             
! ----------------------------------------------------------------------
!                                                                       
      character*1 cs 
!                                                                       
!                                                                       
      pisd=acos(-1.)/180. 
!                                                                       
   10 continue 
!                                                                       
      print 100, 'input coordinates: cartesian (c) or spherical (s)' 
      read *, cs 
      print 100, cs 
!                                                                       
      if(cs.ne.'s'.and.cs.ne.'c') then 
                                  print*, 'only c or s' 
                                  print*, 'again...' 
                                  go to 10 
                                  endif 
      if(cs.eq.'s') then 
                    print 100, 'R,tet, phi (RE,deg.) ?' 
                    read *,  r,tet, phi 
                    print 110,  r,tet, phi 
!                                                                       
                    tet=tet*pisd 
                    phi=phi*pisd 
                    call t_sph_to_car(r,tet,phi,x,y,z) 
                    print 120, 'then x,y,z = ',x,y,z,'(RE)' 
!                                                                       
                    else 
                    print*, 'x,y,z ? (RE) ' 
                    read *,  x,y,z 
                    print*,  x,y,z 
!                                                                       
                    call t_car_to_sph(x,y,z,r,tet,phi) 
                    tet=tet/pisd 
                    phi=phi/pisd 
                    print 200, 'then R,tet,phi= ',r,tet,phi,' (RE,deg.)' 
                    endif 
!                                                                       
  100 format(a) 
  110 format(1Pe14.6,1x,f8.3,1x,f8.3) 
  120 format(1x,a,3(1Pe11.3),a) 
  200 format(1x,a,1Pe11.3,f8.3,f8.3,a) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine r_coordinate_system(csys) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : read modules of Rocotlib Software                        
! *   Object : read coordinate system from input and check validity     
! *   Author : P. Robert, CRPE, 2002                                    
!                                                                       
! *   Comment: read csys string variable and check validity             
! *            (only gei, geo, mag, sma, gsm, gse, gsq)                 
! *            print error if csys is not valid, and ask again          
!                                                                       
! *   Input  : csys on standard input                                   
! *   Output : csys                                                     
! ----------------------------------------------------------------------
!                                                                       
      character*3 csys 
!                                                                       
   10 continue 
!                                                                       
      print 100, 'repere ? (gei, geo, mag, sma, gsm, gse, gsq)' 
      read *,csys 
      print 100,csys 
  100 format(a) 
!                                                                       
      if(csys.eq.'gei') return 
      if(csys.eq.'geo') return 
      if(csys.eq.'mag') return 
      if(csys.eq.'sma') return 
      if(csys.eq.'gsm') return 
      if(csys.eq.'gse') return 
      if(csys.eq.'gsq') return 
!                                                                       
      print*, 'Only gei, geo, mag, sma, gsm, gse or gsq please...' 
      print*, 'again...' 
      go to 10 
!                                                                       
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine r_date(iyear,imonth,iday) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : read modules of Rocotlib Software                        
! *   Object : read_date from input and check validity                  
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Comment: test if year is gt 1900                                  
! *            test if imonth is not greater than 12                    
!              test if iday is not greather than lengh of the month,    
!                        takink into account the leap years.            
!               print error if date is not valid, and ask again         
!                                                                       
! *   Input  : iyear,imonth,iday given in standard input                
! *   Output : iyear,imonth,iday                                        
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      dimension month(12) 
      data month/31,28,31,30,31,30,31,31,30,31,30,31/ 
      data ily /-1/ 
!                                                                       
   10 continue 
      print 100, ' iyear,imonth,iday ? (ex: 1990,10,17)' 
      read *,  iyear,imonth,iday 
      print 110,  iyear,imonth,iday 
!                                                                       
  100 format(a) 
  110 format(1x,i4,1x,i2,1x,i2) 
!                                                                       
      if(iyear.lt.1900) then 
         print*, '*** Rocotlib/r_date: iyear must be greater than 1900' 
         print*, '                     again ...' 
                        go to 10 
                        endif 
!                                                                       
      call cp_leap_year(iyear,ily) 
!                                                                       
!                                                                       
      if(ily.eq.1) then 
                   month(2)=29 
                   else 
                   month(2)=28 
                   endif 
!                                                                       
      if(imonth.lt.1.or.imonth.gt.12) then 
          print*, '*** Rocotlib/r_date: imonth must be between 1 or 12 ' 
          print*, '                     again...' 
                     go to 10 
                                      endif 
!                                                                       
      if(iday.gt.month(imonth)) then 
                     print*, '*** Rocotlib/r_date: this month has only',&
     &                            month(imonth),'days'                  
                     print*, '                     again...' 
                go to 10 
                                endif 
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine r_time(ih,im,is) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : read modules of Rocotlib Software                        
! *   Object : read_time from input and check validity                  
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: read hour, minute and second and verifie validity        
!              ih must be between 1 and 23, im and is between 1 and 59  
!              print error if time is not valid, and ask again          
!                                                                       
! *   Input  : ih,im,is on standard input                               
! *   Output : ih,im,is                                                 
! ----------------------------------------------------------------------
!                                                                       
   10 continue 
      print 100, ' hour, minute, second ? (ex: 10,45,50)' 
      read *, ih,im,is 
      print 110, ih,im,is 
!                                                                       
  100 format(a) 
  110 format(1x,i2,1x,i2,1x,i2) 
!                                                                       
      if(ih.lt.0.or.ih.gt.23) then 
                              print*, 'hour between 0 and 23 please' 
                              print*, 'again...' 
                              go to 10 
                              endif 
!                                                                       
      if(im.lt.0.or.im.gt.59) then 
                              print*, 'minute between 0 and 59 please' 
                              print*, 'again...' 
                              go to 10 
                              endif 
!                                                                       
      if(is.lt.0.or.is.gt.59) then 
                              print*, 'second between 0 and 59 please' 
                              print*, 'again...' 
                              go to 10 
                              endif 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_car_to_sph(x,y,z,r,teta,phi) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_car_to_sph: CAR -> SPH  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: none                                                     
!                                                                       
! *   Input  :   x,y,z        cartesian coordinates                     
! *   Output :   r,teta,phi   spherical coordinates (angles in radians) 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      teta=0. 
      phi=0. 
      sq=x**2+y**2 
      r=sqrt(sq+z**2) 
      pi=acos(-1.)
      pisd=pi/180.
!                                                                       
      if(r.lt.1.e-30) return 
!                                                                       
! *** en dessous de 1/10000 degres, on considere que teta =0 ou 180     
!     et phi indefini, mis a zero                                       
!                                                                       
      if(sq/r.gt.1.7e-6) then 
                   phi=atan2(y,x) 
                 if(phi.lt.-179.999*pisd.and.phi.gt.-180.*pisd) phi=pi
                 if(phi.lt. 0.0001*pisd.and.phi.gt.-0.0001*pisd) phi=0.
                   teta=acos(z/r) 
                   return 
                   endif 
!                                                                       
      if (z.lt.0.) then 
                   teta=acos(-1.) 
                   else 
                   teta=0. 
                   endif 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_dm_to_geo(xdme,ydme,zdme,rlat,rlong,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_dme_to_geo: DM  -> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xdme,ydme,zdme   cartesian dme coordinates               
!              rlat,rlong       latitude and longitude of the point     
!                               of observation (radians)                
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
!                                                                       
      q=cos(rlat) 
      r1=q*cos(rlong) 
      r2=q*sin(rlong) 
      r3=sin(rlat) 
!                                                                       
      y1= real(dprod(gd2,r3) - dprod(gd3,r2)) 
      y2= real(dprod(gd3,r1) - dprod(gd1,r3)) 
      y3= real(dprod(gd1,r2) - dprod(gd2,r1)) 
!                                                                       
      q=sqrt(y1**2 + y2**2 + y3**2) 
!                                                                       
      y1=y1/q 
      y2=y2/q 
      y3=y3/q 
!                                                                       
      x1= real(dprod(y2,gd3) - dprod(y3,gd2)) 
      x2= real(dprod(y3,gd1) - dprod(y1,gd3)) 
      x3= real(dprod(y1,gd2) - dprod(y2,gd1)) 
!                                                                       
!                                                                       
      xgeo= x1*xdme + y1*ydme + gd1*zdme 
      ygeo= x2*xdme + y2*ydme + gd2*zdme 
      zgeo= x3*xdme + y3*ydme + gd3*zdme 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gei_to_geo(xgei,ygei,zgei,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gei_to_geo: GEI -> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgei,ygei,zgei   cartesian gei coordinates               
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq 
!                                                                       
      xgeo=  cgst*xgei + sgst*ygei 
      ygeo= -sgst*xgei + cgst*ygei 
      zgeo=  zgei 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gei_to_gse(xgei,ygei,zgei,xgse,ygse,zgse) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gei_to_gse: GEI -> GSE  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgei,ygei,zgei cartesian gei coordinates                 
! *   Output : xgse,ygse,zgse cartesian gse coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3 
!                                                                       
      xgse=   gs1*xgei +   gs2*ygei +   gs3*zgei 
      ygse= gegs1*xgei + gegs2*ygei + gegs3*zgei 
      zgse=   ge1*xgei +   ge2*ygei +   ge3*zgei 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gei_to_gsm(xgei,ygei,zgei,xgsm,ygsm,zgsm) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gei_to_gsm: GEI -> GSM  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgei,ygei,zgei   cartesian gei coordinates               
! *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3 
!                                                                       
      xgsm=    gs1*xgei +    gs2*ygei +    gs3*zgei 
      ygsm= yeigm1*xgei + yeigm2*ygei + yeigm3*zgei 
      zgsm= zeigm1*xgei + zeigm2*ygei + zeigm3*zgei 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gei_to_gseq(xgei,ygei,zgei,xgsq,ygsq,zgsq) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gei_to_gsq: GEI -> GSEQ system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgei,ygei,zgei   cartesian gei coordinates               
! *   Output : xgsq,ygsq,zgsq   cartesian gsq coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3 
!                                                                       
      xgsq=    gs1*xgei +    gs2*ygei +    gs3*zgei 
      ygsq= yeigq1*xgei + yeigq2*ygei + yeigq3*zgei 
      zgsq= zeigq1*xgei + zeigq2*ygei + zeigq3*zgei 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gei_to_mag(xgei,ygei,zgei,xmag,ymag,zmag) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gei_to_mag: GEI -> MAG  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgei,ygei,zgei   cartesian gei coordinates               
! *   Output : xmag,ymag,zmag   cartesian mag coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3 
!                                                                       
      xmag= xeima1*xgei + xeima2*ygei + xeima3*zgei 
      ymag= yeima1*xgei + yeima2*ygei + yeima3*zgei 
      zmag=    gm1*xgei +    gm2*ygei +    gm3*zgei 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gei_to_sm(xgei,ygei,zgei,xsma,ysma,zsma) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gei_to_sma: GEI -> SM   system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgei,ygei,zgei   cartesian gei coordinates               
! *   Output : xsma,ysma,zsma   cartesian sma coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3 
!                                                                       
      xsma= xeism1*xgei + xeism2*ygei + xeism3*zgei 
      ysma= yeism1*xgei + yeism2*ygei + yeism3*zgei 
      zsma=    gm1*xgei +    gm2*ygei +    gm3*zgei 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_dm(xgeo,ygeo,zgeo,rlat,rlong,xdme,ydme,zdme) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_dme: GEO -> DM   system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates               
!              rlat,rlong       latitude and longitude of the point     
!                               of observation (radians)                
! *   Output : xdme,ydme,zdme   cartesian dme coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
!                                                                       
      q=cos(rlat) 
      r1=q*cos(rlong) 
      r2=q*sin(rlong) 
      r3=sin(rlat) 
!                                                                       
      y1= real(dprod(gd2,r3) - dprod(gd3,r2)) 
      y2= real(dprod(gd3,r1) - dprod(gd1,r3)) 
      y3= real(dprod(gd1,r2) - dprod(gd2,r1)) 
!                                                                       
      q=sqrt(y1**2 + y2**2 + y3**2) 
!                                                                       
      y1=y1/q 
      y2=y2/q 
      y3=y3/q 
!                                                                       
      x1= real(dprod(y2,gd3) - dprod(y3,gd2)) 
      x2= real(dprod(y3,gd1) - dprod(y1,gd3)) 
      x3= real(dprod(y1,gd2) - dprod(y2,gd1)) 
!                                                                       
      xdme=  x1*xgeo +  x2*ygeo +  x3*zgeo 
      ydme=  y1*xgeo +  y2*ygeo +  y3*zgeo 
      zdme= gd1*xgeo + gd2*ygeo + gd3*zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_gei(xgeo,ygeo,zgeo,xgei,ygei,zgei) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_gei: GEO -> GEI  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgeo,ygeo,zgeo cartesian geo coordinates                 
! *   Output : xgei,ygei,zgei cartesian gei coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq 
!                                                                       
      xgei= cgst*xgeo - sgst*ygeo 
      ygei= sgst*xgeo + cgst*ygeo 
      zgei= zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_gse(xgeo,ygeo,zgeo,xgse,ygse,zgse) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_gse: GEO -> GSE  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates               
! *   Output : xgse,ygse,zgse   cartesian gse coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3 
!                                                                       
      xgse=   ps1*xgeo +   ps2*ygeo +   ps3*zgeo 
      ygse= peps1*xgeo + peps2*ygeo + peps3*zgeo 
      zgse=   pe1*xgeo +   pe2*ygeo +   pe3*zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_gsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_gsm: GEO -> GSM  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates               
! *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3 
!                                                                       
      xgsm=    ps1*xgeo +    ps2*ygeo +    ps3*zgeo 
      ygsm= yeogm1*xgeo + yeogm2*ygeo + yeogm3*zgeo 
      zgsm= zeogm1*xgeo + zeogm2*ygeo + zeogm3*zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_gseq(xgeo,ygeo,zgeo,xgsq,ygsq,zgsq) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_gsq: GEO -> GSEQ system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates               
! *   Output : xgsq,ygsq,zgsq   cartesian gsq coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3 
!                                                                       
      xgsq=    ps1*xgeo +    ps2*ygeo +    ps3*zgeo 
      ygsq= yeogq1*xgeo + yeogq2*ygeo + yeogq3*zgeo 
      zgsq= zeogq1*xgeo + zeogq2*ygeo + zeogq3*zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_mag(xgeo,ygeo,zgeo,xmag,ymag,zmag) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_mag: GEO -> MAG  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates               
! *   Output : xmag,ymag,zmag   cartesian mag coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3 
!                                                                       
      xmag= xeoma1*xgeo + xeoma2*ygeo + xeoma3*zgeo 
      ymag= yeoma1*xgeo + yeoma2*ygeo + yeoma3*zgeo 
      zmag=    gd1*xgeo +    gd2*ygeo +    gd3*zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_sm(xgeo,ygeo,zgeo,xsma,ysma,zsma) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_sma: GEO -> SM   system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates               
! *   Output : xsma,ysma,zsma   cartesian sma coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3 
!                                                                       
      xsma= xeosm1*xgeo + xeosm2*ygeo + xeosm3*zgeo 
      ysma= yeosm1*xgeo + yeosm2*ygeo + yeosm3*zgeo 
      zsma=    gd1*xgeo +    gd2*ygeo +    gd3*zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_geo_to_vdh(xgeo,ygeo,zgeo,rlat,rlong,xvdh,yvdh,zvdh) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_geo_to_vdh: GEO -> VDH  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates               
!              rlat,rlong       latitude and longitude of the point     
!                               of observation (radians)                
! *   Output : xvdh,yvdh,zvdh   cartesian vdh coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      q=cos(rlat) 
      r1=q*cos(rlong) 
      r2=q*sin(rlong) 
      r3=sin(rlat) 
!                                                                       
      v1=r1 
      v2=r2 
      v3=r3 
!                                                                       
      q12=sqrt(r1**2 + r2**2) 
!                                                                       
      d1= -r2/q12 
      d2=  r1/q12 
      d3=  0. 
!                                                                       
      h1= -r1*r3/q12 
      h2= -r2*r3/q12 
      h3=  q12 
!                                                                       
      xvdh= v1*xgeo + v2*ygeo + v3*zgeo 
      yvdh= d1*xgeo + d2*ygeo + d3*zgeo 
      zvdh= h1*xgeo + h2*ygeo + h3*zgeo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gse_to_gei(xgse,ygse,zgse,xgei,ygei,zgei) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gse_to_gei: GSE -> GEI  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgse,ygse,zgse    cartesian gse coordinates              
! *   Output : xgei,ygei,zgei    cartesian gei coordinates              
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3 
!                                                                       
      xgei= gs1*xgse + gegs1*ygse + ge1*zgse 
      ygei= gs2*xgse + gegs2*ygse + ge2*zgse 
      zgei= gs3*xgse + gegs3*ygse + ge3*zgse 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gse_to_geo(xgse,ygse,zgse,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gse_to_geo: GSE -> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgse,ygse,zgse   cartesian gse coordinates               
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3 
!                                                                       
      xgeo= ps1*xgse + peps1*ygse + pe1*zgse 
      ygeo= ps2*xgse + peps2*ygse + pe2*zgse 
      zgeo= ps3*xgse + peps3*ygse + pe3*zgse 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gse_to_gsm(xgse,ygse,zgse,xgsm,ygsm,zgsm) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gse_to_gsm: GSE -> GSM  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgse,ygse,zgse   cartesian gse coordinates               
! *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xgsm= xgse 
      ygsm=  cdze*ygse + sdze*zgse 
      zgsm= -sdze*ygse + cdze*zgse 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gse_to_gseq(xgse,ygse,zgse,xgsq,ygsq,zgsq) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gse_to_gsq: GSE -> GSEQ system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgse,ygse,zgse   cartesian gse coordinates               
! *   Output : xgsq,ygsq,zgsq   cartesian gsq coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xgsq= xgse 
      ygsq= ctetq*ygse - stetq*zgse 
      zgsq= stetq*ygse + ctetq*zgse 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gse_to_mfa(xgse,ygse,zgse,bx,by,bz,xmfa,ymfa,zmfa) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gse_to_mfa: GSE -> MFA  system                
! *   Author : P. Robert, LPP , 2016                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xgse,ygse,gsez cartesian gse coordinates                 
!              bx,  by,  bz   cartesian gse coordinates of DC mag field 
! *   Output : xmfa,ymfa,zmfa cartesian mfa coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
! *** tranform the vector from gse to mfa                               
!                                                                       
      bperp= sqrt(bx*bx + by*by) 
      b0= sqrt(bx*bx + by*by + bz*bz) 
!                                                                       
! **  first rotation                                                    
!                                                                       
      sinphi=by/bperp 
      cosphi=bx/bperp 
!                                                                       
      xp=  cosphi*xgse + sinphi*ygse 
      yp= -sinphi*xgse + cosphi*ygse 
      zp=  zgse 
!                                                                       
! **  second rotation                                                   
!                                                                       
      sintet= bperp/b0 
      costet= bz/b0 
!                                                                       
      xmag= costet*xp - sintet*zp 
      ymag= yp 
      zmag= sintet*xp + costet*zp 
!                                                                       
! **  third rotation                                                    
!                                                                       
      sxm= costet*cosphi 
      sym=-sinphi 
!                                                                       
      smperp= sqrt(sxm*sxm + sym*sym) 
!                                                                       
      sinpsi= sym/smperp 
      cospsi= sxm/smperp 
!                                                                       
      xmfa=  cospsi*xmag + sinpsi*ymag 
      ymfa= -sinpsi*xmag + cospsi*ymag 
      zmfa=  zmag 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gse_to_sr2(xgse,ygse,zgse,rotx,roty,rotz,            &
     &                        xsr2,ysr2,zsr2)                           
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gse_to_sr2: GSE -> SR2  system                
! *   Author : P. Robert, CETP, 2001                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xgse,ygse,zgse cartesian gse coordinates                 
!              rotx,roty,rotz cartesian gse coordinates of rotation axis
! *   Output : xsr2,ysr2,zsr2 cartesian sr2 coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
! *** set transform matrix with Spin axis terms                         
!                                                                       
      rmod= sqrt(rotx**2 + roty**2 + rotz**2) 
!                                                                       
      rx=rotx/rmod 
      ry=roty/rmod 
      rz=rotz/rmod 
!                                                                       
      a= 1./sqrt(ry*ry + rz*rz) 
!                                                                       
      x1= (ry*ry + rz*rz)*a 
      x2= -rx*ry*a 
      x3= -rx*rz*a 
!                                                                       
      y1=  0. 
      y2=  rz*a 
      y3= -ry*a 
!                                                                       
      z1= rx 
      z2= ry 
      z3= rz 
!                                                                       
! *** tranform the input vector from gse to sr2                         
!                                                                       
      xsr2= x1*xgse + x2*ygse + x3*zgse 
      ysr2= y1*xgse + y2*ygse + y3*zgse 
      zsr2= z1*xgse + z2*ygse + z3*zgse 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gse_to_tpn(xgse,ygse,zgse,xo,yo,zo,xs,xtpn,ytpn,ztpn) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gse_to_tpn: GSE -> TPN  system                
! *   Author : P. Robert, LPP , 2016                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xgse,ygse,zgse cartesian gse vector                      
!              xo,yo,zo position of the S/C in gse                      
!              xs subsolar point, submit of the paraboloid              
!              from Earth to Sun                                        
! *   Output : xtpn,ytpn,ztpn cartesian tpn coordinates                 
! ----------------------------------------------------------------------
!                                                                       
      real Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz 
!                                                                       
!     computation of the TPN axis in gse system                         
!                                                                       
      call cp_tpn_param(xo,yo,zo,xs, Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz) 
!                                                                       
      xtpn= Tx*xgse + Ty*ygse +Tz*zgse 
      ytpn= Px*xgse + Py*ygse +Pz*zgse 
      ztpn= Nx*xgse + Ny*ygse +Nz*zgse 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gsm_to_gei(xgsm,ygsm,zgsm,xgei,ygei,zgei) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsm_to_gei: GSM -> GEI  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! *   Output : xgei,ygei,zgei   cartesian gei coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3 
!                                                                       
      xgei= gs1*xgsm + yeigm1*ygsm + zeigm1*zgsm 
      ygei= gs2*xgsm + yeigm2*ygsm + zeigm2*zgsm 
      zgei= gs3*xgsm + yeigm3*ygsm + zeigm3*zgsm 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gsm_to_geo(xgsm,ygsm,zgsm,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsm_to_geo: GSM -> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3 
!                                                                       
      xgeo= ps1*xgsm + yeogm1*ygsm + zeogm1*zgsm 
      ygeo= ps2*xgsm + yeogm2*ygsm + zeogm2*zgsm 
      zgeo= ps3*xgsm + yeogm3*ygsm + zeogm3*zgsm 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gsm_to_gse(xgsm,ygsm,zgsm,xgse,ygse,zgse) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsm_to_gse: GSM -> GSE  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! *   Output : xgse,ygse,zgse   cartesian gse coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xgse= xgsm 
      ygse= cdze*ygsm - sdze*zgsm 
      zgse= sdze*ygsm + cdze*zgsm 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gsm_to_gseq(xgsm,ygsm,zgsm,xgsq,ygsq,zgsq) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsm_to_gsq: GSM -> GSQ  system                
! *   Author : P. Robert, CRPE, 2002                                    
!                                                                       
! *   Input  : xgsm,ygsm,zgsm cartesian gsm coordinates                 
! *   Output : xgsq,ygsq,zgsq cartesian gsq coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      call t_gsm_to_gse(xgsm,ygsm,zgsm,ax,ay,az) 
      call t_gse_to_gseq(ax,ay,az,xgsq,ygsq,zgsq) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gsm_to_mag(xgsm,ygsm,zgsm,xmag,ymag,zmag) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsm_to_mag: GSM -> MAG  system                
! *   Author : P. Robert, CRPE, 2002                                    
!                                                                       
! *   Input  : xgsm,ygsm,zgsm cartesian gsm coordinates                 
! *   Output : xmag,ymag,zmag cartesian mag coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      call t_gsm_to_sm(xgsm,ygsm,zgsm,ax,ay,az) 
      call t_sm_to_mag(ax,ay,az,xmag,ymag,zmag) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gsm_to_sm(xgsm,ygsm,zgsm,xsma,ysma,zsma) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsm_to_sma: GSM -> SM   system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! *   Output : xsma,ysma,zsma   cartesian sma coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xsma= cmu*xgsm - smu*zgsm 
      ysma= ygsm 
      zsma= smu*xgsm + cmu*zgsm 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gsm_to_tpn(xgsm,ygsm,zgsm,xo,yo,zo,xs,xtpn,ytpn,ztpn) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsm_to_tpn: GSM -> TPN  system                
! *   Author : P. Robert, LPP , 2016                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xgsm,ygsm,zgsm cartesian gsm vector                      
!              xo,yo,zo position of the S/C in gsm                      
!              xs subsolar point, submit of the paraboloid              
!              from Earth to Sun                                        
! *   Output : xtpn,ytpn,ztpn cartesian tpn coordinates                 
! ----------------------------------------------------------------------
!                                                                       
      real Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz 
!                                                                       
!     computation of the TPN axis in gsm system                         
!                                                                       
      call cp_tpn_param(xo,yo,zo,xs, Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz) 
!                                                                       
      xtpn= Tx*xgsm + Ty*ygsm +Tz*zgsm 
      ytpn= Px*xgsm + Py*ygsm +Pz*zgsm 
      ztpn= Nx*xgsm + Ny*ygsm +Nz*zgsm 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gseq_to_gei(xgsq,ygsq,zgsq,xgei,ygei,zgei) 
!                                                                       
! ----------------------------------------------------------------------
!                                                                       
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsq_to_gei: GSEQ-> GEI  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgsq,ygsq,zgsq   cartesian gsq coordinates               
! *   Output : xgei,ygei,zgei   cartesian gei coordinates               
!                                                                       
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3 
!                                                                       
      xgei= gs1*xgsq + yeigq1*ygsq + zeigq1*zgsq 
      ygei= gs2*xgsq + yeigq2*ygsq + zeigq2*zgsq 
      zgei= gs3*xgsq + yeigq3*ygsq + zeigq3*zgsq 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gseq_to_geo(xgsq,ygsq,zgsq,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsq_to_geo: GSEQ-> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgsq,ygsq,zgsq   cartesian gsq coordinates               
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3 
!                                                                       
      xgeo= ps1*xgsq + yeogq1*ygsq + zeogq1*zgsq 
      ygeo= ps2*xgsq + yeogq2*ygsq + zeogq2*zgsq 
      zgeo= ps3*xgsq + yeogq3*ygsq + zeogq3*zgsq 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gseq_to_gse(xgsq,ygsq,zgsq,xgse,ygse,zgse) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsq_to_gse: GSEQ-> GSE  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xgsq,ygsq,zgsq   cartesian gsq coordinates               
! *   Output : xgse,ygse,zgse   cartesian gse coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xgse= xgsq 
      ygse= ctetq*ygsq + stetq*zgsq 
      zgse=-stetq*ygsq + ctetq*zgsq 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_gseq_to_gsm(xgsq,ygsq,zgsq,xgsm,ygsm,zgsm) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_gsq_to_gsm: GSQ -> GSM  system                
! *   Author : P. Robert, CRPE, 2002                                    
!                                                                       
! *   Input  : xgsq,ygsq,zgsq cartesian gsq coordinates                 
! *   Output : xgsm,ygsm,zgsm cartesian gsm coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      call t_gseq_to_gse(xgsq,ygsq,zgsq,xx,yy,zz) 
      call t_gse_to_gsm(xx,yy,zz,xgsm,ygsm,zgsm) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_mag_to_gei(xmag,ymag,zmag,xgei,ygei,zgei) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_mag_to_gei: MAG -> GEI  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xmag,ymag,zmag cartesian mag coordinates                 
! *   Output : xgei,ygei,zgei cartesian gei coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3 
!                                                                       
      xgei= xeima1*xmag + yeima1*ymag + gm1*zmag 
      ygei= xeima2*xmag + yeima2*ymag + gm2*zmag 
      zgei= xeima3*xmag + yeima3*ymag + gm3*zmag 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_mag_to_geo(xmag,ymag,zmag,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_mag_to_geo: MAG -> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xmag,ymag,zmag   cartesian mag coordinates               
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3 
!                                                                       
      xgeo= xeoma1*xmag + yeoma1*ymag + gd1*zmag 
      ygeo= xeoma2*xmag + yeoma2*ymag + gd2*zmag 
      zgeo= xeoma3*xmag + yeoma3*ymag + gd3*zmag 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_mag_to_gsm(xmag,ymag,zmag,xgsm,ygsm,zgsm) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_mag_to_gsm: MAG -> GSM  system                
! *   Author : P. Robert, CRPE, 2002                                    
!                                                                       
! *   Input  : xmag,ymag,zmag cartesian mag coordinates                 
! *   Output : xgsm,ygsm,zgsm cartesian gsm coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      call t_mag_to_sm(xmag,ymag,zmag,xx,yy,zz) 
      call t_sm_to_gsm(xx,yy,zz,xgsm,ygsm,zgsm) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_mag_to_sm(xmag,ymag,zmag,xsma,ysma,zsma) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_mag_to_sma: MAG -> SM   system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xmag,ymag,zmag   cartesian mag coordinates               
! *   Output : xsma,ysma,zsma   cartesian sma coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xsma= cphi*xmag - sphi*ymag 
      ysma= sphi*xmag + cphi*ymag 
      zsma= zmag 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sm_to_gei(xsma,ysma,zsma,xgei,ygei,zgei) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sma_to_gei: SM  -> GEI  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xsma,ysma,zsma   cartesian sma coordinates               
! *   Output : xgei,ygei,zgei   cartesian gei coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3 
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3 
!                                                                       
      xgei=xeism1*xsma + yeism1*ysma + gm1*zsma 
      ygei=xeism2*xsma + yeism2*ysma + gm2*zsma 
      zgei=xeism3*xsma + yeism3*ysma + gm3*zsma 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sm_to_geo(xsma,ysma,zsma,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sma_to_geo: SM  -> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xsma,ysma,zsma   cartesian sma coordinates               
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3 
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3 
!                                                                       
      xgeo= xeosm1*xsma + yeosm1*ysma + gd1*zsma 
      ygeo= xeosm2*xsma + yeosm2*ysma + gd2*zsma 
      zgeo= xeosm3*xsma + yeosm3*ysma + gd3*zsma 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sm_to_gsm(xsma,ysma,zsma,xgsm,ygsm,zgsm) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sma_to_gsm: SM  -> GSM  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xsma,ysma,zsma   cartesian sma coordinates               
! *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xgsm=  cmu*xsma + smu*zsma 
      ygsm=  ysma 
      zgsm= -smu*xsma + cmu*zsma 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sm_to_mag(xsma,ysma,zsma,xmag,ymag,zmag) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sma_to_mag: SM  -> MAG  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: terms of transformation matrix are given in common       
!                                                                       
! *   Input  : xsma,ysma,zsma   cartesian sma coordinates               
! *   Output : xmag,ymag,zmag   cartesian mag coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi 
!                                                                       
      xmag=  cphi*xsma + sphi*ysma 
      ymag= -sphi*xsma + cphi*ysma 
      zmag=  zsma 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sph_to_car(r,teta,phi,x,y,z) 
!                                                                       
! ----------------------------------------------------------------------
!                                                                       
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sph_to_car: SPH -> CAR  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: none                                                     
!                                                                       
! *   Input  :   r,teta,phi  spherical coordinates (angles in radians)  
! *   Output :   x,y,z       cartesian coordinates                      
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      sq=r*sin(teta) 
      x=sq*cos(phi) 
      y=sq*sin(phi) 
      z=r*cos(teta) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sr2_to_gse(xsr2,ysr2,zsr2,rotx,roty,rotz,            &
     &                        xgse,ygse,zgse)                           
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sr2_to_gse: SR2 -> GSE  system                
! *   Author : P. Robert, CETP, 2001                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xsr2,ysr2,sr2z cartesian sr2 coordinates                 
!              rotx,roty,rotz cartesian gse coordinates of rotation axis
! *   Output : xgse,ygse,zgse cartesian gse coordinates                 
! ----------------------------------------------------------------------
!                                                                       
! *** set transform matrix with Spin axis terms                         
!                                                                       
      rmod= sqrt(rotx**2 + roty**2 + rotz**2) 
!                                                                       
      rx=rotx/rmod 
      ry=roty/rmod 
      rz=rotz/rmod 
!                                                                       
      a= 1./sqrt(ry*ry + rz*rz) 
!                                                                       
      x1= (ry*ry + rz*rz)*a 
      x2= -rx*ry*a 
      x3= -rx*rz*a 
!                                                                       
      y1=  0. 
      y2=  rz*a 
      y3= -ry*a 
!                                                                       
      z1= rx 
      z2= ry 
      z3= rz 
!                                                                       
! *** tranform the input vector from sr2 to gse                         
!                                                                       
      xgse= x1*xsr2 + y1*ysr2 + z1*zsr2 
      ygse= x2*xsr2 + y2*ysr2 + z2*zsr2 
      zgse= x3*xsr2 + y3*ysr2 + z3*zsr2 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sr2_to_mfa(xsr2,ysr2,zsr2,bx,by,bz,rox,roy,roz,      &
     &                        xm,ym,zm)                                 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sr2_to_mfa: SR2 -> MFA  system                
! *   Author : P. Robert, CETP, 2001                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xsr2,ysr2,sr2z cartesian sr2 coordinates                 
!              bx,  by,  bz   cartesian sr2 coordinates of DC mag field 
!              rox, roy, roz  cartesian gse coordinates of rotation axis
! *   Output : xm,  ym,  zm   cartesian mfa coordinates                 
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
! *** normalise R supposed undependant of time                          
!                                                                       
      rmod= sqrt(rox**2 + roy**2 + roz**2) 
!                                                                       
      rx=rox/rmod 
      ry=roy/rmod 
      rz=roz/rmod 
!                                                                       
! *** tranform the vector from sr2 to mfa                               
!                                                                       
      bperp= sqrt(bx*bx + by*by) 
      b0= sqrt(bx*bx + by*by + bz*bz) 
!                                                                       
! **  first rotation                                                    
!                                                                       
      sinphi=by/bperp 
      cosphi=bx/bperp 
!                                                                       
      xp=  cosphi*xsr2 + sinphi*ysr2 
      yp= -sinphi*xsr2 + cosphi*ysr2 
      zp=  zsr2 
!                                                                       
! **  second rotation                                                   
!                                                                       
      sintet= bperp/b0 
      costet= bz/b0 
!                                                                       
      xmag= costet*xp - sintet*zp 
      ymag= yp 
      zmag= sintet*xp + costet*zp 
!                                                                       
! **  third rotation                                                    
!                                                                       
      b= sqrt(ry*ry + rz*rz) 
!                                                                       
      sxm=  b*costet*cosphi - rx*sintet 
      sym= -b*sinphi 
!                                                                       
      smperp= sqrt(sxm*sxm + sym*sym) 
!                                                                       
      sinpsi= sym/smperp 
      cospsi= sxm/smperp 
!                                                                       
      xm=  cospsi*xmag + sinpsi*ymag 
      ym= -sinpsi*xmag + cospsi*ymag 
      zm=  zmag 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sr2_to_sr(xsr2,ysr2,spifre,spipha,deltaT,xsre,ysre) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sr2_to_sre: SR2 -> SRef system                
! *   Author : P. Robert, CRPE, 2001                                    
! *   Comment: Z component is unchanged (spin axis)                     
!                                                                       
! *   Input  : xsr2, ysr2 cartesian sr2 coordinates                     
!              spifre     spin frequency in Hz                          
!              spipha     spin phase in radians, growing with time      
!                         spipha= positive angle between the xsr axis   
!                         and the component of the direction of the Sun 
!                         in the xsr-ysr plane.                         
!              deltaT     (T -To) time (sec.), between the current time 
!                         and the time where is measured the spin phase.
! *   Output : xsre,ysre  cartesian sr coordinates                      
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      pi2= 2.*3.1415927 
!                                                                       
      depift= mod(pi2*spifre*deltaT,pi2) 
      phicr= spipha -depift 
      sinphi=sin(phicr) 
      cosphi=cos(phicr) 
!                                                                       
      xsre=  cosphi*xsr2 -sinphi*ysr2 
      ysre=  sinphi*xsr2 +cosphi*ysr2 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_sr_to_sr2(xsre,ysre,spifre,spipha,deltaT,xsr2,ysr2) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_sre_to_sr2: SRef-> SR2 system                 
! *   Author : P. Robert, CRPE, 2001                                    
! *   Comment: Z component is unchanged (spin axis)                     
!                                                                       
! *   Input  : xsre, ysre cartesian sr coordinates                      
!              spifre     spin frequency in Hz                          
!              spipha     spin phase in radians, growing with time      
!                         spipha= positive angle between the xsr axis   
!                         and the component of the direction of the Sun 
!                         in the xsr-ysr plane.                         
!              deltaT     (T -To) time (sec.), between the current time 
!                         and the time where is measured the spin phase.
! *   Output : xsr2,ysr2  cartesian sr coordinates                      
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      pi2= 2.*3.1415927 
!                                                                       
      depift= mod(pi2*spifre*deltaT,pi2) 
      phicr= spipha -depift 
      sinphi=sin(phicr) 
      cosphi=cos(phicr) 
!                                                                       
      xsr2=  cosphi*xsre +sinphi*ysre 
      ysr2= -sinphi*xsre +cosphi*ysre 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_vdh_to_geo(xvdh,yvdh,zvdh,rlat,rlong,xgeo,ygeo,zgeo) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_vdh_to_geo: VDH -> GEO  system                
! *   Author : P. Robert, CRPE, 1992                                    
! *   Comment: local system, non time dependent                         
!                                                                       
! *   Input  : xvdh,yvdh,zvdh   cartesian vdh coordinates               
!              rlat,rlong       latitude and longitude of the point     
!                               of observation (radians)                
! *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates               
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      q=cos(rlat) 
      r1=q*cos(rlong) 
      r2=q*sin(rlong) 
      r3=sin(rlat) 
!                                                                       
      v1=r1 
      v2=r2 
      v3=r3 
!                                                                       
      q12=sqrt(r1**2 + r2**2) 
!                                                                       
      d1= -r2/q12 
      d2=  r1/q12 
      d3=  0. 
!                                                                       
      h1= -r1*r3/q12 
      h2= -r2*r3/q12 
      h3=  q12 
!                                                                       
      xgeo= v1*xvdh + d1*yvdh + h1*zvdh 
      ygeo= v2*xvdh + d2*yvdh + h2*zvdh 
      zgeo= v3*xvdh + d3*yvdh + h3*zvdh 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_xyz_to_vdh(x,y,z,a1,a2,a3,v,d,h) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_xyz_to_vdh: xyz spinning -> VDH               
! *   Author : P. Robert, SDev, 2020                                    
! *   Comment: use Euler angles in degrees; VDH is the fixed system.    
!                                                                       
! *   Input : x,y,z cartesian xyz coordinates                           
! *   Output: v,d,h cartesian VDH coordinates                           
! ----------------------------------------------------------------------
                                                                        
      pisd=3.1415926/180. 
!                                                                       
      a1r=a1*pisd 
      a2r=a2*pisd 
      a3r=a3*pisd 
!                                                                       
      c1=cos(a1r) 
      c2=cos(a2r) 
      c3=cos(a3r) 
                                                                        
      s1=sin(a1r) 
      s2=sin(a2r) 
      s3=sin(a3r) 
!                                                                       
      v=( c1*c3 -s1*c2*s3)*x  +(-s1*c3 -c1*c2*s3)*y  +( s2*s3)*z 
      d=( c1*s3 +s1*c2*c3)*x  +(-s1*s3 +c1*c2*c3)*y  +(-s2*c3)*z 
      h=           (s1*s2)*x  +           (c1*s2)*y  +(    c2)*z 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine t_vdh_to_xyz(v,d,h,a1,a2,a3,x,y,z) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : transform modules of Rocotlib Software                   
! *   Object : transforms_vdh_to_xyz: VDH -> xyz spinning               
! *   Author : P. Robert, SDev, 2020                                    
! *   Comment: use Euler angles in degrees; VDH is the fixed system.    
!                                                                       
! *   Input  : v,d,h cartesian VDH coordinates                          
! *   Output : x,y,z cartesian xyz coordinates                          
! ----------------------------------------------------------------------
                                                                        
      pisd=3.1415926/180. 
!                                                                       
      a1r=a1*pisd 
      a2r=a2*pisd 
      a3r=a3*pisd 
!                                                                       
      c1=cos(a1r) 
      c2=cos(a2r) 
      c3=cos(a3r) 
                                                                        
      s1=sin(a1r) 
      s2=sin(a2r) 
      s3=sin(a3r) 
!                                                                       
      x=( c1*c3 -s1*c2*s3)*v  +( c1*s3 +s1*c2*c3)*d +( s1*s2)*h 
      y=(-s1*c3 -c1*c2*s3)*v  +(-s1*s3 +c1*c2*c3)*d +( c1*s2)*h 
      z=          ( s2*s3)*v  +          (-s2*c3)*d +(    c2)*h 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
!     Supplements for matrix operations.                                
!     Allows the transformation of an input 3D signal into the          
!     Minimum variance coordinate (MVA Analysis).                       
!                                                                       
!     All codes in this part are extracted of Roproc Software, V 4.5    
!     and rewrited properly in f77, with Rocotlib V2.2 conventions.     
!     P. Robert, LPP , April 2016.                                      
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_cp_varmin(ifc,Vx,Vy,Vz,n,irep,covar,lambda,eigvec) 
!                                                                       
      real Vx(n),Vy(n),Vz(n) 
      real covar(3,3),lambda(3),eigvec(3,3) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : compute variance minimum coordinate of a signal Vx,Vy,Vz 
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!     Comment: check orthogonality of the MVA matrix                    
!                                                                       
! *   Input  : ifc (unit for writing results)                           
!              Vx(n),Vy(n),Vz(n) input signal                           
!              n number of point of the signal                          
!              irep=1 : passe input vector in MVA coordinates           
!                                                                       
! *   Output : covar covariance matrix                                  
!              lamda(3), eigvev(3,3) eigen values and eigen vectors     
!              results are also writted on unit ifc                     
!     ---------------------------------------------------------------+--
!                                                                       
      write(ifc,100) 
      write(ifc,100) '------------------------------------------------' 
      write(ifc,100) 'Compute Minimum Variance Analaysis coord. system' 
      write(ifc,100) '------------------------------------------------' 
!                                                                       
! *** compute covariance matrix                                         
!                                                                       
      call mat_cp_covariance(Vx,Vy,Vz,n,covar) 
      call mat_write(ifc,'covariance matrix of input signal:',covar) 
!                                                                       
! *** diagonalisation of covariance matrix                              
!                                                                       
      call mat_diagonalise(covar,lambda,eigvec) 
      call mat_write(ifc,'covariance matrix diagonalized:',covar) 
      call mat_write(ifc,'eigen vectors matrix:',eigvec) 
!                                                                       
! *** check orthogonality and direct sens of the 3 vectors              
!                                                                       
      call mat_check_ortho(ifc,eigvec) 
!                                                                       
! *** print eigen vectors and eigen values                              
!                                                                       
      call mat_write_eigen_vec(ifc,lambda,eigvec) 
!                                                                       
! *** transform input signal in MVA coordinates, ie. eigen vectors      
!                                                                       
      if(irep.ne.0) then 
        call mat_change_coord(eigvec,Vx,Vy,Vz,n) 
        write(ifc,100) 
        write(ifc,100) 'input signal is passed into MVA coordinates' 
      endif 
!                                                                       
  100 format(1x,a) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_cp_covariance(Vx,Vy,Vz,n,covar) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : compute covariance matrix for a vector series V(n)       
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : Vx(n),Vy(n),Vz(n) vector series                          
! *   Output : covar(3,3)                                               
!     ---------------------------------------------------------------+--
!                                                                       
      dimension Vx(n),Vy(n),Vz(n) 
      dimension covar(3,3) 
      double precision sx, sy, sz, sx2, sy2, sz2, sxy, syz, szx 
!                                                                       
! *   set summation to zero                                             
!                                                                       
      sx=0.d0 
      sy=0.d0 
      sz=0.d0 
                                                                        
      sx2=0.d0 
      sy2=0.d0 
      sz2=0.d0 
                                                                        
      sxy=0.d0 
      syz=0.d0 
      szx=0.d0 
                                                                        
!                                                                       
! *   components summation and square components summation              
!                                                                       
      do k=1,n 
!                                                                       
        sx= sx +dble(Vx(k)) 
        sy= sy +dble(Vy(k)) 
        sz= sz +dble(Vz(k)) 
!                                                                       
        sx2= sx2 +dble(Vx(k)**2) 
        sy2= sy2 +dble(Vy(k)**2) 
        sz2= sz2 +dble(Vz(k)**2) 
!                                                                       
        sxy= sxy +dble(Vx(k)*Vy(k)) 
        syz= syz +dble(Vy(k)*Vz(k)) 
        szx= szx +dble(Vz(k)*Vx(k)) 
!                                                                       
      enddo 
!                                                                       
! *   average values                                                    
!                                                                       
      sx = sx/dble(n) 
      sy = sy/dble(n) 
      sz = sz/dble(n) 
!                                                                       
      sx2= sx2/dble(n) 
      sy2= sy2/dble(n) 
      sz2= sz2/dble(n) 
!                                                                       
      sxy= sxy/dble(n) 
      syz= syz/dble(n) 
      szx= szx/dble(n) 
!                                                                       
! *   variance set to diagonals terms                                   
!                                                                       
      covar(1,1)= sngl(sx2-sx**2) 
      covar(2,2)= sngl(sy2-sy**2) 
      covar(3,3)= sngl(sz2-sz**2) 
!                                                                       
! *   covariances set to semi \diagonal terms                           
!                                                                       
      covar(2,1)= sngl(sxy-sx*sy) 
      covar(3,2)= sngl(syz-sy*sz) 
      covar(3,1)= sngl(szx-sz*sx) 
!                                                                       
! *   second half triangle                                              
!                                                                       
      covar(1,2)= covar(2,1) 
      covar(1,3)= covar(3,1) 
      covar(2,3)= covar(3,2) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_diagonalise(mat,lambda,eigvec) 
!                                                                       
      real mat(3,3),lambda(3),eigvec(3,3) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : diagonalise the given matrix mat(3,3)                    
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : real mat(3,3)                                            
! *   Output : lambda(3)   eigen values                                 
!              eigvec(3,3) eigen vectors normalised to 1.               
!     ---------------------------------------------------------------+--
!                                                                       
!     compute eigen vectors and eigen values of real mat(3,3)           
!                                                                       
      call mat_cp_eigen_vec(mat,lambda,eigvec) 
!                                                                       
!     normalise to 1 eigen vectors                                      
                                                                        
      call mat_normalize_vec(eigvec) 
                                                                        
! *   set positive direction for \diag                                  
!                                                                       
      if (eigvec(3,1).lt.0.) then 
          eigvec(1,1)= -eigvec(1,1) 
          eigvec(2,1)= -eigvec(2,1) 
          eigvec(3,1)= -eigvec(3,1) 
      endif 
                                                                        
      if (eigvec(1,3).lt.0.) then 
          eigvec(1,3)= -eigvec(1,3) 
          eigvec(2,3)= -eigvec(2,3) 
          eigvec(3,3)= -eigvec(3,3) 
      endif 
!                                                                       
! *   sort eigen vector with eigen value (z=min)                        
!                                                                       
      do k=1,4 
         i=k 
         if (k.gt.2) i=k-2 
                                                                        
         if (lambda(i).lt.lambda(i+1)) then 
            dint=lambda(i) 
            lambda(i)=lambda(i+1) 
            lambda(i+1)=dint 
            do j=1,3 
               dint=eigvec(j,i) 
               eigvec(j,i)=eigvec(j,i+1) 
               eigvec(j,i+1)=dint 
            enddo 
         endif 
      enddo 
!                                                                       
! *   determinant computation:                                          
!                                                                       
      call mat_cp_determin(eigvec,det) 
!                                                                       
      if (det.lt.0.) then 
         eigvec(1,2)= -eigvec(1,2) 
         eigvec(2,2)= -eigvec(2,2) 
         eigvec(3,2)= -eigvec(3,2) 
      endif 
!                                                                       
!                                                                       
! *   modif P.R. 2004: the vector corresponding to lambda min has z >0  
!                                                                       
      if(eigvec(3,3).lt.0.) then 
         eigvec(1,2)= -eigvec(1,2) 
         eigvec(2,2)= -eigvec(2,2) 
         eigvec(3,2)= -eigvec(3,2) 
                                                                        
         eigvec(1,3)= -eigvec(1,3) 
         eigvec(2,3)= -eigvec(2,3) 
         eigvec(3,3)= -eigvec(3,3) 
      endif 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_check_ortho(ifc,mat) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : check orthogonality of matrix components                 
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : ifc (unit for writing results), real mat(3,3)            
! *   Output : writing result on unit ifc                               
!     ---------------------------------------------------------------+--
!                                                                       
      real mat(3,3) 
!                                                                       
      v1x=mat(1,1) 
      v1y=mat(2,1) 
      v1z=mat(3,1) 
!                                                                       
      v2x=mat(1,2) 
      v2y=mat(2,2) 
      v2z=mat(3,2) 
!                                                                       
      v3x=mat(1,3) 
      v3y=mat(2,3) 
      v3z=mat(3,3) 
!                                                                       
      write(ifc,100) 
      write(ifc,100) 'Check orthogonality of given matrix:' 
      write(ifc,100) '----------------------------------- ' 
      write(ifc,100) 
!                                                                       
      write(ifc,100) '1) dot product of any two rows or any two colums '&
     &               //'must be equal to zero:'                         
                                                                        
      v1v2= v1x*v2x + v1y*v2y + v1z*v2z 
      v2v3= v2x*v3x + v2y*v3y + v2z*v3z 
      v3v1= v3x*v1x + v3y*v1y + v3z*v1z 
!                                                                       
      write(ifc,100) 
      write(ifc,100) 'V1.V2 =',v1v2 
      write(ifc,100) 'V2.V3 =',v2v3 
      write(ifc,100) 'V3.V1 =',v3v1 
!                                                                       
      write(ifc,100) 
      write(ifc,100) '2) cross product of any two rows or colums '//    &
     &                  'must be equal to the third '                   
      write(ifc,100) '   row or column or its negative):' 
!                                                                       
! *   v1 X v2                                                           
!                                                                       
      v1cv2x= v1y*v2z - v1z*v2y 
      v1cv2y= v1z*v2x - v1x*v2z 
      v1cv2z= v1x*v2y - v1y*v2x 
!                                                                       
! *   v2 X v3                                                           
!                                                                       
      v2cv3x= v2y*v3z - v2z*v3y 
      v2cv3y= v2z*v3x - v2x*v3z 
      v2cv3z= v2x*v3y - v2y*v3x 
!                                                                       
! *   v3 X v1                                                           
!                                                                       
      v3cv1x= v3y*v1z - v3z*v1y 
      v3cv1y= v3z*v1x - v3x*v1z 
      v3cv1z= v3x*v1y - v3y*v1x 
!                                                                       
      write(ifc,100) 
      write(ifc,200) 'V1XV2 =',v1cv2x,v1cv2y,v1cv2z 
      write(ifc,200) '   V3 =',v3x,v3y,v3z 
!                                                                       
      write(ifc,100) 
      write(ifc,200) 'V2XV3 =',v2cv3x,v2cv3y,v2cv3z 
      write(ifc,200) '   V1 =',v1x,v1y,v1z 
!                                                                       
      write(ifc,100) 
      write(ifc,200) 'V3XV1 =',v3cv1x,v3cv1y,v3cv1z 
      write(ifc,200) '   V2 =',v2x,v2y,v2z 
!                                                                       
  100 format(3(1x,a,1Pe14.5)) 
  200 format(1x,a,3(1Pe16.5)) 
                                                                        
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_cp_determin(mat,det) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : compute determinant of the given matrix                  
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : real mat(3,3)                                            
! *   Output : det                                                      
!     ---------------------------------------------------------------+--
!                                                                       
      real mat(3,3) 
!                                                                       
      det=    mat(1,1)*mat(2,2)*mat(3,3)+mat(1,2)*mat(2,3)*mat(3,1) 
      det=det+mat(2,1)*mat(3,2)*mat(1,3)-mat(3,1)*mat(2,2)*mat(1,3) 
      det=det-mat(1,1)*mat(3,2)*mat(2,3)-mat(2,1)*mat(1,2)*mat(3,3) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_cp_eigen_vec(mat,lambda,eigvec) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : compute eigen vectors and eigen values of real mat(3,3)  
! *   Author : unknown, CETP, 2001, rev. PR 2016                        
! *   Comment: mat(3,3) must be real and symmetrical;                   
!              Method used is Householder.                              
!                                                                       
! *   Input  : real mat(3,3)                                            
! *   Output : lambda(3)   eigen values                                 
!              eigvec(3,3) eigen vectors                                
!     ---------------------------------------------------------------+--
!                                                                       
      real mat(3,3),lambda(3),eigvec(3,3),rr(3,3),e(3) 
!                                                                       
! *   set eigen vectors to unity matrix                                 
!                                                                       
      do i=1,3 
         do j=1,3 
            eigvec(i,j)=0.0 
         enddo 
         eigvec(i,i)=1.0 
      enddo 
!                                                                       
!     set eigen values                                                  
!                                                                       
      lambda(1)=mat(1,3) 
      lambda(2)=mat(2,3) 
!                                                                       
      signe=1.0 
      if (lambda(2).lt.0.0) signe=-1.0 
      lambda(2)=lambda(2)+signe*sqrt(lambda(1)**2+lambda(2)**2) 
      lambda(3)=0.0 
!                                                                       
      h=(lambda(1)**2+lambda(2)**2)/2.0 
!                                                                       
      epsi=1.e-30 
                                                                        
      if (abs(h).gt.epsi) then 
!                                                                       
!        determination of the Householder matrix (results in  eigvec).  
!                                                                       
         do i=1,3 
            do j=1,3 
               eigvec(i,j)=eigvec(i,j)-(lambda(i)*lambda(j))/h 
            enddo 
         enddo 
!                                                                       
!        tridiagonalisation of mat matrix                               
!                                                                       
         call mat_product(eigvec,mat,rr) 
         call mat_product(rr,eigvec,mat) 
!                                                                       
      endif 
                                                                        
!     set eigen values to diagonal terms                                
!                                                                       
      do i=1,3 
         lambda(i)=mat(i,i) 
      enddo 
!                                                                       
      do i=1,2 
         e(i)=mat(i+1,i) 
      enddo 
      e(3)=0.0 
!                                                                       
! *   search for eigen values (max iteration = 30)                      
!                                                                       
      do l=1,3 
         iter=0 
    1    do m=l,2 
            dd=abs(lambda(m))+abs(lambda(m+1)) 
! PR        if (abs(e(m))+dd.eq.dd) goto 2                              
            if (abs(e(m)).le.abs(dd)*1E-6) goto 2 
         enddo 
         m=3 
    2    if (m.ne.l) then 
!                                                                       
            if (iter.eq.30) then 
               stop '*** Rocotlib/mat_cp_eigen_vec: ABORTED ! '//       &
     &                   'TOO MANY ITERATIONS'                          
            endif 
                                                                        
            iter=iter+1 
            g=(lambda(l+1)-lambda(l))/(2.0*e(l)) 
            call mat_cp_pythag_func(g,1.0,r) 
            g=lambda(m)-lambda(l)+e(l)/(g+sign(r,g)) 
            s=1.0 
            c=1.0 
            p=0.0 
                                                                        
            do i=m-1,l,-1 
               f=s*e(i) 
               b=c*e(i) 
               call mat_cp_pythag_func(f,g,r) 
               e(i+1)=r 
                                                                        
               if (abs(r).lt.epsi) then 
                  lambda(i+1)=lambda(i+1)-p 
                  e(m)=0.0 
                  goto 1 
               endif 
                                                                        
               s=f/r 
               c=g/r 
               g=lambda(i+1)-p 
               r=(lambda(i)-g)*s+2.0*c*b 
               p=s*r 
               lambda(i+1)=g+p 
               g=c*r-b 
!                                                                       
               do k=1,3 
                  f=eigvec(k,i+1) 
                  eigvec(k,i+1)=s*eigvec(k,i)+c*f 
                  eigvec(k,i)=c*eigvec(k,i)-s*f 
               enddo 
!                                                                       
            enddo 
                                                                        
            lambda(l)=lambda(l)-p 
            e(l)=g 
            e(m)=0.0 
            goto 1 
!                                                                       
         endif 
      enddo 
!                                                                       
! *   modif PR  sep 11 2001 to load  mat with the result of the diag.   
!                                                                       
      do i=1,3 
         mat(i,i)=lambda(i) 
      enddo 
!                                                                       
      do i=1,2 
         mat(i+1,i)=e(i) 
         mat(i,i+1)=e(i) 
      enddo 
!                                                                       
      mat(1,3)=e(3) 
      mat(3,1)=e(3) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_cp_pythag_func(a,b,fpyth) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : Pythagore function of two real (used by mat_cp_eigen_vec)
! *   Author : unknown, CETP, 2001, rev. PR 2016                        
!                                                                       
! *   Input  : a,b                                                      
! *   Output : fpyth                                                    
!     ---------------------------------------------------------------+--
!                                                                       
      absa=abs(a) 
      absb=abs(b) 

      epsi=1.e-37
                                                                        
      if (absa.gt.absb)  then 
             fpyth=absa*sqrt(1.0+(absb/absa)**2) 
                         else 
         if (absb.lt.epsi) then 
             fpyth=0.0 
                          else 
             fpyth=absb*sqrt(1.0+(absa/absb)**2) 
         endif 
      endif 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_normalize_vec(mat) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : normalize to 1. the vectors of the input matrix          
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : real mat(3,3)                                            
! *   Output : real mat(3,3) with vectors normalized to 1.              
!     ---------------------------------------------------------------+--
!                                                                       
      real mat(3,3),matmod 

      epsi=1.e-37
!                                                                       
      do j=1,3 
         matmod=0. 
         do i=1,3 
           matmod=matmod+mat(i,j)*mat(i,j) 
         enddo 
         do i=1,3 
            if (abs(matmod).gt.epsi) then 
               mat(i,j)=mat(i,j)/sqrt(matmod) 
            else 
                print*, '*** Rocotlib/mat_normalize_vec: module is = 0' 
            endif 
         enddo 
      enddo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_product(mat1,mat2,mat3) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : matrix product of two given matrix of dim. 3             
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : n, mat1(3,3), mat2(3,3)                                  
! *   Output : mat3(3,3)                                                
!     ---------------------------------------------------------------+--
!                                                                       
      real mat1(3,3),mat2(3,3),mat3(3,3) 
!                                                                       
      do i=1,3 
         do j=1,3 
            mat3(i,j)=0.0 
            do k=1,3 
               mat3(i,j)=mat3(i,j)+mat1(i,k)*mat2(k,j) 
            enddo 
         enddo 
      enddo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_somme(mat1,mat2,mat3) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : matrix somme of two given matrix of dim. 3               
! *   Author : P. Robert, LPP , 2016                                    
!                                                                       
! *   Input  : n, mat1(3,3), mat2(3,3)                                  
! *   Output : mat3(3,3)                                                
!     ---------------------------------------------------------------+--
!                                                                       
      real mat1(3,3),mat2(3,3),mat3(3,3) 
!                                                                       
      do i=1,3 
         do j=1,3 
            mat3(i,j)=mat1(i,j)+mat2(i,j) 
         enddo 
      enddo 
                                                                        
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_diff(mat1,mat2,mat3) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : matrix difference of two given matrix                    
! *   Author : P. Robert, LPP , 2016                                    
!                                                                       
! *   Input  : n, mat1(3,3), mat2(3,3)                                  
! *   Output : mat3(3,3)                                                
!     ---------------------------------------------------------------+--
!                                                                       
      real mat1(3,3),mat2(3,3),mat3(3,3) 
!                                                                       
      do i=1,3 
         do j=1,3 
            mat3(i,j)=mat1(i,j)-mat2(i,j) 
         enddo 
      enddo 
                                                                        
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_transpose(mat) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : transpose input matrix                                   
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : real mat(3,3)                                            
! *   Output : real mat(3,3) transposed                                 
!     ---------------------------------------------------------------+--
!                                                                       
      real mat(3,3) 
!                                                                       
      r12= mat(1,2) 
      r13= mat(1,3) 
      r23= mat(2,3) 
!                                                                       
      mat(1,2)=mat(2,1) 
      mat(1,3)=mat(3,1) 
      mat(2,3)=mat(3,2) 
!                                                                       
      mat(2,1)=r12 
      mat(3,1)=r13 
      mat(3,2)=r23 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_change_coord(mat,Vx,Vy,Vz,n) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : change coordinate of a vector serie with a given matrix  
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
! *   Comment: compute V(n)= mat*V(n)                                   
!                                                                       
! *   Input  : ifc (unit for writing results), real mat(3,3)            
! *   Output : writing result on unit ifc                               
!     ---------------------------------------------------------------+--
!                                                                       
      real Vx(n),Vy(n),Vz(n) 
      real mat(3,3) 
!                                                                       
      do i=1,n 
!                                                                       
         Vxr= mat(1,1)*Vx(i) +mat(2,1)*Vy(i) +mat(3,1)*Vz(i) 
         Vyr= mat(1,2)*Vx(i) +mat(2,2)*Vy(i) +mat(3,2)*Vz(i) 
         Vzr= mat(1,3)*Vx(i) +mat(2,3)*Vy(i) +mat(3,3)*Vz(i) 
!                                                                       
         Vx(i)=Vxr 
         Vy(i)=Vyr 
         Vz(i)=Vzr 
!                                                                       
      enddo 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_write(ifc,com,mat) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : print on ifc unit mat(3,3) with a comment                
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : ifc (unit for writing results), com, real mat(3,3)       
! *   Output : writing result on unit ifc                               
!     ---------------------------------------------------------------+--
!                                                                       
      real mat(3,3) 
      character*(*) com 
!                                                                       
      write(ifc,100) 
      write(ifc,100) com 
      write(ifc,100) 
      write(ifc,200) (mat(1,j),j=1,3) 
      write(ifc,200) (mat(2,j),j=1,3) 
      write(ifc,200) (mat(3,j),j=1,3) 
!                                                                       
  100 format(1x,a) 
  200 format(11x,3(1Pe15.5)) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
      subroutine mat_write_eigen_vec(ifc,lambda,mat) 
!                                                                       
!     ---------------------------------------------------------------+--
! *   Class  : matrix operation of Rocotlib Software                    
! *   Object : print on ifc unit eigen values & vectors of mat(3,3)     
! *   Author : P. Robert, CETP, 2001, rev. PR 2016                      
!                                                                       
! *   Input  : ifc (unit for writing results), lambda(3), real mat(3,3) 
! *   Output : writing result on unit ifc                               
!     ---------------------------------------------------------------+--
!                                                                       
      real lambda(3),mat(3,3),r(3),theta(3),phi(3) 
!                                                                       
      dpi=180./3.14159 
!                                                                       
      do i=1,3 
         call t_car_to_sph(mat(1,i),mat(2,i),mat(3,i),r(i),teta,phir) 
         theta(i)=teta*dpi 
         phi(i)  =phir*dpi 
      enddo 
                                                                        
      write(ifc,100) 
      write(ifc,100) 'Eigen vectors and eigen values' 
      write(ifc,100) '------------------------------' 
      write(ifc,100) 
      write(ifc,110) '   ', 'V1','V2','V3' 
      write(ifc,100) 'x     ', mat(1,1),mat(1,2),mat(1,3) 
      write(ifc,100) 'y     ', mat(2,1),mat(2,2),mat(2,3) 
      write(ifc,100) 'z     ', mat(3,1),mat(3,2),mat(3,3) 
      write(ifc,100) 
      write(ifc,100) 'r     ', r(1),r(2),r(3) 
      write(ifc,100) 'theta ', theta(1),theta(2),theta(3) 
      write(ifc,100) 'phi   ', phi(1),phi(2),phi(3) 
      write(ifc,100) 
      write(ifc,120) 'Lambda', lambda(1),lambda(2),lambda(3) 
!                                                                       
  100 format(1x,a,3f15.4) 
  110 format(1x,a,3a15  ) 
  120 format(1x,a,3(1Pe15.4)) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!                                                                       
! BEGIN V2.2 compatibility                                              
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!     P. Robert, ScientiDev, Janvier 2019                               
!     subroutine for compatibility with previous V2.2 versions          
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
                                                                        
      subroutine cangrat(ux,uy,uz,vx,vy,vz,angle,ratio) 
      call cp_angle_and_ratio(ux,uy,uz,vx,vy,vz,angle,ratio) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cdatdoy(idoy,iyear,imonth,iday) 
      call cv_doty_to_date(idoy,iyear,imonth,iday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cdatj00(jd00,iyear,imonth,iday) 
      call cv_jul2000_to_date(jd00,iyear,imonth,iday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cdatj50(jd50,iyear,imonth,iday) 
      call cv_jul1950_to_date(jd50,iyear,imonth,iday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cdatwee(iweek,iyear,imonth,iday) 
      call cv_weekn_to_date(iweek,iyear,imonth,iday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cdipdir(iyear,idoy,d1,d2,d3) 
      call cp_geo_dipole_dir(iyear,idoy,d1,d2,d3) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cdoweek(iyear,imonth,iday,idow) 
      call cv_date_to_dotw(iyear,imonth,iday,idow) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cdoyear(iyear,imonth,iday,idoy) 
      call cv_date_to_doty(iyear,imonth,iday,idoy) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cfrdayn(iday,cday,nbcha) 
      character*(*) cday 
      call cp_fr_day_name(iday,cday,nbcha) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
                                                                        
      subroutine cfrmonn(imonth,cmonth,nchar) 
      character*(*) cmonth 
      call cp_fr_month_name(imonth,cmonth,nchar) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine chouday(ih,im,is,houday) 
      call cv_hms_to_dech(ih,im,is,houday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cjd1950(iyear,imonth,iday,jd50) 
      call cv_date_to_jul1950(iyear,imonth,iday,jd50) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cjd2000(iyear,imonth,iday,jd00) 
      call cv_date_to_jul2000(iyear,imonth,iday,jd00) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cmilday(ih,im,is,ims,milday) 
      call cv_dhms_to_msotd(ih,im,is,ims,milday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cnbdmon(iyear,imonth,nbday) 
      call cp_nbday_in_month(iyear,imonth,nbday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine coleapy(iyear,ileap) 
      call cp_leap_year(iyear,ileap) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine csundir(iyear,idoy,ih,im,is,gst,slong,sra,sdec,obliq) 
      call cp_gei_sun_dir(iyear,idoy,ih,im,is,gst,slong,sra,sdec,obliq) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine csunset(iyear,imon,iday,rlat,rlon,tmer,tris,tset,      &
     &             durd,elemer,azimer,eleris,aziris,eleset,aziset,icor) 
                                                                        
                                                                        
      common /sunVDH/ suntim(1441),sunele(1441),sunazi(1441) 
                                                                        
      character*(*) tmer,tris,tset,durd 
                                                                        
      call cp_sunrise_sunset(iyear,imon,iday,rlat,rlon,tmer,tris,tset,  &
     &             durd,elemer,azimer,eleris,aziris,eleset,aziset,icor) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine ctimhou(houday,ih,im,is) 
      call cv_dech_to_hms(houday,ih,im,is) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine ctimmil(milday,ih,im,is,ims) 
      call cv_msotd_to_hmsms(milday,ih,im,is,ims) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine ctimpa2(jd1950,houday) 
      call cp_time_param2(jd1950,houday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine ctimpa3(jd2000,houday) 
      call cp_time_param3(jd2000,houday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine ctimpar(iyear,imonth,iday,ih,im,is) 
      call cp_time_param(iyear,imonth,iday,ih,im,is) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cusdayn(iday,cday,nbcha) 
      character*(*) cday 
      call cp_en_day_name(iday,cday,nbcha) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cusmonn(imonth,cmonth,nchar) 
      character*(*) cmonth 
      call cp_en_month_name(imonth,cmonth,nchar) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine cweedoy(iyear,imonth,iday,iweek) 
      call cv_date_to_weekn(iyear,imonth,iday,iweek) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine gdipdir(dxgei,dygei,dzgei,dxgeo,dygeo,dzgeo) 
      call g_gei_geo_dipole_dir(dxgei,dygei,dzgei,dxgeo,dygeo,dzgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine gdiptan(diptan) 
      call g_gsm_dipole_tilt_angle(diptan) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine gecldir(exgei,eygei,ezgei,exgeo,eygeo,ezgeo) 
      call g_gei_geo_ecliptic_dir(exgei,eygei,ezgei,exgeo,eygeo,ezgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine gsrodir(rxgei,rygei,rzgei,rxgeo,rygeo,rzgeo) 
      call g_gei_geo_sun_rot(rxgei,rygei,rzgei,rxgeo,rygeo,rzgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine gsunpar(gmst,slon,sras,sdec,obli) 
      call g_gei_sun_param(gmst,slon,sras,sdec,obli) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine gvernum(vernum,verdat) 
      character*14 verdat 
      call g_rocot_version_number(vernum,verdat) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine plibinf 
      call print_rocot_info 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine recoor(x,y,z,cs) 
      character*1 cs 
      call r_coordinate_values(x,y,z,cs) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine recsys(csys) 
      character*(*) csys 
      call r_coordinate_system(csys) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine redate(iyear,imonth,iday) 
      call r_date(iyear,imonth,iday) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine retime(ih,im,is) 
      call r_time(ih,im,is) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tcarsph(x,y,z,r,teta,phi) 
      call t_car_to_sph(x,y,z,r,teta,phi) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tdmegeo(xdme,ydme,zdme,rlat,rlong,xgeo,ygeo,zgeo) 
      call t_dm_to_geo(xdme,ydme,zdme,rlat,rlong,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeigeo(xgei,ygei,zgei,xgeo,ygeo,zgeo) 
      call t_gei_to_geo(xgei,ygei,zgei,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeigse(xgei,ygei,zgei,xgse,ygse,zgse) 
      call t_gei_to_gse(xgei,ygei,zgei,xgse,ygse,zgse) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeigsm(xgei,ygei,zgei,xgsm,ygsm,zgsm) 
      call t_gei_to_gsm(xgei,ygei,zgei,xgsm,ygsm,zgsm) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeigsq(xgei,ygei,zgei,xgsq,ygsq,zgsq) 
      call t_gei_to_gseq(xgei,ygei,zgei,xgsq,ygsq,zgsq) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeimag(xgei,ygei,zgei,xmag,ymag,zmag) 
      call t_gei_to_mag(xgei,ygei,zgei,xmag,ymag,zmag) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeisma(xgei,ygei,zgei,xsma,ysma,zsma) 
      call t_gei_to_sm(xgei,ygei,zgei,xsma,ysma,zsma) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeodme(xgeo,ygeo,zgeo,rlat,rlong,xdme,ydme,zdme) 
      call t_geo_to_dm(xgeo,ygeo,zgeo,rlat,rlong,xdme,ydme,zdme) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeogei(xgeo,ygeo,zgeo,xgei,ygei,zgei) 
      call t_geo_to_gei(xgeo,ygeo,zgeo,xgei,ygei,zgei) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeogse(xgeo,ygeo,zgeo,xgse,ygse,zgse) 
      call t_geo_to_gse(xgeo,ygeo,zgeo,xgse,ygse,zgse) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeogsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm) 
      call t_geo_to_gsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeogsq(xgeo,ygeo,zgeo,xgsq,ygsq,zgsq) 
      call t_geo_to_gseq(xgeo,ygeo,zgeo,xgsq,ygsq,zgsq) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeomag(xgeo,ygeo,zgeo,xmag,ymag,zmag) 
      call t_geo_to_mag(xgeo,ygeo,zgeo,xmag,ymag,zmag) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeosma(xgeo,ygeo,zgeo,xsma,ysma,zsma) 
      call t_geo_to_sm(xgeo,ygeo,zgeo,xsma,ysma,zsma) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgeovdh(xgeo,ygeo,zgeo,rlat,rlong,xvdh,yvdh,zvdh) 
      call t_geo_to_vdh(xgeo,ygeo,zgeo,rlat,rlong,xvdh,yvdh,zvdh) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsegei(xgse,ygse,zgse,xgei,ygei,zgei) 
      call t_gse_to_gei(xgse,ygse,zgse,xgei,ygei,zgei) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsegeo(xgse,ygse,zgse,xgeo,ygeo,zgeo) 
      call t_gse_to_geo(xgse,ygse,zgse,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsegsm(xgse,ygse,zgse,xgsm,ygsm,zgsm) 
      call t_gse_to_gsm(xgse,ygse,zgse,xgsm,ygsm,zgsm) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsegsq(xgse,ygse,zgse,xgsq,ygsq,zgsq) 
      call t_gse_to_gseq(xgse,ygse,zgse,xgsq,ygsq,zgsq) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsesr2(xgse,ygse,zgse,rotx,roty,rotz,xsr2,ysr2,zsr2) 
      call t_gse_to_sr2(xgse,ygse,zgse,rotx,roty,rotz,xsr2,ysr2,zsr2) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsmgei(xgsm,ygsm,zgsm,xgei,ygei,zgei) 
      call t_gsm_to_gei(xgsm,ygsm,zgsm,xgei,ygei,zgei) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsmgeo(xgsm,ygsm,zgsm,xgeo,ygeo,zgeo) 
      call t_gsm_to_geo(xgsm,ygsm,zgsm,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsmgse(xgsm,ygsm,zgsm,xgse,ygse,zgse) 
      call t_gsm_to_gse(xgsm,ygsm,zgsm,xgse,ygse,zgse) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsmgsq(xgsm,ygsm,zgsm,xgsq,ygsq,zgsq) 
      call t_gsm_to_gseq(xgsm,ygsm,zgsm,xgsq,ygsq,zgsq) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsmmag(xgsm,ygsm,zgsm,xmag,ymag,zmag) 
      call t_gsm_to_mag(xgsm,ygsm,zgsm,xmag,ymag,zmag) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsmsma(xgsm,ygsm,zgsm,xsma,ysma,zsma) 
      call t_gsm_to_sm(xgsm,ygsm,zgsm,xsma,ysma,zsma) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsqgei(xgsq,ygsq,zgsq,xgei,ygei,zgei) 
      call t_gseq_to_gei(xgsq,ygsq,zgsq,xgei,ygei,zgei) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsqgeo(xgsq,ygsq,zgsq,xgeo,ygeo,zgeo) 
      call t_gseq_to_geo(xgsq,ygsq,zgsq,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsqgse(xgsq,ygsq,zgsq,xgse,ygse,zgse) 
      call t_gseq_to_gse(xgsq,ygsq,zgsq,xgse,ygse,zgse) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tgsqgsm(xgsq,ygsq,zgsq,xgsm,ygsm,zgsm) 
      call t_gseq_to_gsm(xgsq,ygsq,zgsq,xgsm,ygsm,zgsm) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tmaggei(xmag,ymag,zmag,xgei,ygei,zgei) 
      call t_mag_to_gei(xmag,ymag,zmag,xgei,ygei,zgei) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tmaggeo(xmag,ymag,zmag,xgeo,ygeo,zgeo) 
      call t_mag_to_geo(xmag,ymag,zmag,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tmaggsm(xmag,ymag,zmag,xgsm,ygsm,zgsm) 
      call t_mag_to_gsm(xmag,ymag,zmag,xgsm,ygsm,zgsm) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tmagsma(xmag,ymag,zmag,xsma,ysma,zsma) 
      call t_mag_to_sm(xmag,ymag,zmag,xsma,ysma,zsma) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsmagei(xsma,ysma,zsma,xgei,ygei,zgei) 
      call t_sm_to_gei(xsma,ysma,zsma,xgei,ygei,zgei) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsmageo(xsma,ysma,zsma,xgeo,ygeo,zgeo) 
      call t_sm_to_geo(xsma,ysma,zsma,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsmagsm(xsma,ysma,zsma,xgsm,ygsm,zgsm) 
      call t_sm_to_gsm(xsma,ysma,zsma,xgsm,ygsm,zgsm) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsmamag(xsma,ysma,zsma,xmag,ymag,zmag) 
      call t_sm_to_mag(xsma,ysma,zsma,xmag,ymag,zmag) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsphcar(r,teta,phi,x,y,z) 
      call t_sph_to_car(r,teta,phi,x,y,z) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsr2gse(xsr2,ysr2,zsr2,rotx,roty,rotz,xgse,ygse,zgse) 
      call t_sr2_to_gse(xsr2,ysr2,zsr2,rotx,roty,rotz,xgse,ygse,zgse) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsr2mfa(xsr2,ysr2,zsr2,bx,by,bz,rox,roy,roz,xm,ym,zm) 
      call t_sr2_to_mfa(xsr2,ysr2,zsr2,bx,by,bz,rox,roy,roz,xm,ym,zm) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsr2sre(xsr2,ysr2,spifre,spipha,deltaT,xsre,ysre) 
      call t_sr2_to_sr(xsr2,ysr2,spifre,spipha,deltaT,xsre,ysre) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tsresr2(xsre,ysre,spifre,spipha,deltaT,xsr2,ysr2) 
      call t_sr_to_sr2(xsre,ysre,spifre,spipha,deltaT,xsr2,ysr2) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
                                                                        
      subroutine tvdhgeo(xvdh,yvdh,zvdh,rlat,rlong,xgeo,ygeo,zgeo) 
      call t_vdh_to_geo(xvdh,yvdh,zvdh,rlat,rlong,xgeo,ygeo,zgeo) 
      return 
      END                                           
                                                                        
!     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
!                                                                       
! END ROCOTLIB 3.2                                                    
