c BEGIN ROCOTLIB 3.3
c
c XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c XX                                                                  XX
c XX   Robert's Coordinates Transformation Library  ROCOTLIB          XX
c XX                                                                  XX
c XX            Revised  version 3.3, November 2021                   XX
c XX                    ___________________                           XX
c XX                                                                  XX
c XX             Version 1.0 initially supported by                   XX
c XX                                                                  XX
c XX                  EUROPEAN  SPACE  AGENCY                         XX
c XX                                                                  XX
c XX               Study of the Cluster Mission                       XX
c XX                 Planning Related Aspects                         XX
c XX          within the Numerical Simulations Network                XX
c XX                                                                  XX
c XX            Patrick ROBERT, CRPE, November 1992                   XX
c XX                    ___________________                           XX
c XX                                                                  XX
c XX Further versions developped by the Author:                       XX
c XX                                                                  XX
c XX    version 1.0, November 1992                                    xx
c XX    version 1.1, July     1993                                    xx
c XX    version 1.2, January  1995                                    xx
c XX    version 1.3, July     2000 (Julian day 2000 et sun dir)       xx
c XX    version 1.4, June     2001 (for automatic documentation)      xx
c XX    version 1.5, December 2001 (adding cp_sunrise_sunset)         xx
c XX    version 1.6, Juin     2002 (upgrade IGRF -> 2005)             xx
c XX    version 1.7, December 2002 (Version for CDPP)                 xx
c XX    version 1.8, November 2003 (adding t_sr2_to_sr)               xx
c XX    version 1.9, March    2004 (compatibility with IDL)           xx
c XX    version 2.0, November 2006 (Update  IGRF -> 2010)             xx
c XX    version 2.1, November 2010 (Update  IGRF -> 2015)             xx
c XX    version 2.2, December 2011 (cp_sunrise_sunset on polar zone   xx
c XX    version 3.0, February 2016 (IGRF -> table of geomagnetic pole XX
c XX                                -> 2020 + some other coord. sys.) XX
c XX    version 3.1, January  2019 (add name compatibility with V2.2) XX
c XX    version 3.2, February 2020 (add trans. with Euler angles)     XX
c XX    version 3.3, November 2021 (update geo_dipole_dir 1900-2025)  XX
c XX                                                                  XX
c XX Copyright 1992, Patrick ROBERT, CNRS-ESA, All Rights reserved    XX
c XX                    ___________________                           XX
c XX                                                                  XX
c XX         For details, see the orginal document untitled:          XX
c XX CLUSTER Software Tools, part I: Coordinate Transformation LibraryXX
c XX             Document de travail DT/CRPE/1231                     XX
c XX             Patrick Robert, CRPE/TID, Juillet 1993               XX
c XX                                                                  XX
c XX         Available in CDPP:                                       XX
c XX             ROCOTLIB: a coordinate Transformation Library        XX
c XX             for Solar-Terrestrial studies                        XX
c XX             Patrick ROBERT, version 1.7 - January 2003,          XX
c XX             Rapport Interne no RI-CETP/02/2003                   XX
c XX                    ___________________                           XX
c XX                                                                  XX
c XX         The previous version 3.0 contains new transformations    XX
c XX         (TPN, MVA), and the cp_sunrise_sunset module allowing    XX
c XX         computation of sunrise and sunset anywhere on Earth,     XX
c XX         including polar zones.                                   XX
c XX                                                                  XX
c XX         The determination of the dipole axis direction           XX
c XX         has been replaced by a table available from years        XX
c XX         1900 to 2020, in order to have an easy way to code       XX
c XX         maintenance in the coming years.                         XX
c XX                                                                  XX
c XX         Operations on matrix (somme, product, inversion,         XX
c XX         computation of eigen vectors, diagonalization...)        XX
c XX         required for Minimum Variance Analysis coordinates       XX
c XX         have been added with prefix "mat_".                      XX
c XX                                                                  XX
c XX         The previous and temporary V3.1 is the same as 3.0       XX
c XX         but is compatible with V2.2 version concerning           XX
c XX         subroutine name (ex: tmaggsm -> t_mag_to_gsm)            XX
c XX                                                                  XX
c XX         The V3.2 version contains transformations  with Euler    XX
c XX         angles and their interpolation.                          XX
c XX                                                                  XX
c XX         Fortran 90, C, and IDL codes are also available.         XX
c XX                                                                  XX
c XX                                                                  XX
c XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
c
c     ***************************************************************0**
c
c     Reminder for the developper:
c     Any changes leadind to a new version implies following actions:
c        - Update the file header above
c        - Update the subroutine g_rocot_version_number
c        - Update the subroutine print_rocot_info
c        - Update the readme.txt file in the full package
c
c     ***************************************************************0**
c
c
      subroutine cp_angle_and_ratio(ux,uy,uz,vx,vy,vz,angle,ratio)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_angle_and_ratio beetween U and V vectors
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : ux,uy,uz
c              vx,vy,vz
c
c *   Output : sp=U.V
c              angle=angle beetween U and V (radians)
c              ratio= mod(U)/mod(V)
c ----------------------------------------------------------------------
c
c
      double precision u1,u2,u3,v1,v2,v3,dp,ru,rv,cot
c
      u1=dble(ux)
      u2=dble(uy)
      u3=dble(uz)
      v1=dble(vx)
      v2=dble(vy)
      v3=dble(vz)
c
      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
c
      angle=real(dacos(cot))
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_Euler_interp(a1,b1,c1,a2,b2,c2,ti,dt,ai,bi,ci)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_Euler_angles_interpolation
c *   Author : P. Robert, SDev, 2020
c
c *   Input  :  a1,b1,c1 Euler angles at time t1
c               a2,b2,c2 Euler angles at time t2
c               ti time for interpolation t1 < ti < t2
c               dt = t2 - t1
c
c *   Output :  ai,bi,ci Euler angles interpolated at time ti
c ----------------------------------------------------------------------
c

      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= Rmodulo(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 = Rmodulo(a2,360.)
      bi = Rmodulo(b2,360.)
      ci = Rmodulo(c2,360.)

      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_geo_dipole_dir(iyear,idoy,d1,d2,d3)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_dipole_direction in GEO system
c *   Author : P. Robert, LPP , 2016
c
c *   Input  :  iyear (1900 - 2020), idoy= day of year (1/1=1)
c *   Output :  d1,d2,d3  cartesian dipole components in GEO system
c ----------------------------------------------------------------------
c
      dimension kyear(132),rlat(132),rlon(132)
c
c *** table of geomagnetic Nort pole

c     issued from GEOPACK_08.f library
c     by   N.A. TSYGANENKO, DEC.1, 1991
c     revised JANUARY 01, 2020

c     from also Alken, P., Thébault, E., Beggan, C.D. et al. 
c     International Geomagnetic Reference Field: the thirteenth generation. 
c     Earth Planets Space 73, 49 (2021). 
c     https://doi.org/10.1186/s40623-020-01288-x
c     (warning: in this paper data are in WGS84 geodetic system)

      data kyear(  1),rlat(  1),rlon(  1) / 1900 , 78.68 , -68.79 /
      data kyear(  2),rlat(  2),rlon(  2) / 1901 , 78.68 , -68.78 /
      data kyear(  3),rlat(  3),rlon(  3) / 1902 , 78.68 , -68.77 /
      data kyear(  4),rlat(  4),rlon(  4) / 1903 , 78.68 , -68.77 /
      data kyear(  5),rlat(  5),rlon(  5) / 1904 , 78.68 , -68.76 /
      data kyear(  6),rlat(  6),rlon(  6) / 1905 , 78.68 , -68.75 /
      data kyear(  7),rlat(  7),rlon(  7) / 1906 , 78.68 , -68.74 /
      data kyear(  8),rlat(  8),rlon(  8) / 1907 , 78.67 , -68.74 /
      data kyear(  9),rlat(  9),rlon(  9) / 1908 , 78.67 , -68.73 /
      data kyear( 10),rlat( 10),rlon( 10) / 1909 , 78.66 , -68.73 /
      data kyear( 11),rlat( 11),rlon( 11) / 1910 , 78.66 , -68.72 /
      data kyear( 12),rlat( 12),rlon( 12) / 1911 , 78.66 , -68.69 /
      data kyear( 13),rlat( 13),rlon( 13) / 1912 , 78.65 , -68.66 /
      data kyear( 14),rlat( 14),rlon( 14) / 1913 , 78.65 , -68.63 /
      data kyear( 15),rlat( 15),rlon( 15) / 1914 , 78.64 , -68.60 /
      data kyear( 16),rlat( 16),rlon( 16) / 1915 , 78.64 , -68.57 /
      data kyear( 17),rlat( 17),rlon( 17) / 1916 , 78.64 , -68.53 /
      data kyear( 18),rlat( 18),rlon( 18) / 1917 , 78.64 , -68.49 /
      data kyear( 19),rlat( 19),rlon( 19) / 1918 , 78.63 , -68.46 /
      data kyear( 20),rlat( 20),rlon( 20) / 1919 , 78.63 , -68.42 /
      data kyear( 21),rlat( 21),rlon( 21) / 1920 , 78.63 , -68.38 /
      data kyear( 22),rlat( 22),rlon( 22) / 1921 , 78.63 , -68.36 /
      data kyear( 23),rlat( 23),rlon( 23) / 1922 , 78.63 , -68.34 /
      data kyear( 24),rlat( 24),rlon( 24) / 1923 , 78.62 , -68.31 /
      data kyear( 25),rlat( 25),rlon( 25) / 1924 , 78.62 , -68.29 /
      data kyear( 26),rlat( 26),rlon( 26) / 1925 , 78.62 , -68.27 /
      data kyear( 27),rlat( 27),rlon( 27) / 1926 , 78.62 , -68.27 /
      data kyear( 28),rlat( 28),rlon( 28) / 1927 , 78.61 , -68.27 /
      data kyear( 29),rlat( 29),rlon( 29) / 1928 , 78.61 , -68.26 /
      data kyear( 30),rlat( 30),rlon( 30) / 1929 , 78.60 , -68.26 /
      data kyear( 31),rlat( 31),rlon( 31) / 1930 , 78.60 , -68.26 /
      data kyear( 32),rlat( 32),rlon( 32) / 1931 , 78.59 , -68.28 /
      data kyear( 33),rlat( 33),rlon( 33) / 1932 , 78.59 , -68.30 /
      data kyear( 34),rlat( 34),rlon( 34) / 1933 , 78.58 , -68.32 /
      data kyear( 35),rlat( 35),rlon( 35) / 1934 , 78.58 , -68.34 /
      data kyear( 36),rlat( 36),rlon( 36) / 1935 , 78.57 , -68.36 /
      data kyear( 37),rlat( 37),rlon( 37) / 1936 , 78.57 , -68.39 /
      data kyear( 38),rlat( 38),rlon( 38) / 1937 , 78.56 , -68.42 /
      data kyear( 39),rlat( 39),rlon( 39) / 1938 , 78.56 , -68.45 /
      data kyear( 40),rlat( 40),rlon( 40) / 1939 , 78.55 , -68.48 /
      data kyear( 41),rlat( 41),rlon( 41) / 1940 , 78.55 , -68.51 /
      data kyear( 42),rlat( 42),rlon( 42) / 1941 , 78.55 , -68.51 /
      data kyear( 43),rlat( 43),rlon( 43) / 1942 , 78.55 , -68.52 /
      data kyear( 44),rlat( 44),rlon( 44) / 1943 , 78.55 , -68.52 /
      data kyear( 45),rlat( 45),rlon( 45) / 1944 , 78.55 , -68.53 /
      data kyear( 46),rlat( 46),rlon( 46) / 1945 , 78.55 , -68.53 /
      data kyear( 47),rlat( 47),rlon( 47) / 1946 , 78.55 , -68.59 /
      data kyear( 48),rlat( 48),rlon( 48) / 1947 , 78.55 , -68.66 /
      data kyear( 49),rlat( 49),rlon( 49) / 1948 , 78.55 , -68.72 /
      data kyear( 50),rlat( 50),rlon( 50) / 1949 , 78.55 , -68.79 /
      data kyear( 51),rlat( 51),rlon( 51) / 1950 , 78.55 , -68.85 /
      data kyear( 52),rlat( 52),rlon( 52) / 1951 , 78.55 , -68.91 /
      data kyear( 53),rlat( 53),rlon( 53) / 1952 , 78.55 , -68.97 /
      data kyear( 54),rlat( 54),rlon( 54) / 1953 , 78.54 , -69.04 /
      data kyear( 55),rlat( 55),rlon( 55) / 1954 , 78.54 , -69.10 /
      data kyear( 56),rlat( 56),rlon( 56) / 1955 , 78.54 , -69.16 /
      data kyear( 57),rlat( 57),rlon( 57) / 1956 , 78.55 , -69.22 /
      data kyear( 58),rlat( 58),rlon( 58) / 1957 , 78.56 , -69.28 /
      data kyear( 59),rlat( 59),rlon( 59) / 1958 , 78.56 , -69.35 /
      data kyear( 60),rlat( 60),rlon( 60) / 1959 , 78.57 , -69.41 /
      data kyear( 61),rlat( 61),rlon( 61) / 1960 , 78.58 , -69.47 /
      data kyear( 62),rlat( 62),rlon( 62) / 1961 , 78.58 , -69.55 /
      data kyear( 63),rlat( 63),rlon( 63) / 1962 , 78.59 , -69.62 /
      data kyear( 64),rlat( 64),rlon( 64) / 1963 , 78.59 , -69.70 /
      data kyear( 65),rlat( 65),rlon( 65) / 1964 , 78.60 , -69.77 /
      data kyear( 66),rlat( 66),rlon( 66) / 1965 , 78.60 , -69.85 /
      data kyear( 67),rlat( 67),rlon( 67) / 1966 , 78.61 , -69.92 /
      data kyear( 68),rlat( 68),rlon( 68) / 1967 , 78.62 , -69.98 /
      data kyear( 69),rlat( 69),rlon( 69) / 1968 , 78.64 , -70.05 /
      data kyear( 70),rlat( 70),rlon( 70) / 1969 , 78.65 , -70.11 /
      data kyear( 71),rlat( 71),rlon( 71) / 1970 , 78.66 , -70.18 /
      data kyear( 72),rlat( 72),rlon( 72) / 1965 , 78.53 , -69.85 /
      data kyear( 73),rlat( 73),rlon( 73) / 1966 , 78.55 , -69.92 /
      data kyear( 74),rlat( 74),rlon( 74) / 1967 , 78.56 , -69.98 /
      data kyear( 75),rlat( 75),rlon( 75) / 1968 , 78.57 , -70.05 /
      data kyear( 76),rlat( 76),rlon( 76) / 1969 , 78.58 , -70.11 /
      data kyear( 77),rlat( 77),rlon( 77) / 1970 , 78.59 , -70.18 /
      data kyear( 78),rlat( 78),rlon( 78) / 1971 , 78.61 , -70.24 /
      data kyear( 79),rlat( 79),rlon( 79) / 1972 , 78.63 , -70.29 /
      data kyear( 80),rlat( 80),rlon( 80) / 1973 , 78.65 , -70.35 /
      data kyear( 81),rlat( 81),rlon( 81) / 1974 , 78.67 , -70.41 /
      data kyear( 82),rlat( 82),rlon( 82) / 1975 , 78.69 , -70.47 /
      data kyear( 83),rlat( 83),rlon( 83) / 1976 , 78.71 , -70.53 /
      data kyear( 84),rlat( 84),rlon( 84) / 1977 , 78.73 , -70.58 /
      data kyear( 85),rlat( 85),rlon( 85) / 1978 , 78.76 , -70.64 /
      data kyear( 86),rlat( 86),rlon( 86) / 1979 , 78.78 , -70.70 /
      data kyear( 87),rlat( 87),rlon( 87) / 1980 , 78.81 , -70.76 /
      data kyear( 88),rlat( 88),rlon( 88) / 1981 , 78.84 , -70.79 /
      data kyear( 89),rlat( 89),rlon( 89) / 1982 , 78.87 , -70.81 /
      data kyear( 90),rlat( 90),rlon( 90) / 1983 , 78.91 , -70.84 /
      data kyear( 91),rlat( 91),rlon( 91) / 1984 , 78.94 , -70.87 /
      data kyear( 92),rlat( 92),rlon( 92) / 1985 , 78.97 , -70.90 /
      data kyear( 93),rlat( 93),rlon( 93) / 1986 , 79.01 , -70.94 /
      data kyear( 94),rlat( 94),rlon( 94) / 1987 , 79.04 , -70.99 /
      data kyear( 95),rlat( 95),rlon( 95) / 1988 , 79.07 , -71.03 /
      data kyear( 96),rlat( 96),rlon( 96) / 1989 , 79.11 , -71.08 /
      data kyear( 97),rlat( 97),rlon( 97) / 1990 , 79.14 , -71.13 /
      data kyear( 98),rlat( 98),rlon( 98) / 1991 , 79.18 , -71.18 /
      data kyear( 99),rlat( 99),rlon( 99) / 1992 , 79.21 , -71.24 /
      data kyear(100),rlat(100),rlon(100) / 1993 , 79.25 , -71.30 /
      data kyear(101),rlat(101),rlon(101) / 1994 , 79.29 , -71.36 /
      data kyear(102),rlat(102),rlon(102) / 1995 , 79.32 , -71.42 /
      data kyear(103),rlat(103),rlon(103) / 1996 , 79.37 , -71.45 /
      data kyear(104),rlat(104),rlon(104) / 1997 , 79.41 , -71.48 /
      data kyear(105),rlat(105),rlon(105) / 1998 , 79.46 , -71.51 /
      data kyear(106),rlat(106),rlon(106) / 1999 , 79.50 , -71.54 /
      data kyear(107),rlat(107),rlon(107) / 2000 , 79.54 , -71.57 /
      data kyear(108),rlat(108),rlon(108) / 2001 , 79.58 , -71.62 /
      data kyear(109),rlat(109),rlon(109) / 2002 , 79.63 , -71.66 /
      data kyear(110),rlat(110),rlon(110) / 2003 , 79.67 , -71.71 /
      data kyear(111),rlat(111),rlon(111) / 2004 , 79.71 , -71.76 /
      data kyear(112),rlat(112),rlon(112) / 2005 , 79.75 , -71.81 /
      data kyear(113),rlat(113),rlon(113) / 2006 , 79.80 , -71.88 /
      data kyear(114),rlat(114),rlon(114) / 2007 , 79.86 , -71.96 /
      data kyear(115),rlat(115),rlon(115) / 2008 , 79.91 , -72.05 /
      data kyear(116),rlat(116),rlon(116) / 2009 , 79.96 , -72.13 /
      data kyear(117),rlat(117),rlon(117) / 2010 , 80.02 , -72.21 /
      data kyear(118),rlat(118),rlon(118) / 2011 , 80.08 , -72.29 /
      data kyear(119),rlat(119),rlon(119) / 2012 , 80.13 , -72.37 /
      data kyear(120),rlat(120),rlon(120) / 2013 , 80.19 , -72.45 /
      data kyear(121),rlat(121),rlon(121) / 2014 , 80.25 , -72.53 /
      data kyear(122),rlat(122),rlon(122) / 2015 , 80.31 , -72.61 /
      data kyear(123),rlat(123),rlon(123) / 2016 , 80.37 , -72.63 /
      data kyear(124),rlat(124),rlon(124) / 2017 , 80.42 , -72.64 /
      data kyear(125),rlat(125),rlon(125) / 2018 , 80.48 , -72.65 /
      data kyear(126),rlat(126),rlon(126) / 2019 , 80.53 , -72.67 /
      data kyear(127),rlat(127),rlon(127) / 2020 , 80.59 , -72.68 /
      data kyear(128),rlat(128),rlon(128) / 2021 , 80.64 , -72.67 /
      data kyear(129),rlat(129),rlon(129) / 2022 , 80.69 , -72.66 /
      data kyear(130),rlat(130),rlon(130) / 2023 , 80.74 , -72.66 /
      data kyear(131),rlat(131),rlon(131) / 2024 , 80.79 , -72.65 /
      data kyear(132),rlat(132),rlon(132) / 2025 , 80.84 , -72.64 /
c
c                        ------------------------
c
      iy=iyear
      id=idoy
c
c *** Check date interval of validity
c
c     we are restricted by the interval 1900-2025, for which the igrf
c     coefficients are known;
c     if iyear is outside this interval, then the subroutine uses the
c     nearest limiting value and prints a warning:
c
      if(iy.le.1900      ) then
                     print 10, iy,1900
                     iy=1900
                     alat=rlat(1)
                     alon=rlon(1)
                     endif

      if(iy.ge.2025) then
                     print 10, iy,2025
                     iy=2025
                     alat=rlat(132)
                     alon=rlon(132)
                     endif
c
c *** decimal year
c
      year=float(iy) +float(id)/365.25
c
c *** find two known intervals

      do i=2,132
       if((iy .ge.kyear(i-1)) .and. (iy .lt. kyear(i))) then
         alat= rlat(i-1) + (rlat(i)-rlat(i-1))*(year-float(kyear(i-1)))
         alon= rlon(i-1) + (rlon(i)-rlon(i-1))*(year-float(kyear(i-1)))
         go to 20
       endif
      enddo

   20 continue

c
c *** direction of dipole axis in GEO system:
c
      clat=(90. -alat)*3.141593/180.
      clon=alon*3.141593/180.

      call t_sph_to_car(1.,clat,clon,d1,d2,d3)
c
      d1p=d1
      d2p=d2
      d3p=d3
c

   10 format(' * ROCOTLIB/cp_geo_dipole_dir: Warning! year=',i4.4,    
     &       ' [1900-2025] only, result for ',i4.4) 
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_gei_sun_dir(iyear,idoy,ih,im,is,
     &                          gst,slong,sra,sdec,obliq)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_sun_direction in GEI system
c *   Author : CT.Russel, CE-D, 1971, rev. P.R., 1992,2001,2002
c
c *   Comment: calculates four quantities in gei system necessary for
c              coordinate transformations dependent on sun position
c              (and, hence, on universal time and season)
c              Initial code from C.T. Russel, cosmic electro-dynamics,
c              v.2, 184-196, 1971.
c              accuracy: 0.006 degree.
c              Adaptation P.Robert, November 1992.
c              Revised and F90 compatibility, P. Robert June 2001.
c              Optimisation of DBLE computations and comments,
c              P. Robert, December 2002
c
c *   Input  : iyear : year (1901-2099)
c              idoy : day of the year (1 for january 1)
c              ih,im,is : hours, minutes, seconds U.T.
c
c *   Output : gst      greenwich mean sideral time (radians)
c              slong    longitude along ecliptic (radians)
c              sra      right ascension (radians)
c              sdec     declination of the sun (radians)
c              obliq    inclination of Earth's axis (radians)
c ----------------------------------------------------------------------
c
      double precision dj,fday
c
c
      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
c
      pi=acos(-1.)
      pisd= pi/180.
c
c *** Julian day and greenwich mean sideral time
c
      fday=dble(ih*3600+im*60+is)/86400.d0
      ileap=(iyear-1901)/4
c     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
c
c *** longitude along ecliptic
c
      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
c
c *** inclination of Earth's axis
c
      obliq=(23.45229 -0.0130125*t)*pisd
      sob=sin(obliq)
      cob=cos(obliq)
c
c     precession of declination (about 0.0056 deg., ajou P. Robert)
c
      pre= (0.0055686 - 0.025e-4*t)*pisd
c     (in TSY pre=9.924E-5)
c
c *** declination of the sun
c
      slp=slong -pre
      
c     the last constant is a correction for the angular aberration
c     due to earth's orbital motion    
      sind=sob*sin(slp)
      cosd=sqrt(1. -sind**2)
      sc=sind/cosd
      sdec=atan(sc)
c
c *** right ascension of the sun
c
      sra=pi -atan2((cob/sob)*sc, -cos(slp)/cosd)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      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)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_sunset_time and others
c *   Author : P. Robert, CRPE, 2001 Revised Dec. 2011
c *   Comment: Spherical Earth assumed
c
c *   Input  : iyear,imon,iday, rlat, rlong (rad)
c              icor: type de correction
c                  0 pas de correction
c                  1: correction de la refraction atmospherique seule
c                     (calcul comparable a l'IMC, ancien BDL)
c                  2: correction de la refraction et du bord superieur
c                     du soleil (calcul de la SAF).
c                  3: correction de la refraction avec l'elevation
c                     pour le passage au meridien (utile pour les zones
c                     polaires, ou le soleil peut raser l'horizon
c                     a midi)
c
c                  voir explications plus precises ci-dessous
c
c *   Output : tmer,tris,tset,durd on char*8 format, as '23:42:37'
c              tmer: Sun Meridian time (local noon)
c              tris: Sunrise time
c              tset: Sunset  time
c              durd: Duration of the day
c
c              elemer; elevation, en degres, du soleil au meridin
c              azimer: azimuth, en degres, a partir du nord vers l'est
c                      du soleil au meridien
c              eleris,aziris: meme chose pour le lever de soleil
c              eleset,aziset: meme chose pour le coucher
c
c              /sunVDH/ trajectoire du Soleil pendant la journée
c              Accuracy: 10 sec.
c
c    Les valeures pour Paris (latitude moyenne) et Tromso (zone polaire)
c    ont ete comparees aux ephemerides de l'institut de mecanique 
c    celeste (anciennement bureau des longitudes) avec icor=1
c    Les resultats sont comparables a la minute près pour les levers
c    et couchers, et a quelques secondes d'arc pour les elevations et 
c    azimuts, y compris dans les zones polaires et pour les transitions
c    avec les nuits polaires ou jours polaires (pas de lever ni coucher)
c
c ----------------------------------------------------------------------
c
c *** calcul des heures de lever et coucher de soleil:
c     correction du diametre apparent du soleil et de la refraction
c     a l'horizon. Les valeurs prises sont extraites de l'ouvrage
c     "introduction aux ephemerides astronomiques", publie par le
c     Bureau Des Longitudes.
c     Diametre apparent du soleil: 32'
c     Refraction a l'horizon     : 34' (36.6 pour le BDL)
c     correction= 32'/2 + 34' = 50'
c
c *** Attention :
c     Les Ephemrides Astronomiques du Bureau des Longitudes ne tiennent
c     compte que de la refraction a laquelle elles donnent la valeur
c     r = 36,6'.
c     Les Ephemrides Astronomiques de la S.A.F. considere le lever comme
c     l'apparition de son bord superieur.
c     Elles utilisent la valeur ht = -50'.
c
c *   variation avec la temperature et la pression:
c     on peut utiliser la formule de Bennet, pour corriger la refraction
c     en la multipliant par le coefficient:
c     (P/1010)(283/(273+T))
c     avec P en millibar et T en Celsius.
c     Note: Le BDL prend 1 pour cette correction, donc
c     suppose une presion P=1010 et une temperature T=10 deg.
c
c *   variation de R avec la hauteur apparente:
c     Une formule assez simple est donne par Bennet :
c     R = 1/Tan[h+7,31/(h+4,4)]
c     avec h en degre et R en minute d'arc
c     soit pour h=0 on trouve R= 34.5 
c          pour h=45          R=  0.99' 
c     (34' et 1' etant les valeurs couramment admises)
c
c *   Crepuscules : par definition la fin (le soir) ou le debut 
c     (le matin) des crepuscules civil, nautique et astronomique 
c     se produit quand le centre du Soleil est abaisse de 6, 12 et 18 
c     degres sous l'horizon. 
c ----------------------------------------------------------------------
c
      character*(*) tmer,tris,tset,durd
      dimension mimer(2), vemer(2)
c
      lmer=len(tmer)
      lris=len(tris)
      lset=len(tset)
      ldur=len(durd)
c
      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
c
      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
c
      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
c
      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
c
      milris= -1
      milset= -1
      milmer= -1
      milmin= -1
      
      imer=0
      vemer(1)= -2.
      vemer(2)= -2.
c
      sdp=0.
      svp=0.

      deg=180./acos(-1.)
c
c *   tolerance de 1 degre sur le meridien quand le soleil est 
c     au ras de l'horizon et donc sujet a la refraction
c
      cormax=1.
      svmin= -sin(cormax/deg)
      
c *** 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)

c *   correction de la refraction (voir + loin)

      call t_car_to_sph(sd,sh,sv,r,tetnc,phi)
c *   pas de correction de refraction par defaut
      tet=tetnc
c *   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
c
c     ----------------------------------------------
c     boucle sur le temps de la journee 
c     ----------------------------------------------

      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
c
c *** calcul de la direction du Soleil dans le VDH pour chaque pas de 
c     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)
c
c *** calcul des heures de lever et coucher de soleil:
c     correction du diametre apparent du soleil et de la refraction
c     a l'horizon. Les valeurs prises sont extraites de l'ouvrage
c     "introduction aux ephemerides astronomiques", publie par le
c     Bureau Des Longitudes.
c     Diametre apparent du soleil: 32'
c     Refraction a l'horizon     : 34'
c     correction= 32'/2 + 34' = 50'
c
c *** Attention :
c     Les Ephemrides Astronomiques du Bureau des Longitudes ne tiennent
c     compte que de la refraction a laquelle elles donnent la valeur
c     r = 36,6'.
c     Les Ephemrides Astronomiques de la S.A.F. considere le lever comme
c     l'apparition de son bord superieur.
c     Elles utilisent la valeur ht = -50'.
c
c
c *   dans le reperes DHV, teta est l'angle entre la verticale sortante 
c     et la direction du Soleil

      call t_car_to_sph(sd,sh,sv,r,tetnc,phi)

c *   pas de correction de refraction par defaut

      tet=tetnc

c *   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)

c *   calcul de la refraction a teta pour le calcul de la trajectoire
c     Quand le Soleil est vertical (tet=0), la correction est minime;
c     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

c *** chargement du common /sunVDH/
c     resolution en temps: 1mn 
c     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
c
c *** calcul du midi local, quand D (vers l'est) change de signe
c     attention: dans le cas du soleil de minuit, il y a deux meridien
c     ou le soleil est visible.
c     On prend celui ou le soleil est au plus haut 
c     (a midi, et non pas a minuit)
c
      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
c
c *** calcul du temps ou l'elevation est minimum (soleil de minuit)
c
      if(sv.ge.0. .and. svp.ge.0. .and. sv.lt.svp) milmin= isec*1000
c
c *** on a une transition (lever ou couche) quand V (vertical sortante)
c     change de signe (unite: millisec du jour)
c
c *   cas ou sv est nul (le soleil est a l'horizon)
c
      epsi=1.e-6
      
ccc   if(sv.eq.0.) then
      if(abs(sv).lt.epsi) then
                   sv=0.
                   if(svp.lt.0.) then
c                    juste avant, le soleil etait sous l'horizon
c                    la transition est donc un lever de soleil
                                 milris= isec*1000
                                 else
c                    juste avant, le soleil etait au dessus de l'horizon
c                    la transition est donc un coucher de soleil
                                 milset= isec*1000
                    endif
ccc   print*, '          cas 1 sv,svp=',sv,svp,' ris,set=',milris,milset
      go to 20
      endif

c *   cas ou svp est nul (le soleil etait a l'horizon)
c
      if(abs(svp).lt.epsi) then
                   if(sv.gt.0.) then
c                    le soleil est au dessus de l'horizon
c                    la transition est donc un lever de soleil
                                 milris= isec*1000
                                 else
c                    le soleil est au dessous de l'horizon
c                    la transition est donc un coucher de soleil
                                 milset= isec*1000
                    endif
ccc   print*, '          cas 2 sv,svp=',sv,svp,' ris,set=',milris,milset
      go to 20
      endif
c
c *   sv et svp sont non nuls
c
      if(sv*svp.lt.0.) then
c              si sv et svp sont de signe contraire, on a une transition
c              (lever ou coucher)
                       if(sv.gt.0.) then
c                    le soleil passe au passe au dessus de l'horizon
c                    la transition est donc un lever de soleil
                                    milris= isec*1000
                                    else
c                    le soleil passe au passe au dessous de l'horizon
c                    la transition est donc un coucher de soleil
                                    milset= isec*1000
                       endif
ccc   print*, '         cas 3 sv,svp=',sv,svp,' ris,set=',milris,milset     
      endif

      
   20 continue   
      sdp=sd
      svp=sv
      
    8 continue     
   10 continue
c
c     ----------------------------------------------
c     la boucle sur le temps de la journee est finie
c     ----------------------------------------------
      
c *   on restaure le jour en cours

      iyear=iyearp
      imon =imonp
      iday =idayp
      
c *** encodage des resultats sous la forme HH:MM:SS
c     et calcul de l'elevation et de l'azimuth pour chacun des 3 cas
c
      tmer='night   '
      tris='no SR   '
      tset='no SS   '
      durd='00:00:00'
      
      elemer=999.
      eleris=999.
      eleset=999.
      
      azimer=999.
      aziris=999.
      aziset=999.
c
c *   choix du midi local si soleil de minuit
c
      if(imer.gt.0) then
                    if(vemer(1).gt.vemer(2)) then
                                             milmer=mimer(1)
                                             else
                                             milmer=mimer(2)
                    endif
      endif
c
c *   calcul des positions au temps du meridien
c
      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)
c
            if(icor.eq.3) then
c                         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
c
c *   duree du jour calculee seuleument si le lever et le coucher
c     sont definis
c
      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
c
  100 format(i2.2,':',i2.2,':',i2.2)
  200 format(i2  ,':',i2.2,':',i2.2)
c
c *** cas du soleil de minuit: lever = coucher= temp elevation min
c
      if(durd.eq.'24:00:00' .and. milmin.gt.-1)  then
                      milris=milmin
                      milset=milmin
      endif
c
      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
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_time_param(iyear,imonth,iday,ih,im,is)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_time_parameters and time-dependent matrix
c *   Author : P. Robert, CRPE, 1992
c
c *   Comment: Prepare all time varying quantities for computations of
c              coordinate transforms of the library.
c
c *   Input  : iyear,imonth,iday (1901<year<2099)
c              ih,im,is : hours, minutes, seconds U.T.
c *   Output : in common statements
c ----------------------------------------------------------------------
c
      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
c
      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
c
      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
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      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
c
      data idoy /-1/
c
c *** day of year
c
      call cv_date_to_doty(iyear,imonth,iday,idoy)
c
c *** Dipole direction in GEO system
c
      call cp_geo_dipole_dir(iyear,idoy,gd1,gd2,gd3)
c
      qd12=sqrt(gd1**2 + gd2**2)
c
c *** Sun direction in GEI (from right ascension and declination),
c     greenwich mean sideral time and longitude along ecliptic
c
      call cp_gei_sun_dir(iyear,idoy,ih,im,is,
     &                    gst,slong,srasn,sdecl,obliq)
c
      gs1=cos(srasn)*cos(sdecl)
      gs2=sin(srasn)*cos(sdecl)
      gs3=sin(sdecl)
c
c *** sin and cos of GMST
c
      sgst=sin(gst)
      cgst=cos(gst)
c
c *** ecliptic pole in GEI system
c
      ge1=  0.
      ge2= -sin(obliq)
      ge3=  cos(obliq)
c
c *** direction of the rotation axis of the sun in GEI system
c     (from C.T. Russell, Cosmic Electro-dynamics, V.2, 184-196, 1971.)
c
      rad=  57.29578
c
      ras= -74.0/rad
      dec=  63.8/rad
c
      gr1=cos(ras)*cos(dec)
      gr2=sin(ras)*cos(dec)
      gr3=sin(dec)
c
c *** dipole direction in GEI system
c
      gm1= gd1*cgst - gd2*sgst
      gm2= gd1*sgst + gd2*cgst
      gm3= gd3
c
      qm12=sqrt(gm1**2 + gm2**2)
c
c *** direction of the sun in GEO system
c
      ps1=  gs1*cgst + gs2*sgst
      ps2= -gs1*sgst + gs2*cgst
      ps3=  gs3
c
c *** direction of the ecliptic in GEO system
c
      pe1=  ge1*cgst + ge2*sgst
      pe2= -ge1*sgst + ge2*cgst
      pe3=  ge3
c
c *** direction of the rotation axis of the sun in GEO system
c
      pr1=  gr1*cgst + gr2*sgst
      pr2= -gr1*sgst + gr2*cgst
      pr3=  gr3
c
c *** cross product MxS in GEI system
c
      gmgs1= gm2*gs3 - gm3*gs2
      gmgs2= gm3*gs1 - gm1*gs3
      gmgs3= gm1*gs2 - gm2*gs1
c
      rgmgs=sqrt(gmgs1**2 + gmgs2**2 + gmgs3**2)
c
c *** cross product ExS in GEI system
c
      gegs1= ge2*gs3 - ge3*gs2
      gegs2= ge3*gs1 - ge1*gs3
      gegs3= ge1*gs2 - ge2*gs1
c
c *** cross product RxS in GEI system
c
      grgs1= gr2*gs3 - gr3*gs2
      grgs2= gr3*gs1 - gr1*gs3
      grgs3= gr1*gs2 - gr2*gs1
c
      rgrgs=sqrt(grgs1**2 + grgs2**2 + grgs3**2)
c
c *** cross product RxE in GEI system
c
      grge1= gr2*ge3 - gr3*ge2
      grge2= gr3*ge1 - gr1*ge3
      grge3= gr1*ge2 - gr2*ge1
c
c *** cross product DxS in GEO system
c
      gdps1= gd2*ps3 - gd3*ps2
      gdps2= gd3*ps1 - gd1*ps3
      gdps3= gd1*ps2 - gd2*ps1
c
      rgdps=sqrt(gdps1**2 + gdps2**2 + gdps3**2)
c
c *** cross product ExS in GEO system
c
      peps1= pe2*ps3 - pe3*ps2
      peps2= pe3*ps1 - pe1*ps3
      peps3= pe1*ps2 - pe2*ps1
c
c *** cross product RxS in GEO system
c
      prps1= pr2*ps3 - pr3*ps2
      prps2= pr3*ps1 - pr1*ps3
      prps3= pr1*ps2 - pr2*ps1
c
      rprps=sqrt(prps1**2 + prps2**2 + prps3**2)
c
c *** computation of gei to mag vectors
c
      if(qm12.lt.1.e-30) stop '*** Rocotlib error qm12'
      xeima1=  gm1*gm3/qm12
      xeima2=  gm2*gm3/qm12
      xeima3= -qm12
c
      yeima1= -gm2/qm12
      yeima2=  gm1/qm12
      yeima3=  0.
c
c *** computation of gei to sm vectors
c
      if(rgmgs.lt.1.e-30) stop '*** Rocotlib error rgmgs'
      yeism1=gmgs1/rgmgs
      yeism2=gmgs2/rgmgs
      yeism3=gmgs3/rgmgs
c
      xeism1= yeism2*gm3 - yeism3*gm2
      xeism2= yeism3*gm1 - yeism1*gm3
      xeism3= yeism1*gm2 - yeism2*gm1
c
c *** computation of gei to gsm vectors
c
      yeigm1= gmgs1/rgmgs
      yeigm2= gmgs2/rgmgs
      yeigm3= gmgs3/rgmgs
c
      zeigm1= gs2*yeigm3 - gs3*yeigm2
      zeigm2= gs3*yeigm1 - gs1*yeigm3
      zeigm3= gs1*yeigm2 - gs2*yeigm1
c
c *** computation of gei to gseq vectors
c
      if(rgrgs.lt.1.e-30) stop '*** Rocotlib error rgrgs'
      yeigq1= grgs1/rgrgs
      yeigq2= grgs2/rgrgs
      yeigq3= grgs3/rgrgs
c
      zeigq1= gs2*yeigq3 - gs3*yeigq2
      zeigq2= gs3*yeigq1 - gs1*yeigq3
      zeigq3= gs1*yeigq2 - gs2*yeigq1
c
c *** computation of tetq angle
c
      stetq= (grge1*gs1 + grge2*gs2 + grge3*gs3)/rgrgs
      ctetq= sqrt(1.-stetq**2)
c
c *** computation of mu angle
c
      smu= ps1*gd1 + ps2*gd2 + ps3*gd3
      cmu= sqrt(1.-smu*smu)
c
c *** computation of dzeta angle
c
      cdze= (ge1*gm1   + ge2*gm2   + ge3*gm3)/rgmgs
      sdze= (ge1*gmgs1 + ge2*gmgs2 + ge3*gmgs3)/rgmgs
c     accuracy low on this angle
      epsi=1.e-4
      if(abs(sdze**2 +cdze**2-1.).gt.epsi) stop '*** Rocotlib error 3'
c
c *** computation of phi angle
c
      q=qd12*rgdps
      if(q.lt.1.e-30) stop '*** Rocotlib error q'
c
      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) stop '*** Rocotlib error 4'
c
c *** computation of geo to mag vectors
c
      if(qd12.lt.1.e-30) stop '*** Rocotlib error qd12'
      yeoma1= -gd2/qd12
      yeoma2=  gd1/qd12
      yeoma3=  0.
c
      xeoma1=  yeoma2*gd3
      xeoma2= -yeoma1*gd3
      xeoma3=  yeoma1*gd2 - yeoma2*gd1
c
c *** computation of geo to sm vectors
c
      if(rgdps.lt.1.e-30) stop '*** Rocotlib error rgdps'
      yeosm1= gdps1/rgdps
      yeosm2= gdps2/rgdps
      yeosm3= gdps3/rgdps
c
      xeosm1= yeosm2*gd3 - yeosm3*gd2
      xeosm2= yeosm3*gd1 - yeosm1*gd3
      xeosm3= yeosm1*gd2 - yeosm2*gd1
c
c *** computation of geo to gsm vectors
c
      yeogm1=yeosm1
      yeogm2=yeosm2
      yeogm3=yeosm3
c
      zeogm1= ps2*yeogm3 - ps3*yeogm2
      zeogm2= ps3*yeogm1 - ps1*yeogm3
      zeogm3= ps1*yeogm2 - ps2*yeogm1
c
c *** computation of geo to gsq vectors
c
      if(rprps.lt.1.e-30) stop '*** Rocotlib error rprps'
      yeogq1= prps1/rprps
      yeogq2= prps2/rprps
      yeogq3= prps3/rprps
c
      zeogq1= ps2*yeogq3 - ps3*yeogq2
      zeogq2= ps3*yeogq1 - ps1*yeogq3
      zeogq3= ps1*yeogq2 - ps2*yeogq1
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_time_param2(jd1950,houday)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_time_parameters and time-dependent matrix
c *   Author : P. Robert, CRPE, 2001
c
c *   Comment: Prepare all time varying quantities for computations of
c              coordinate transforms of the library.
c              Same as cp_time_param, only input arguments are changed.
c
c *   Input  : jd1950: Julian day 1950 (0= 1/1/1950)
c              houday: decimal hour of the day (U.T.)
c *   Output : in common statements
c ----------------------------------------------------------------------
c
      data iyear,imonth,iday /-1,-1,-1/
      data ih,im,is /-1,-1,-1/
c
      call cv_jul1950_to_date(jd1950,iyear,imonth,iday)
      call cv_dech_to_hms(houday,ih,im,is)
c
      call cp_time_param(iyear,imonth,iday,ih,im,is)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_time_param3(jd2000,houday)
c
c ----------------------------------------------------------------------
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_time_parameters and time-dependent matrix
c *   Author : P. Robert, CRPE, 2001
c
c *   Comment: Prepare all time varying quantities for computations of
c              coordinate transforms of the library.
c              Same as cp_time_param, only input arguments are changed.
c
c *   Input  : jd2000: Julian day 2000 (0= 1/1/2000)
c              houday: decimal hour of the day (U.T.)
c *   Output : in common statements
c ----------------------------------------------------------------------
c
      data iyear,imonth,iday /-1,-1,-1/
      data ih,im,is /-1,-1,-1/
c
      call cv_jul1950_to_date(jd2000,iyear,imonth,iday)
      call cv_dech_to_hms(houday,ih,im,is)
c
      call cp_time_param(iyear,imonth,iday,ih,im,is)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_tpn_param(xo,yo,zo,xs, Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz)
c
      real Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz
      
c
c     ---------------------------------------------------------------+--
c *   Class  : basic compute modules of Rocotlib Software
c *   Object : compute_TPN_system
c *   Author : P. Robert, CETP, 2004
c *   Comment: Compute TPN vector in GSE system or any system having X
c              axis towards the SUN.
c
c         N: Output normal to the paraboloid
c         T: tengente to the paraboloid, towards the summit
c         P: tengente to the paraboloid, P=N X T
c
c     The paraboloid is defined by its summit Xs and the local point 
c     Xo, Yo, Zo
c
c     Note: the paraboloid is close to the magnetopause if the summit
c           is defined as the subsolar point (by T87,T89 model or other)
c           and if the local point Xo,Yo,Zo correspond to magnetopause 
c           crossing.
c
c *   Input : xo,yo,zo,xs
c *   Output: Nx,Ny,Nz,Tx,Ty,Tz,Px,Py,Pz
c     ---------------------------------------------------------------+--
c      
c
      r= sqrt(yo**2 +zo**2)
      dx=xs-xo
c
      Nx= r
      Ny= 2.*dx*yo/r
      Nz= 2.*dx*zo/r
c
      Tx= 2.*dx
      Ty= -yo
      Tz= -zo
c
      Px= Ny*Tz -Nz*Ty
      Py= Nz*Tx -Nx*Tz
      Pz= Nx*Ty -Ny*Tx
c
c     normalisation
c
      rnoN= sqrt(Nx**2 +Ny**2 +Nz**2)
      rnoT= sqrt(Tx**2 +Ty**2 +Tz**2)
      rnoP= sqrt(Px**2 +Py**2 +Pz**2)
c
      Nx=Nx/rnoN
      Ny=Ny/rnoN
      Nz=Nz/rnoN
c
      Tx=Tx/rnoT
      Ty=Ty/rnoT
      Tz=Tz/rnoT
c
      Px=Px/rnoP
      Py=Py/rnoP
      Pz=Pz/rnoP
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_nbday_in_month(iyear,imonth,nbday)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : compute_number_of_day_of_the_month
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : iyear,imonth (1-12)
c *   Output : nbday
c ----------------------------------------------------------------------
c
c
      dimension nday(12)
      data nday/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
c
      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
c
      call cp_leap_year(iyear,ily)
c
      if(ily.eq.1) then
                   nday(2)=29
                   else
                   nday(2)=28
                   endif
c
      nbday=nday(imonth)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_en_day_name(iday,cday,nbcha)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : compute_english_day_name, ex: 'Monday' for iday=1
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : iday (1-7, otherwise modulo is done)
c *   Output : cday,nchar
c ----------------------------------------------------------------------
c
      character*(*) cday
      character*9 days (7)
      integer nbca(7)
c
      data days /'Monday','Tuesday','Wednesday','Thurday','Friday',
     &           'Saturday','Sunday'/
      data nbca /6,7,9,7,6,8,6/
c
      iweek= iday/7
      nday= iday -7*iweek
      if(nday.lt.1) nday=nday+7
c
      cday= days(nday)
      nbcha=nbca(nday)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_en_month_name(imonth,cmonth,nchar)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : compute_english_month_name
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : imonth (1-12)
c *   Output : cmonth,nchar
c ----------------------------------------------------------------------
c
c
      integer nbca(12)
      character*9 cara(12)
      character*(*) cmonth
c
      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/
c
c
      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
c
      cmonth=cara(imonth)
      nchar =nbca(imonth)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_fr_day_name(iday,cday,nbcha)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : compute_french_day_name, ex: 'Lundi' for iday=1
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : iday (1-7, otherwise modulo is done)
c *   Output : cday,nchar
c ----------------------------------------------------------------------
c
      character*(*) cday
      character*8 days (7)
      integer nbca(7)
c
      data days /'Lundi','Mardi','Mercredi','Jeudi','Vendredi',
     &           'Samedi','Dimanche'/
      data nbca /5,5,8,5,8,6,8/
c
      iweek= iday/7
      nday= iday -7*iweek
      if(nday.lt.1) nday=nday+7
c
      cday= days(nday)
      nbcha=nbca(nday)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_fr_month_name(imonth,cmonth,nchar)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : compute_french_month_name
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : imonth (1-12)
c *   Output : cmonth,nchar
c ----------------------------------------------------------------------
c
c
      integer nbca(12)
      character*9 cara(12)
      character*(*) cmonth
c
      data cara /'Janvier','Février','Mars','Avril','Mai','Juin',
     &           'Juillet','Août','Septembre','Octobre','Novembre',
     &           'Décembre'/
c
      data nbca /7,7,4,5,3,4,7,4,9,7,8,8/
c
c
      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
c
      cmonth=cara(imonth)
      nchar =nbca(imonth)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_leap_year(iyear,ileap)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : compute_leap_year with ileap=1 for leap year, 0 if not
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : iyear (ex: 1980)
c *   Output : ileap (1 or 0 if iyear is or not a leap year)
c ----------------------------------------------------------------------
c
      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
c
      ir=iyear-(iyear/4)*4
      if(ir.eq.0) then
                  ileap=1
                  else
                  ileap=0
                  endif
c
      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
c
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cp_seasons(iyear,id_sso,id_wso,id_seq,id_feq,
     &                            ct_sso,ct_wso,ct_seq,ct_feq)
c
      character*5 ct_sso,ct_wso,ct_seq,ct_feq
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : compute_seasons, i.e. solstice & equinox
c *   Author : P. Robert, SDev, 2017
c
c *   Input  : iyear (ex: 1980)
c *   Output : id_sso,id_wso : June and december day of summer and
c                              winter solstice.
c              id_seq,id_feq : same for march and september spring
c                              and fall equinoxes.
c              ct_sso,ct_wso : June and december time of summer and
c                              winter solstice.
c              ct_seq,ct_feq : same for march and september spring
c                             and fall equinoxes; Cha*5 (hh:mm)
c *
c *   Comment: calcul des saisons: solstices & equinoxes, 
c *            precision 1 heure environ 
c ----------------------------------------------------------------------    

      id_sso=0
      id_wso=0
      id_seq=0
      id_feq=0

      ct_sso='?'
      ct_wso='?'
      ct_seq='?'
      ct_feq='?'
c
c *** calcul de l'equinoxe de printemps
c
      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
c
c *** calcul de l'equinoxe d'automne
c
      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
c
c *** calcul du solstice d'ete
c
      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
c
c *** calcul du solstice d'hiver
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_doty_to_date(idoy,iyear,imonth,iday)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_day_of_year for a given year in date
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : iyear,idoy (idoy=1 for january 1)
c *   Output : imonth,iday
c ----------------------------------------------------------------------
c
c
      dimension month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
c
      call cp_leap_year(iyear,ily)
c
      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
c
      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
c
      if(ily.eq.1) then
                   month(2)=29
                   else
                   month(2)=28
                   endif
c
      m=0
c
      do 10 im=1,12
      mp=m
      m=m+month(im)
      if(idoy.le.m) go to 20
   10 continue
   20 continue
c
      imonth=im
      iday=idoy-mp
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_jul2000_to_date(jd00,iyear,imonth,iday)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_julian_day_2000 in date
c *   Author : P. Robert, CRPE, 1992
c *   Comment: compute date as year, month, day from julian day 2000
c
c *   Input  : jd00 julian day 2000 (0= 1/1/2000)
c *   Output : iyear,imonth,iday
c ----------------------------------------------------------------------
c
c
      jd50 = jd00 +18262
c
      call cv_jul1950_to_date(jd50,iyear,imonth,iday)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_jul1950_to_date(jd50,iyear,imonth,iday)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_julian_day_1950 in date
c *   Author : P. Robert, CRPE, 1992
c *   Comment: compute date as year, month, day from julian day 1950
c
c *   Input  : jd50  julian day 1950 (0= 1/1/1950)
c *   Output : iyear,imonth,iday
c ----------------------------------------------------------------------
c
      data ily /-1/
c
c
      jd1= -1
c
      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
c
   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'
c
   20 continue
      jd=jd50-jdp
      call cv_doty_to_date(jd,iy,imonth,iday)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_weekn_to_date(iweek,iyear,imonth,iday)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_first_day_of_week_number in date
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : iweek,iyear
c *   Output : imonth,iday
c ----------------------------------------------------------------------
c
c
c  calcul de la date correspondant au premier jour de la semaine
c
c *** premier lundi de l'annee
c
      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
c
c *** semaine correspondant au premier lundi
c
      call cv_date_to_weekn(iyear,1,ipl,iweek1)
c
c *** jour de l'annee correspondant au lundi de la semaine
c
      idoy= ipl +(iweek-iweek1)*7

      if(idoy.lt.1) then
                    iyear=iyear-1
                    imonth=12
                    iday=31+idoy
                    return
                    endif
c
c *** date correspondante
c
      call cv_doty_to_date(idoy,iyear,imonth,iday)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_date_to_dotw(iyear,imonth,iday,idow)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_date in day_of_the_week
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : iyear,imonth,iday
c *   Output : idow (1-7)
c ----------------------------------------------------------------------
c
      data jdref,julday /-100000,-100000/
c
c
c *** jour julien du lundi de reference (le 1/1/2001 etait un lundi)
c
      call cv_date_to_jul2000(2001,1,1,jdref)
c
c *** jour julien de la date demandee
c
      call cv_date_to_jul2000(iyear,imonth,iday,julday)
c
c *** calcul du jour de la semaine
c
      idiff= julday-jdref
c
      iweek= idiff/7
      idow= idiff -7*iweek +1
      if(idow.lt.1) idow=idow+7
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_date_to_doty(iyear,imonth,iday,idoy)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_date in day_of_year with idoy=1 for january 1
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : iyear,imonth,iday   ex: 1990,10,17
c *   Output : idoy                ex: 290
c ----------------------------------------------------------------------
c
c
      dimension month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
c
      call cp_leap_year(iyear,ily)
c
      if(ily.eq.1) then
                   month(2)=29
                   else
                   month(2)=28
                   endif
c
      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
c
      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
c
      idoy=iday
      do 10 i=1,imonth-1
      idoy=idoy+month(i)
   10 continue
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_hms_to_dech(ih,im,is,houday)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_hours_minutes_seconds in decimal hour of the day
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : ih,im,is
c *   Output : houday  decimal hour of the day
c ----------------------------------------------------------------------
c
      houday=float(ih)+float(im)/60.+float(is)/3600.
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_date_to_jul1950(iyear,imonth,iday,jd50)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_date in julian_day_1950 with jd50=0 for jan 1
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : iyear,imonth,iday   ex: 1990,10,17
c *   Output : jd50
c ----------------------------------------------------------------------
c
      data ily,idoy /-1,-1/
c
c
      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
c
      call cv_date_to_doty(iyear,imonth,iday,idoy)
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_date_to_jul2000(iyear,imonth,iday,jd00)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_date in julian_day_2000 with jd00=0 for january 1
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : iyear,imonth,iday   ex: 2001,10,17
c *   Output : jd00 (may be negative)
c ----------------------------------------------------------------------
c
      data jd50 /-1/
c
c
      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
c
      call cv_date_to_jul1950(iyear,imonth,iday,jd50)
c
      jd00= jd50 - 18262
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_dhms_to_msotd(ih,im,is,ims,milday)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_hours_minutes_seconds_ms in millisec_of_day
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : ih,im,is,ims
c *   Output : milday  millisecond of the day
c ----------------------------------------------------------------------
c
      milday= ih*3600000 + im*60000 + is*1000 +ims
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_dech_to_hms(houday,ih,im,is)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_decimal hour of the day in time
c *   Author : P. Robert, CRPE, 1992
c
c *   Input  : houday decimal hour of the day
c *   Output : ih,im,is
c ----------------------------------------------------------------------
c
c
      nj=int(abs(houday)/24.)
      rdech=abs(houday)-float(nj*24)
c
      ih=int(rdech)
      im=int((rdech-float(ih))*60.)
      is=int((rdech-float(ih)-float(im)/60.)*3600. +0.001)
c
      if(is.eq.60) then
                   is=0
                   im=im +1
                   endif
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_msotd_to_hmsms(milday,ih,im,is,ims)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_millisec. of the day in time
c *   Author : P. Robert, CRPE, 2001
c
c *   Input  : milday millisecond of the day
c *   Output : ih,im,is,ims
c ----------------------------------------------------------------------
c
c
      nj= milday/86400000
      mil2= milday-nj*86400000
c
      if(nj.ne.0) then
          print*, '*** Rocotlib/cv_msotd_to_hmsms: milday=',milday,
     &                                           'is > 86400000'
          print*, '                      assumed:',mil2
                  endif
c
      ih=  milday/3600000
      im= (milday-ih*3600000)/60000
      is= (milday-ih*3600000-im*60000)/1000
      ims= milday-ih*3600000 -im*60000 -is*1000
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cv_date_to_weekn(iyear,imonth,iday,iweek)
c
c ----------------------------------------------------------------------
c *   Class  : calendar modules of Rocotlib Software
c *   Object : convert_date in week_of_the_year
c *   Author : P. Robert, CRPE, 2001, revised SD April 2017
c
c *   Input  : iyear,imonth,iday
c *   Output : iweek
c ----------------------------------------------------------------------
c
c     la semaine 1 de toute année est celle qui contient le 4 janvier 
c     ou celle qui contient le 1er jeudi de janvier(norme ISO-8601). 

c *** calcul du jour de la semaine du 1er janvier de l'annee

      call cv_date_to_dotw(iyear,1,1,idow)

c *** 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
c
c *** jour de l'annee de la date demandee
c
      call cv_date_to_doty(iyear,imonth,iday,idoty)
c
c *** jour de la semaine de la date demandee
c
      call cv_date_to_dotw(iyear,imonth,iday,idow)
c
c *** jour de l'annee du lundi de la semaine de la date demandee
c
      idotL= idoty -idow +1

c *** numero de la semaine

      iweek= (idotL -idoty2)/7 +2
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine g_gei_geo_dipole_dir(dxgei,dygei,dzgei,
     &                                dxgeo,dygeo,dzgeo)
c
c ----------------------------------------------------------------------
c *   Class  : give modules of Rocotlib Software
c *   Object : give_dipole_direction in GEI and GEO system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: values are extracted from common
c
c *   Input  : none
c *   Output : dxgei,dygei,dzgei  cartesian dipole GEI coordinates
c              dxgeo,dygeo,dzgeo  cartesian dipole GEO coordinates
c ----------------------------------------------------------------------
c
      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
c
      dxgei=gm1
      dygei=gm2
      dzgei=gm3
c
      dxgeo=gd1
      dygeo=gd2
      dzgeo=gd3
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine g_gsm_dipole_tilt_angle(diptan)
c
c ----------------------------------------------------------------------
c *   Class  : give modules of Rocotlib Software
c *   Object : give_dipole_tilt_angle in radians
c *   Author : P. Robert, CRPE, 1992
c *   Comment: values are extracted from common
c
c *   Input  : none
c *   Output : diptan (radians)
c ----------------------------------------------------------------------
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      diptan=atan2(smu,cmu)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine g_gei_geo_ecliptic_dir(exgei,eygei,ezgei,
     &                                  exgeo,eygeo,ezgeo)
c
c ----------------------------------------------------------------------
c *   Class  : give modules of Rocotlib Software
c *   Object : give_ecliptic_direction in GEI and GEO system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: values are extracted from common
c
c *   Input  : none
c *   Output : exgei,eygei,ezgei  cartesian ecliptic GEI coordinates
c              exgeo,eygeo,ezgeo  cartesian ecliptic GEO coordinates
c ----------------------------------------------------------------------
c
      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
c
      exgei=ge1
      eygei=ge2
      ezgei=ge3
c
      exgeo=pe1
      eygeo=pe2
      ezgeo=pe3
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine g_gei_geo_sun_rot(rxgei,rygei,rzgei,rxgeo,rygeo,rzgeo)
c
c ----------------------------------------------------------------------
c *   Class  : give modules of Rocotlib Software
c *   Object : give_sun_rotation_direction in GEI and GEO system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: values are extracted from common
c
c *   Input  : none
c *   Output : rxgei,rygei,rzgei cartesian sun rotation GEI coordinates
c              rxgeo,rygeo,rzgeo cartesian sun rotation GEO coordinates
c ----------------------------------------------------------------------
c
      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
c
      rxgei=gr1
      rygei=gr2
      rzgei=gr3
c
      rxgeo=pr1
      rygeo=pr2
      rzgeo=pr3
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
c
c ----------------------------------------------------------------------
c *   Class  : give modules of Rocotlib Software
c *   Object : give_sun_direction in GEI and GEO system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: values are extracted from common
c
c *   Input  : none
c *   Output : sxgei,sygei,szgei  cartesian sun GEI coordinates
c              sxgeo,sygeo,szgeo  cartesian sun GEO coordinates
c ----------------------------------------------------------------------
c
      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
c
      sxgei=gs1
      sygei=gs2
      szgei=gs3
c
      sxgeo=ps1
      sygeo=ps2
      szgeo=ps3
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine g_gei_sun_param(gmst,slon,sras,sdec,obli)
c
c ----------------------------------------------------------------------
c *   Class  : give modules of Rocotlib Software
c *   Object : give_sun_parameter dependant of time in GEI system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: values are extracted from common
c
c *   Input  : none
c *   Output : gmst   greenwich mean sideral time (radians)
c              slon   longitude along ecliptic (radians)
c              sras   right ascension (radians)
c              sdec   declination of the sun (radians)
c ----------------------------------------------------------------------
c
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq
c
      gmst= gst
      slon= slong
      sras= srasn
      sdec= sdecl
      obli= obliq
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine g_rocot_version_number(vernum,verdat)
c
      character*14 verdat
c
c ----------------------------------------------------------------------
c *   Class  : give modules of Rocotlib Software
c *   Object : give_version_number and modification date of the library
c *   Author : P. Robert, CRPE, 1992
c *   Comment: values are extracted from common
c
c *   Input  : none
c *   Output : vernum (ex 1.0) and verdat (ex: 'January 1995')
c ----------------------------------------------------------------------
c
c     For previous versions informations, see print_rocot_info subrout.
c
      vernum=3.3
      verdat='Nov 2021'
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine print_rocot_info
c
      character*14 verdat
c
c ----------------------------------------------------------------------
c *   Class  : print modules of Rocotlib Software
c *   Object : print_library_informations
c *   Author : P. Robert, CRPE, 1992
c *   Comment: could be a main program
c
c *   Input  : none
c *   Output : none; print infos on output
c ----------------------------------------------------------------------
c
      call g_rocot_version_number(vernum,verdat)
c
c
      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, ' version 3.3, November 2021 (update geo_dipole_dir )'
      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 V3.2 version contains transformations with Euler'
      print 100, ' 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, ' '
c
  100 format(3a)
  200 format(a,f4.1,2a)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine r_coordinate_values(x,y,z,cs)
c
c ----------------------------------------------------------------------
c *   Class  : read modules of Rocotlib Software
c *   Object : read coordinate values from input
c *   Author : P. Robert, CRPE, 2002
c
c *   Comment: read cs and x,y,z cartesian  or spherical coordinates
c *            print error if cs is not valid, and ask again
c
c *   Input  : cs (c or s) and x,y,z on standard input
c *   Output : cs and x,y,z always in cartesian coordinates
c ----------------------------------------------------------------------
c
      character*1 cs
c
c
      pisd=acos(-1.)/180.
c
   10 continue
c
      print 100, 'input coordinates: cartesian (c) or spherical (s)'
      read *, cs
      print 100, cs
c
      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
c
                    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)'
c
                    else
                    print*, 'x,y,z ? (RE) '
                    read *,  x,y,z
                    print*,  x,y,z
c
                    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
c
  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)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine r_coordinate_system(csys)
c
c ----------------------------------------------------------------------
c *   Class  : read modules of Rocotlib Software
c *   Object : read coordinate system from input and check validity
c *   Author : P. Robert, CRPE, 2002
c
c *   Comment: read csys string variable and check validity
c *            (only gei, geo, mag, sma, gsm, gse, gsq)
c *            print error if csys is not valid, and ask again
c
c *   Input  : csys on standard input
c *   Output : csys
c ----------------------------------------------------------------------
c
      character*3 csys
c
   10 continue
c
      print 100, 'repere ? (gei, geo, mag, sma, gsm, gse, gsq)'
      read *,csys
      print 100,csys
  100 format(a)
c
      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
c
      print*, 'Only gei, geo, mag, sma, gsm, gse or gsq please...'
      print*, 'again...'
      go to 10
c
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine r_date(iyear,imonth,iday)
c
c ----------------------------------------------------------------------
c *   Class  : read modules of Rocotlib Software
c *   Object : read_date from input and check validity
c *   Author : P. Robert, CRPE, 1992
c
c *   Comment: test if year is gt 1900
c *            test if imonth is not greater than 12
c              test if iday is not greather than lengh of the month,
c                        takink into account the leap years.
c               print error if date is not valid, and ask again
c
c *   Input  : iyear,imonth,iday given in standard input
c *   Output : iyear,imonth,iday
c ----------------------------------------------------------------------
c
c
      dimension month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
c
   10 continue
      print 100, ' iyear,imonth,iday ? (ex: 1990,10,17)'
      read *,  iyear,imonth,iday
      print 110,  iyear,imonth,iday
c
  100 format(a)
  110 format(1x,i4,1x,i2,1x,i2)
c
      if(iyear.lt.1900) then
         print*, '*** Rocotlib/r_date: iyear must be greater than 1900'
         print*, '                     again ...'
                        go to 10
                        endif
c
      call cp_leap_year(iyear,ily)
c
c
      if(ily.eq.1) then
                   month(2)=29
                   else
                   month(2)=28
                   endif
c
      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
c
      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
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine r_time(ih,im,is)
c
c ----------------------------------------------------------------------
c *   Class  : read modules of Rocotlib Software
c *   Object : read_time from input and check validity
c *   Author : P. Robert, CRPE, 1992
c *   Comment: read hour, minute and second and verifie validity
c              ih must be between 1 and 23, im and is between 1 and 59
c              print error if time is not valid, and ask again
c
c *   Input  : ih,im,is on standard input
c *   Output : ih,im,is
c ----------------------------------------------------------------------
c
   10 continue
      print 100, ' hour, minute, second ? (ex: 10,45,50)'
      read *, ih,im,is
      print 110, ih,im,is
c
  100 format(a)
  110 format(1x,i2,1x,i2,1x,i2)
c
      if(ih.lt.0.or.ih.gt.23) then
                              print*, 'hour between 0 and 23 please'
                              print*, 'again...'
                              go to 10
                              endif
c
      if(im.lt.0.or.im.gt.59) then
                              print*, 'minute between 0 and 59 please'
                              print*, 'again...'
                              go to 10
                              endif
c
      if(is.lt.0.or.is.gt.59) then
                              print*, 'second between 0 and 59 please'
                              print*, 'again...'
                              go to 10
                              endif
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_car_to_sph(x,y,z,r,teta,phi)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_car_to_sph: CAR -> SPH  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: none
c
c *   Input  :   x,y,z        cartesian coordinates
c *   Output :   r,teta,phi   spherical coordinates (angles in radians)
c ----------------------------------------------------------------------
c
c
      teta=0.
      phi=0.
      sq=x**2+y**2
      r=sqrt(sq+z**2)
      pi=acos(-1.)
      pisd=pi/180.
c
      if(r.lt.1.e-30) return
c
c *** en dessous de 1/10000 degres, on considere que teta =0 ou 180
c     et phi indefini, mis a zero
c
      if(sq/r.gt.1.7e-6) then
              phi=atan2(y,x)
              if(phi.lt.-179.999*pisd.and.phi.gt.-180.0001*pisd) phi=pi
              if(phi.lt. 0.0001*pisd.and.phi.gt.-0.0001*pisd) phi=0.
              teta=acos(z/r)
              return
      endif
c
      if (z.lt.0.) then
                   teta=acos(-1.)
                   else
                   teta=0.
      endif
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_dm_to_geo(xdme,ydme,zdme,rlat,rlong,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_dme_to_geo: DM  -> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xdme,ydme,zdme   cartesian dme coordinates
c              rlat,rlong       latitude and longitude of the point
c                               of observation (radians)
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
c
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
c
      y1= real(dprod(gd2,r3) - dprod(gd3,r2))
      y2= real(dprod(gd3,r1) - dprod(gd1,r3))
      y3= real(dprod(gd1,r2) - dprod(gd2,r1))
c
      q=sqrt(y1**2 + y2**2 + y3**2)
c
      y1=y1/q
      y2=y2/q
      y3=y3/q
c
      x1= real(dprod(y2,gd3) - dprod(y3,gd2))
      x2= real(dprod(y3,gd1) - dprod(y1,gd3))
      x3= real(dprod(y1,gd2) - dprod(y2,gd1))
c
c
      xgeo= x1*xdme + y1*ydme + gd1*zdme
      ygeo= x2*xdme + y2*ydme + gd2*zdme
      zgeo= x3*xdme + y3*ydme + gd3*zdme
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gei_to_geo(xgei,ygei,zgei,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gei_to_geo: GEI -> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgei,ygei,zgei   cartesian gei coordinates
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq
c
      xgeo=  cgst*xgei + sgst*ygei
      ygeo= -sgst*xgei + cgst*ygei
      zgeo=  zgei
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gei_to_gse(xgei,ygei,zgei,xgse,ygse,zgse)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gei_to_gse: GEI -> GSE  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgei,ygei,zgei cartesian gei coordinates
c *   Output : xgse,ygse,zgse cartesian gse coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3
c
      xgse=   gs1*xgei +   gs2*ygei +   gs3*zgei
      ygse= gegs1*xgei + gegs2*ygei + gegs3*zgei
      zgse=   ge1*xgei +   ge2*ygei +   ge3*zgei
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gei_to_gsm(xgei,ygei,zgei,xgsm,ygsm,zgsm)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gei_to_gsm: GEI -> GSM  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgei,ygei,zgei   cartesian gei coordinates
c *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3
c
      xgsm=    gs1*xgei +    gs2*ygei +    gs3*zgei
      ygsm= yeigm1*xgei + yeigm2*ygei + yeigm3*zgei
      zgsm= zeigm1*xgei + zeigm2*ygei + zeigm3*zgei
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gei_to_gseq(xgei,ygei,zgei,xgsq,ygsq,zgsq)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gei_to_gsq: GEI -> GSEQ system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgei,ygei,zgei   cartesian gei coordinates
c *   Output : xgsq,ygsq,zgsq   cartesian gsq coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3
c
      xgsq=    gs1*xgei +    gs2*ygei +    gs3*zgei
      ygsq= yeigq1*xgei + yeigq2*ygei + yeigq3*zgei
      zgsq= zeigq1*xgei + zeigq2*ygei + zeigq3*zgei
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gei_to_mag(xgei,ygei,zgei,xmag,ymag,zmag)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gei_to_mag: GEI -> MAG  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgei,ygei,zgei   cartesian gei coordinates
c *   Output : xmag,ymag,zmag   cartesian mag coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3
c
      xmag= xeima1*xgei + xeima2*ygei + xeima3*zgei
      ymag= yeima1*xgei + yeima2*ygei + yeima3*zgei
      zmag=    gm1*xgei +    gm2*ygei +    gm3*zgei
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gei_to_sm(xgei,ygei,zgei,xsma,ysma,zsma)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gei_to_sma: GEI -> SM   system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgei,ygei,zgei   cartesian gei coordinates
c *   Output : xsma,ysma,zsma   cartesian sma coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3
c
      xsma= xeism1*xgei + xeism2*ygei + xeism3*zgei
      ysma= yeism1*xgei + yeism2*ygei + yeism3*zgei
      zsma=    gm1*xgei +    gm2*ygei +    gm3*zgei
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_dm(xgeo,ygeo,zgeo,rlat,rlong,xdme,ydme,zdme)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_dme: GEO -> DM   system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates
c              rlat,rlong       latitude and longitude of the point
c                               of observation (radians)
c *   Output : xdme,ydme,zdme   cartesian dme coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
c
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
c
      y1= real(dprod(gd2,r3) - dprod(gd3,r2))
      y2= real(dprod(gd3,r1) - dprod(gd1,r3))
      y3= real(dprod(gd1,r2) - dprod(gd2,r1))
c
      q=sqrt(y1**2 + y2**2 + y3**2)
c
      y1=y1/q
      y2=y2/q
      y3=y3/q
c
      x1= real(dprod(y2,gd3) - dprod(y3,gd2))
      x2= real(dprod(y3,gd1) - dprod(y1,gd3))
      x3= real(dprod(y1,gd2) - dprod(y2,gd1))
c
      xdme=  x1*xgeo +  x2*ygeo +  x3*zgeo
      ydme=  y1*xgeo +  y2*ygeo +  y3*zgeo
      zdme= gd1*xgeo + gd2*ygeo + gd3*zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_gei(xgeo,ygeo,zgeo,xgei,ygei,zgei)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_gei: GEO -> GEI  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgeo,ygeo,zgeo cartesian geo coordinates
c *   Output : xgei,ygei,zgei cartesian gei coordinates
c ----------------------------------------------------------------------
c
c
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq
c
      xgei= cgst*xgeo - sgst*ygeo
      ygei= sgst*xgeo + cgst*ygeo
      zgei= zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_gse(xgeo,ygeo,zgeo,xgse,ygse,zgse)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_gse: GEO -> GSE  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates
c *   Output : xgse,ygse,zgse   cartesian gse coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3
c
      xgse=   ps1*xgeo +   ps2*ygeo +   ps3*zgeo
      ygse= peps1*xgeo + peps2*ygeo + peps3*zgeo
      zgse=   pe1*xgeo +   pe2*ygeo +   pe3*zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_gsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_gsm: GEO -> GSM  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates
c *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3
c
      xgsm=    ps1*xgeo +    ps2*ygeo +    ps3*zgeo
      ygsm= yeogm1*xgeo + yeogm2*ygeo + yeogm3*zgeo
      zgsm= zeogm1*xgeo + zeogm2*ygeo + zeogm3*zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_gseq(xgeo,ygeo,zgeo,xgsq,ygsq,zgsq)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_gsq: GEO -> GSEQ system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates
c *   Output : xgsq,ygsq,zgsq   cartesian gsq coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3
c
      xgsq=    ps1*xgeo +    ps2*ygeo +    ps3*zgeo
      ygsq= yeogq1*xgeo + yeogq2*ygeo + yeogq3*zgeo
      zgsq= zeogq1*xgeo + zeogq2*ygeo + zeogq3*zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_mag(xgeo,ygeo,zgeo,xmag,ymag,zmag)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_mag: GEO -> MAG  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates
c *   Output : xmag,ymag,zmag   cartesian mag coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3
c
      xmag= xeoma1*xgeo + xeoma2*ygeo + xeoma3*zgeo
      ymag= yeoma1*xgeo + yeoma2*ygeo + yeoma3*zgeo
      zmag=    gd1*xgeo +    gd2*ygeo +    gd3*zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_sm(xgeo,ygeo,zgeo,xsma,ysma,zsma)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_sma: GEO -> SM   system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates
c *   Output : xsma,ysma,zsma   cartesian sma coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3
c
      xsma= xeosm1*xgeo + xeosm2*ygeo + xeosm3*zgeo
      ysma= yeosm1*xgeo + yeosm2*ygeo + yeosm3*zgeo
      zsma=    gd1*xgeo +    gd2*ygeo +    gd3*zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_geo_to_vdh(xgeo,ygeo,zgeo,rlat,rlong,xvdh,yvdh,zvdh)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_geo_to_vdh: GEO -> VDH  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: local system, non time dependent
c
c *   Input  : xgeo,ygeo,zgeo   cartesian geo coordinates
c              rlat,rlong       latitude and longitude of the point
c                               of observation (radians)
c *   Output : xvdh,yvdh,zvdh   cartesian vdh coordinates
c ----------------------------------------------------------------------
c
c
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
c
      v1=r1
      v2=r2
      v3=r3
c
      q12=sqrt(r1**2 + r2**2)
c
      d1= -r2/q12
      d2=  r1/q12
      d3=  0.
c
      h1= -r1*r3/q12
      h2= -r2*r3/q12
      h3=  q12
c
      xvdh= v1*xgeo + v2*ygeo + v3*zgeo
      yvdh= d1*xgeo + d2*ygeo + d3*zgeo
      zvdh= h1*xgeo + h2*ygeo + h3*zgeo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gse_to_gei(xgse,ygse,zgse,xgei,ygei,zgei)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gse_to_gei: GSE -> GEI  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgse,ygse,zgse    cartesian gse coordinates
c *   Output : xgei,ygei,zgei    cartesian gei coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3
c
      xgei= gs1*xgse + gegs1*ygse + ge1*zgse
      ygei= gs2*xgse + gegs2*ygse + ge2*zgse
      zgei= gs3*xgse + gegs3*ygse + ge3*zgse
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gse_to_geo(xgse,ygse,zgse,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gse_to_geo: GSE -> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgse,ygse,zgse   cartesian gse coordinates
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3
c
      xgeo= ps1*xgse + peps1*ygse + pe1*zgse
      ygeo= ps2*xgse + peps2*ygse + pe2*zgse
      zgeo= ps3*xgse + peps3*ygse + pe3*zgse
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gse_to_gsm(xgse,ygse,zgse,xgsm,ygsm,zgsm)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gse_to_gsm: GSE -> GSM  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgse,ygse,zgse   cartesian gse coordinates
c *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xgsm= xgse
      ygsm=  cdze*ygse + sdze*zgse
      zgsm= -sdze*ygse + cdze*zgse
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gse_to_gseq(xgse,ygse,zgse,xgsq,ygsq,zgsq)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gse_to_gsq: GSE -> GSEQ system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgse,ygse,zgse   cartesian gse coordinates
c *   Output : xgsq,ygsq,zgsq   cartesian gsq coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xgsq= xgse
      ygsq= ctetq*ygse - stetq*zgse
      zgsq= stetq*ygse + ctetq*zgse
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gse_to_mfa(xgse,ygse,zgse,bx,by,bz,xmfa,ymfa,zmfa)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gse_to_mfa: GSE -> MFA  system
c *   Author : P. Robert, LPP , 2016
c *   Comment: local system, non time dependent
c
c *   Input  : xgse,ygse,gsez cartesian gse coordinates
c              bx,  by,  bz   cartesian gse coordinates of DC mag field
c *   Output : xmfa,ymfa,zmfa cartesian mfa coordinates
c ----------------------------------------------------------------------
c
c
c *** tranform the vector from gse to mfa
c
      bperp= sqrt(bx*bx + by*by)
      b0= sqrt(bx*bx + by*by + bz*bz)
c
c **  first rotation
c
      sinphi=by/bperp
      cosphi=bx/bperp
c
      xp=  cosphi*xgse + sinphi*ygse
      yp= -sinphi*xgse + cosphi*ygse
      zp=  zgse
c
c **  second rotation
c
      sintet= bperp/b0
      costet= bz/b0
c
      xmag= costet*xp - sintet*zp
      ymag= yp
      zmag= sintet*xp + costet*zp
c
c **  third rotation
c
      sxm= costet*cosphi
      sym=-sinphi
c
      smperp= sqrt(sxm*sxm + sym*sym)
c
      sinpsi= sym/smperp
      cospsi= sxm/smperp
c
      xmfa=  cospsi*xmag + sinpsi*ymag
      ymfa= -sinpsi*xmag + cospsi*ymag
      zmfa=  zmag
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gse_to_sr2(xgse,ygse,zgse,rotx,roty,rotz,
     &                        xsr2,ysr2,zsr2)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gse_to_sr2: GSE -> SR2  system
c *   Author : P. Robert, CETP, 2001
c *   Comment: local system, non time dependent
c
c *   Input  : xgse,ygse,zgse cartesian gse coordinates
c              rotx,roty,rotz cartesian gse coordinates of rotation axis
c *   Output : xsr2,ysr2,zsr2 cartesian sr2 coordinates
c ----------------------------------------------------------------------
c
c
c *** set transform matrix with Spin axis terms
c
      rmod= sqrt(rotx**2 + roty**2 + rotz**2)
c
      rx=rotx/rmod
      ry=roty/rmod
      rz=rotz/rmod
c
      a= 1./sqrt(ry*ry + rz*rz)
c
      x1= (ry*ry + rz*rz)*a
      x2= -rx*ry*a
      x3= -rx*rz*a
c
      y1=  0.
      y2=  rz*a
      y3= -ry*a
c
      z1= rx
      z2= ry
      z3= rz
c
c *** tranform the input vector from gse to sr2
c
      xsr2= x1*xgse + x2*ygse + x3*zgse
      ysr2= y1*xgse + y2*ygse + y3*zgse
      zsr2= z1*xgse + z2*ygse + z3*zgse
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gse_to_tpn(xgse,ygse,zgse,xo,yo,zo,xs,xtpn,ytpn,ztpn)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gse_to_tpn: GSE -> TPN  system
c *   Author : P. Robert, LPP , 2016
c *   Comment: local system, non time dependent
c
c *   Input  : xgse,ygse,zgse cartesian gse vector
c              xo,yo,zo position of the S/C in gse
c              xs subsolar point, submit of the paraboloid
c              from Earth to Sun
c *   Output : xtpn,ytpn,ztpn cartesian tpn coordinates
c ----------------------------------------------------------------------
c
      real Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz
c
c     computation of the TPN axis in gse system
c
      call cp_tpn_param(xo,yo,zo,xs, Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz)
c
      xtpn= Tx*xgse + Ty*ygse +Tz*zgse
      ytpn= Px*xgse + Py*ygse +Pz*zgse
      ztpn= Nx*xgse + Ny*ygse +Nz*zgse
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gsm_to_gei(xgsm,ygsm,zgsm,xgei,ygei,zgei)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsm_to_gei: GSM -> GEI  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
c *   Output : xgei,ygei,zgei   cartesian gei coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3
c
      xgei= gs1*xgsm + yeigm1*ygsm + zeigm1*zgsm
      ygei= gs2*xgsm + yeigm2*ygsm + zeigm2*zgsm
      zgei= gs3*xgsm + yeigm3*ygsm + zeigm3*zgsm
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gsm_to_geo(xgsm,ygsm,zgsm,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsm_to_geo: GSM -> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3
c
      xgeo= ps1*xgsm + yeogm1*ygsm + zeogm1*zgsm
      ygeo= ps2*xgsm + yeogm2*ygsm + zeogm2*zgsm
      zgeo= ps3*xgsm + yeogm3*ygsm + zeogm3*zgsm
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gsm_to_gse(xgsm,ygsm,zgsm,xgse,ygse,zgse)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsm_to_gse: GSM -> GSE  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
c *   Output : xgse,ygse,zgse   cartesian gse coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xgse= xgsm
      ygse= cdze*ygsm - sdze*zgsm
      zgse= sdze*ygsm + cdze*zgsm
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gsm_to_gseq(xgsm,ygsm,zgsm,xgsq,ygsq,zgsq)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsm_to_gsq: GSM -> GSQ  system
c *   Author : P. Robert, CRPE, 2002
c
c *   Input  : xgsm,ygsm,zgsm cartesian gsm coordinates
c *   Output : xgsq,ygsq,zgsq cartesian gsq coordinates
c ----------------------------------------------------------------------
c
c
      call t_gsm_to_gse(xgsm,ygsm,zgsm,ax,ay,az)
      call t_gse_to_gseq(ax,ay,az,xgsq,ygsq,zgsq)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gsm_to_mag(xgsm,ygsm,zgsm,xmag,ymag,zmag)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsm_to_mag: GSM -> MAG  system
c *   Author : P. Robert, CRPE, 2002
c
c *   Input  : xgsm,ygsm,zgsm cartesian gsm coordinates
c *   Output : xmag,ymag,zmag cartesian mag coordinates
c ----------------------------------------------------------------------
c
c
      call t_gsm_to_sm(xgsm,ygsm,zgsm,ax,ay,az)
      call t_sm_to_mag(ax,ay,az,xmag,ymag,zmag)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gsm_to_sm(xgsm,ygsm,zgsm,xsma,ysma,zsma)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsm_to_sma: GSM -> SM   system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
c *   Output : xsma,ysma,zsma   cartesian sma coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xsma= cmu*xgsm - smu*zgsm
      ysma= ygsm
      zsma= smu*xgsm + cmu*zgsm
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gsm_to_tpn(xgsm,ygsm,zgsm,xo,yo,zo,xs,xtpn,ytpn,ztpn)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsm_to_tpn: GSM -> TPN  system
c *   Author : P. Robert, LPP , 2016
c *   Comment: local system, non time dependent
c
c *   Input  : xgsm,ygsm,zgsm cartesian gsm vector
c              xo,yo,zo position of the S/C in gsm
c              xs subsolar point, submit of the paraboloid
c              from Earth to Sun
c *   Output : xtpn,ytpn,ztpn cartesian tpn coordinates
c ----------------------------------------------------------------------
c
      real Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz
c
c     computation of the TPN axis in gsm system
c
      call cp_tpn_param(xo,yo,zo,xs, Tx,Ty,Tz, Px,Py,Pz, Nx,Ny,Nz)
c
      xtpn= Tx*xgsm + Ty*ygsm +Tz*zgsm
      ytpn= Px*xgsm + Py*ygsm +Pz*zgsm
      ztpn= Nx*xgsm + Ny*ygsm +Nz*zgsm
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gseq_to_gei(xgsq,ygsq,zgsq,xgei,ygei,zgei)
c
c ----------------------------------------------------------------------
c
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsq_to_gei: GSEQ-> GEI  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgsq,ygsq,zgsq   cartesian gsq coordinates
c *   Output : xgei,ygei,zgei   cartesian gei coordinates
c
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3
c
      xgei= gs1*xgsq + yeigq1*ygsq + zeigq1*zgsq
      ygei= gs2*xgsq + yeigq2*ygsq + zeigq2*zgsq
      zgei= gs3*xgsq + yeigq3*ygsq + zeigq3*zgsq
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gseq_to_geo(xgsq,ygsq,zgsq,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsq_to_geo: GSEQ-> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgsq,ygsq,zgsq   cartesian gsq coordinates
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3
c
      xgeo= ps1*xgsq + yeogq1*ygsq + zeogq1*zgsq
      ygeo= ps2*xgsq + yeogq2*ygsq + zeogq2*zgsq
      zgeo= ps3*xgsq + yeogq3*ygsq + zeogq3*zgsq
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gseq_to_gse(xgsq,ygsq,zgsq,xgse,ygse,zgse)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsq_to_gse: GSEQ-> GSE  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xgsq,ygsq,zgsq   cartesian gsq coordinates
c *   Output : xgse,ygse,zgse   cartesian gse coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xgse= xgsq
      ygse= ctetq*ygsq + stetq*zgsq
      zgse=-stetq*ygsq + ctetq*zgsq
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_gseq_to_gsm(xgsq,ygsq,zgsq,xgsm,ygsm,zgsm)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_gsq_to_gsm: GSQ -> GSM  system
c *   Author : P. Robert, CRPE, 2002
c
c *   Input  : xgsq,ygsq,zgsq cartesian gsq coordinates
c *   Output : xgsm,ygsm,zgsm cartesian gsm coordinates
c ----------------------------------------------------------------------
c
c
      call t_gseq_to_gse(xgsq,ygsq,zgsq,xx,yy,zz)
      call t_gse_to_gsm(xx,yy,zz,xgsm,ygsm,zgsm)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_mag_to_gei(xmag,ymag,zmag,xgei,ygei,zgei)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_mag_to_gei: MAG -> GEI  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xmag,ymag,zmag cartesian mag coordinates
c *   Output : xgei,ygei,zgei cartesian gei coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3
c
      xgei= xeima1*xmag + yeima1*ymag + gm1*zmag
      ygei= xeima2*xmag + yeima2*ymag + gm2*zmag
      zgei= xeima3*xmag + yeima3*ymag + gm3*zmag
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_mag_to_geo(xmag,ymag,zmag,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_mag_to_geo: MAG -> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xmag,ymag,zmag   cartesian mag coordinates
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3
c
      xgeo= xeoma1*xmag + yeoma1*ymag + gd1*zmag
      ygeo= xeoma2*xmag + yeoma2*ymag + gd2*zmag
      zgeo= xeoma3*xmag + yeoma3*ymag + gd3*zmag
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_mag_to_gsm(xmag,ymag,zmag,xgsm,ygsm,zgsm)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_mag_to_gsm: MAG -> GSM  system
c *   Author : P. Robert, CRPE, 2002
c
c *   Input  : xmag,ymag,zmag cartesian mag coordinates
c *   Output : xgsm,ygsm,zgsm cartesian gsm coordinates
c ----------------------------------------------------------------------
c
c
      call t_mag_to_sm(xmag,ymag,zmag,xx,yy,zz)
      call t_sm_to_gsm(xx,yy,zz,xgsm,ygsm,zgsm)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_mag_to_sm(xmag,ymag,zmag,xsma,ysma,zsma)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_mag_to_sma: MAG -> SM   system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xmag,ymag,zmag   cartesian mag coordinates
c *   Output : xsma,ysma,zsma   cartesian sma coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xsma= cphi*xmag - sphi*ymag
      ysma= sphi*xmag + cphi*ymag
      zsma= zmag
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sm_to_gei(xsma,ysma,zsma,xgei,ygei,zgei)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sma_to_gei: SM  -> GEI  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xsma,ysma,zsma   cartesian sma coordinates
c *   Output : xgei,ygei,zgei   cartesian gei coordinates
c ----------------------------------------------------------------------
c
c
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3
c
      xgei=xeism1*xsma + yeism1*ysma + gm1*zsma
      ygei=xeism2*xsma + yeism2*ysma + gm2*zsma
      zgei=xeism3*xsma + yeism3*ysma + gm3*zsma
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sm_to_geo(xsma,ysma,zsma,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sma_to_geo: SM  -> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xsma,ysma,zsma   cartesian sma coordinates
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3
c
      xgeo= xeosm1*xsma + yeosm1*ysma + gd1*zsma
      ygeo= xeosm2*xsma + yeosm2*ysma + gd2*zsma
      zgeo= xeosm3*xsma + yeosm3*ysma + gd3*zsma
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sm_to_gsm(xsma,ysma,zsma,xgsm,ygsm,zgsm)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sma_to_gsm: SM  -> GSM  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xsma,ysma,zsma   cartesian sma coordinates
c *   Output : xgsm,ygsm,zgsm   cartesian gsm coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xgsm=  cmu*xsma + smu*zsma
      ygsm=  ysma
      zgsm= -smu*xsma + cmu*zsma
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sm_to_mag(xsma,ysma,zsma,xmag,ymag,zmag)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sma_to_mag: SM  -> MAG  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: terms of transformation matrix are given in common
c
c *   Input  : xsma,ysma,zsma   cartesian sma coordinates
c *   Output : xmag,ymag,zmag   cartesian mag coordinates
c ----------------------------------------------------------------------
c
c
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
c
      xmag=  cphi*xsma + sphi*ysma
      ymag= -sphi*xsma + cphi*ysma
      zmag=  zsma
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sph_to_car(r,teta,phi,x,y,z)
c
c ----------------------------------------------------------------------
c
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sph_to_car: SPH -> CAR  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: none
c
c *   Input  :   r,teta,phi  spherical coordinates (angles in radians)
c *   Output :   x,y,z       cartesian coordinates
c ----------------------------------------------------------------------
c
c
      sq=r*sin(teta)
      x=sq*cos(phi)
      y=sq*sin(phi)
      z=r*cos(teta)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sr2_to_gse(xsr2,ysr2,zsr2,rotx,roty,rotz,
     &                        xgse,ygse,zgse)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sr2_to_gse: SR2 -> GSE  system
c *   Author : P. Robert, CETP, 2001
c *   Comment: local system, non time dependent
c
c *   Input  : xsr2,ysr2,sr2z cartesian sr2 coordinates
c              rotx,roty,rotz cartesian gse coordinates of rotation axis
c *   Output : xgse,ygse,zgse cartesian gse coordinates
c ----------------------------------------------------------------------
c
c *** set transform matrix with Spin axis terms
c
      rmod= sqrt(rotx**2 + roty**2 + rotz**2)
c
      rx=rotx/rmod
      ry=roty/rmod
      rz=rotz/rmod
c
      a= 1./sqrt(ry*ry + rz*rz)
c
      x1= (ry*ry + rz*rz)*a
      x2= -rx*ry*a
      x3= -rx*rz*a
c
      y1=  0.
      y2=  rz*a
      y3= -ry*a
c
      z1= rx
      z2= ry
      z3= rz
c
c *** tranform the input vector from sr2 to gse
c
      xgse= x1*xsr2 + y1*ysr2 + z1*zsr2
      ygse= x2*xsr2 + y2*ysr2 + z2*zsr2
      zgse= x3*xsr2 + y3*ysr2 + z3*zsr2
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sr2_to_mfa(xsr2,ysr2,zsr2,bx,by,bz,rox,roy,roz,
     &                        xm,ym,zm)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sr2_to_mfa: SR2 -> MFA  system
c *   Author : P. Robert, CETP, 2001
c *   Comment: local system, non time dependent
c
c *   Input  : xsr2,ysr2,sr2z cartesian sr2 coordinates
c              bx,  by,  bz   cartesian sr2 coordinates of DC mag field
c              rox, roy, roz  cartesian gse coordinates of rotation axis
c *   Output : xm,  ym,  zm   cartesian mfa coordinates
c ----------------------------------------------------------------------
c
c
c *** normalise R supposed undependant of time
c
      rmod= sqrt(rox**2 + roy**2 + roz**2)
c
      rx=rox/rmod
      ry=roy/rmod
      rz=roz/rmod
c
c *** tranform the vector from sr2 to mfa
c
      bperp= sqrt(bx*bx + by*by)
      b0= sqrt(bx*bx + by*by + bz*bz)
c
c **  first rotation
c
      sinphi=by/bperp
      cosphi=bx/bperp
c
      xp=  cosphi*xsr2 + sinphi*ysr2
      yp= -sinphi*xsr2 + cosphi*ysr2
      zp=  zsr2
c
c **  second rotation
c
      sintet= bperp/b0
      costet= bz/b0
c
      xmag= costet*xp - sintet*zp
      ymag= yp
      zmag= sintet*xp + costet*zp
c
c **  third rotation
c
      b= sqrt(ry*ry + rz*rz)
c
      sxm=  b*costet*cosphi - rx*sintet
      sym= -b*sinphi
c
      smperp= sqrt(sxm*sxm + sym*sym)
c
      sinpsi= sym/smperp
      cospsi= sxm/smperp
c
      xm=  cospsi*xmag + sinpsi*ymag
      ym= -sinpsi*xmag + cospsi*ymag
      zm=  zmag
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sr2_to_sr(xsr2,ysr2,spifre,spipha,deltaT,xsre,ysre)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sr2_to_sre: SR2 -> SRef system
c *   Author : P. Robert, CRPE, 2001
c *   Comment: Z component is unchanged (spin axis)
c
c *   Input  : xsr2, ysr2 cartesian sr2 coordinates
c              spifre     spin frequency in Hz
c              spipha     spin phase in radians, growing with time
c                         spipha= positive angle between the xsr axis
c                         and the component of the direction of the Sun
c                         in the xsr-ysr plane.
c              deltaT     (T -To) time (sec.), between the current time
c                         and the time where is measured the spin phase.
c *   Output : xsre,ysre  cartesian sr coordinates
c ----------------------------------------------------------------------
c
c
      pi2= 2.*3.1415927
c
      depift= mod(pi2*spifre*deltaT,pi2)
      phicr= spipha -depift
      sinphi=sin(phicr)
      cosphi=cos(phicr)
c
      xsre=  cosphi*xsr2 -sinphi*ysr2
      ysre=  sinphi*xsr2 +cosphi*ysr2
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_sr_to_sr2(xsre,ysre,spifre,spipha,deltaT,xsr2,ysr2)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_sre_to_sr2: SRef-> SR2 system
c *   Author : P. Robert, CRPE, 2001
c *   Comment: Z component is unchanged (spin axis)
c
c *   Input  : xsre, ysre cartesian sr coordinates
c              spifre     spin frequency in Hz
c              spipha     spin phase in radians, growing with time
c                         spipha= positive angle between the xsr axis
c                         and the component of the direction of the Sun
c                         in the xsr-ysr plane.
c              deltaT     (T -To) time (sec.), between the current time
c                         and the time where is measured the spin phase.
c *   Output : xsr2,ysr2  cartesian sr coordinates
c ----------------------------------------------------------------------
c
c
      pi2= 2.*3.1415927
c
      depift= mod(pi2*spifre*deltaT,pi2)
      phicr= spipha -depift
      sinphi=sin(phicr)
      cosphi=cos(phicr)
c
      xsr2=  cosphi*xsre +sinphi*ysre
      ysr2= -sinphi*xsre +cosphi*ysre
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_vdh_to_geo(xvdh,yvdh,zvdh,rlat,rlong,xgeo,ygeo,zgeo)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_vdh_to_geo: VDH -> GEO  system
c *   Author : P. Robert, CRPE, 1992
c *   Comment: local system, non time dependent
c
c *   Input  : xvdh,yvdh,zvdh   cartesian vdh coordinates
c              rlat,rlong       latitude and longitude of the point
c                               of observation (radians)
c *   Output : xgeo,ygeo,zgeo   cartesian geo coordinates
c ----------------------------------------------------------------------
c
c
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
c
      v1=r1
      v2=r2
      v3=r3
c
      q12=sqrt(r1**2 + r2**2)
c
      d1= -r2/q12
      d2=  r1/q12
      d3=  0.
c
      h1= -r1*r3/q12
      h2= -r2*r3/q12
      h3=  q12
c
      xgeo= v1*xvdh + d1*yvdh + h1*zvdh
      ygeo= v2*xvdh + d2*yvdh + h2*zvdh
      zgeo= v3*xvdh + d3*yvdh + h3*zvdh
c
      return
      end
c      
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_xyz_to_vdh(x,y,z,a1,a2,a3,v,d,h)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_xyz_to_vdh: xyz spinning -> VDH
c *   Author : P. Robert, SDev, 2020
c *   Comment: use Euler angles in degrees; VDH is the fixed system.
c
c *   Input : x,y,z cartesian xyz coordinates
c *   Output: v,d,h cartesian VDH coordinates
c ----------------------------------------------------------------------

      pisd=3.1415926/180.
c
      a1r=a1*pisd
      a2r=a2*pisd
      a3r=a3*pisd
c
      c1=cos(a1r)
      c2=cos(a2r)
      c3=cos(a3r)

      s1=sin(a1r)
      s2=sin(a2r)
      s3=sin(a3r)
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine t_vdh_to_xyz(v,d,h,a1,a2,a3,x,y,z)
c
c ----------------------------------------------------------------------
c *   Class  : transform modules of Rocotlib Software
c *   Object : transforms_vdh_to_xyz: VDH -> xyz spinning 
c *   Author : P. Robert, SDev, 2020
c *   Comment: use Euler angles in degrees; VDH is the fixed system.
c
c *   Input  : v,d,h cartesian VDH coordinates
c *   Output : x,y,z cartesian xyz coordinates
c ----------------------------------------------------------------------

      pisd=3.1415926/180.
c
      a1r=a1*pisd
      a2r=a2*pisd
      a3r=a3*pisd
c
      c1=cos(a1r)
      c2=cos(a2r)
      c3=cos(a3r)

      s1=sin(a1r)
      s2=sin(a2r)
      s3=sin(a3r)
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
c     Supplements for matrix operations.
c     Allows the transformation of an input 3D signal into the
c     Minimum variance coordinate (MVA Analysis).
c
c     All codes in this part are extracted of Roproc Software, V 4.5
c     and rewrited properly in f77, with Rocotlib V2.2 conventions.
c     P. Robert, LPP , April 2016.
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_cp_varmin(ifc,Vx,Vy,Vz,n,irep,covar,lambda,eigvec)
c
      real Vx(n),Vy(n),Vz(n)
      real covar(3,3),lambda(3),eigvec(3,3)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : compute variance minimum coordinate of a signal Vx,Vy,Vz
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c     Comment: check orthogonality of the MVA matrix
c
c *   Input  : ifc (unit for writing results)
c              Vx(n),Vy(n),Vz(n) input signal
c              n number of point of the signal
c              irep=1 : passe input vector in MVA coordinates
c
c *   Output : covar covariance matrix
c              lamda(3), eigvev(3,3) eigen values and eigen vectors
c              results are also writted on unit ifc
c     ---------------------------------------------------------------+--
c
      write(ifc,100)
      write(ifc,100) '------------------------------------------------'
      write(ifc,100) 'Compute Minimum Variance Analaysis coord. system'
      write(ifc,100) '------------------------------------------------'
c
c *** compute covariance matrix
c
      call mat_cp_covariance(Vx,Vy,Vz,n,covar)
      call mat_write(ifc,'covariance matrix of input signal:',covar)
c
c *** diagonalisation of covariance matrix
c
      call mat_diagonalise(covar,lambda,eigvec)
      call mat_write(ifc,'covariance matrix diagonalized:',covar)
      call mat_write(ifc,'eigen vectors matrix:',eigvec)
c
c *** check orthogonality and direct sens of the 3 vectors
c
      call mat_check_ortho(ifc,eigvec)
c
c *** print eigen vectors and eigen values
c
      call mat_write_eigen_vec(ifc,lambda,eigvec)
c
c *** transform input signal in MVA coordinates, ie. eigen vectors
c
      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
c
  100 format(1x,a)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_cp_covariance(Vx,Vy,Vz,n,covar)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : compute covariance matrix for a vector series V(n)
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : Vx(n),Vy(n),Vz(n) vector series
c *   Output : covar(3,3)
c     ---------------------------------------------------------------+--
c
      dimension Vx(n),Vy(n),Vz(n)
      dimension covar(3,3)
      double precision sx, sy, sz, sx2, sy2, sz2, sxy, syz, szx
c
c *   set summation to zero
c
      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

c
c *   components summation and square components summation
c
      do k=1,n
c
        sx= sx +dble(Vx(k))
        sy= sy +dble(Vy(k))
        sz= sz +dble(Vz(k))
c
        sx2= sx2 +dble(Vx(k)**2)
        sy2= sy2 +dble(Vy(k)**2)
        sz2= sz2 +dble(Vz(k)**2)
c
        sxy= sxy +dble(Vx(k)*Vy(k))
        syz= syz +dble(Vy(k)*Vz(k))
        szx= szx +dble(Vz(k)*Vx(k))
c
      enddo
c
c *   average values
c
      sx = sx/dble(n)
      sy = sy/dble(n)
      sz = sz/dble(n)
c
      sx2= sx2/dble(n)
      sy2= sy2/dble(n)
      sz2= sz2/dble(n)
c
      sxy= sxy/dble(n)
      syz= syz/dble(n)
      szx= szx/dble(n)
c
c *   variance set to diagonals terms
c
      covar(1,1)= sngl(sx2-sx**2)
      covar(2,2)= sngl(sy2-sy**2)
      covar(3,3)= sngl(sz2-sz**2)
c
c *   covariances set to semi \diagonal terms
c
      covar(2,1)= sngl(sxy-sx*sy)
      covar(3,2)= sngl(syz-sy*sz)
      covar(3,1)= sngl(szx-sz*sx)
c
c *   second half triangle
c
      covar(1,2)= covar(2,1)
      covar(1,3)= covar(3,1)
      covar(2,3)= covar(3,2)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_diagonalise(mat,lambda,eigvec)
c
      real mat(3,3),lambda(3),eigvec(3,3)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : diagonalise the given matrix mat(3,3)
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : real mat(3,3)
c *   Output : lambda(3)   eigen values
c              eigvec(3,3) eigen vectors normalised to 1.
c     ---------------------------------------------------------------+--
c
c     compute eigen vectors and eigen values of real mat(3,3)
c
      call mat_cp_eigen_vec(mat,lambda,eigvec)
c
c     normalise to 1 eigen vectors

      call mat_normalize_vec(eigvec)

c *   set positive direction for \diag
c
      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
c
c *   sort eigen vector with eigen value (z=min)
c
      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
c
c *   determinant computation:
c
      call mat_cp_determin(eigvec,det)
c
      if (det.lt.0.) then
         eigvec(1,2)= -eigvec(1,2)
         eigvec(2,2)= -eigvec(2,2)
         eigvec(3,2)= -eigvec(3,2)
      endif
c
c
c *   modif P.R. 2004: the vector corresponding to lambda min has z >0
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_check_ortho(ifc,mat)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : check orthogonality of matrix components
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : ifc (unit for writing results), real mat(3,3)
c *   Output : writing result on unit ifc
c     ---------------------------------------------------------------+--
c
      real mat(3,3)
c
      v1x=mat(1,1)
      v1y=mat(2,1)
      v1z=mat(3,1)
c
      v2x=mat(1,2)
      v2y=mat(2,2)
      v2z=mat(3,2)
c
      v3x=mat(1,3)
      v3y=mat(2,3)
      v3z=mat(3,3)
c
      write(ifc,100)
      write(ifc,100) 'Check orthogonality of given matrix:'
      write(ifc,100) '----------------------------------- '
      write(ifc,100)
c
      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
c
      write(ifc,100)
      write(ifc,100) 'V1.V2 =',v1v2
      write(ifc,100) 'V2.V3 =',v2v3
      write(ifc,100) 'V3.V1 =',v3v1
c
      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):'
c
c *   v1 X v2
c
      v1cv2x= v1y*v2z - v1z*v2y
      v1cv2y= v1z*v2x - v1x*v2z
      v1cv2z= v1x*v2y - v1y*v2x
c
c *   v2 X v3
c
      v2cv3x= v2y*v3z - v2z*v3y
      v2cv3y= v2z*v3x - v2x*v3z
      v2cv3z= v2x*v3y - v2y*v3x
c
c *   v3 X v1
c
      v3cv1x= v3y*v1z - v3z*v1y
      v3cv1y= v3z*v1x - v3x*v1z
      v3cv1z= v3x*v1y - v3y*v1x
c
      write(ifc,100)
      write(ifc,200) 'V1XV2 =',v1cv2x,v1cv2y,v1cv2z
      write(ifc,200) '   V3 =',v3x,v3y,v3z
c
      write(ifc,100)
      write(ifc,200) 'V2XV3 =',v2cv3x,v2cv3y,v2cv3z
      write(ifc,200) '   V1 =',v1x,v1y,v1z
c
      write(ifc,100)
      write(ifc,200) 'V3XV1 =',v3cv1x,v3cv1y,v3cv1z
      write(ifc,200) '   V2 =',v2x,v2y,v2z
c
  100 format(3(1x,a,1Pe14.5))
  200 format(1x,a,3(1Pe16.5))

c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_cp_determin(mat,det)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : compute determinant of the given matrix
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : real mat(3,3)
c *   Output : det
c     ---------------------------------------------------------------+--
c
      real mat(3,3)
c
      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)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_cp_eigen_vec(mat,lambda,eigvec)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : compute eigen vectors and eigen values of real mat(3,3)
c *   Author : unknown, CETP, 2001, rev. PR 2016
c *   Comment: mat(3,3) must be real and symmetrical;
c              Method used is Householder.
c
c *   Input  : real mat(3,3)
c *   Output : lambda(3)   eigen values
c              eigvec(3,3) eigen vectors
c     ---------------------------------------------------------------+--
c
      real mat(3,3),lambda(3),eigvec(3,3),rr(3,3),e(3)
c
c *   set eigen vectors to unity matrix
c
      do i=1,3
         do j=1,3
            eigvec(i,j)=0.0
         enddo
         eigvec(i,i)=1.0
      enddo
c
c     set eigen values
c
      lambda(1)=mat(1,3)
      lambda(2)=mat(2,3)
c
      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
c
      h=(lambda(1)**2+lambda(2)**2)/2.0
c
      epsi=1.e-30
      
      if (abs(h).gt.epsi) then
c
c        determination of the Householder matrix (results in  eigvec).
c
         do i=1,3
            do j=1,3
               eigvec(i,j)=eigvec(i,j)-(lambda(i)*lambda(j))/h
            enddo
         enddo
c
c        tridiagonalisation of mat matrix
c
         call mat_product(eigvec,mat,rr)
         call mat_product(rr,eigvec,mat)
c
      endif

c     set eigen values to diagonal terms
c
      do i=1,3
         lambda(i)=mat(i,i)
      enddo
c
      do i=1,2
         e(i)=mat(i+1,i)
      enddo
      e(3)=0.0
c
c *   search for eigen values (max iteration = 30)
c
      do l=1,3
         iter=0
  1      do m=l,2
            dd=abs(lambda(m))+abs(lambda(m+1))
c 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
c
            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
c
               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
c
            enddo

            lambda(l)=lambda(l)-p
            e(l)=g
            e(m)=0.0
            goto 1
c
         endif
      enddo
c
c *   modif PR  sep 11 2001 to load  mat with the result of the diag.
c
      do i=1,3
         mat(i,i)=lambda(i)
      enddo
c
      do i=1,2
         mat(i+1,i)=e(i)
         mat(i,i+1)=e(i)
      enddo
c
      mat(1,3)=e(3)
      mat(3,1)=e(3)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_cp_pythag_func(a,b,fpyth)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : Pythagore function of two real (used by mat_cp_eigen_vec)
c *   Author : unknown, CETP, 2001, rev. PR 2016
c
c *   Input  : a,b
c *   Output : fpyth
c     ---------------------------------------------------------------+--
c
      absa=abs(a)
      absb=abs(b)
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_normalize_vec(mat)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : normalize to 1. the vectors of the input matrix
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : real mat(3,3)
c *   Output : real mat(3,3) with vectors normalized to 1.
c     ---------------------------------------------------------------+--
c
      real mat(3,3),matmod
c
      epsi=1.e-37
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_product(mat1,mat2,mat3)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : matrix product of two given matrix of dim. 3
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : n, mat1(3,3), mat2(3,3)
c *   Output : mat3(3,3)
c     ---------------------------------------------------------------+--
c
      real mat1(3,3),mat2(3,3),mat3(3,3)
c
      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
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_somme(mat1,mat2,mat3)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : matrix somme of two given matrix of dim. 3
c *   Author : P. Robert, LPP , 2016
c
c *   Input  : n, mat1(3,3), mat2(3,3)
c *   Output : mat3(3,3)
c     ---------------------------------------------------------------+--
c
      real mat1(3,3),mat2(3,3),mat3(3,3)
c
      do i=1,3
         do j=1,3
            mat3(i,j)=mat1(i,j)+mat2(i,j)
         enddo
      enddo

      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_diff(mat1,mat2,mat3)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : matrix difference of two given matrix 
c *   Author : P. Robert, LPP , 2016
c
c *   Input  : n, mat1(3,3), mat2(3,3)
c *   Output : mat3(3,3)
c     ---------------------------------------------------------------+--
c
      real mat1(3,3),mat2(3,3),mat3(3,3)
c
      do i=1,3
         do j=1,3
            mat3(i,j)=mat1(i,j)-mat2(i,j)
         enddo
      enddo

      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_transpose(mat)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : transpose input matrix
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : real mat(3,3)
c *   Output : real mat(3,3) transposed
c     ---------------------------------------------------------------+--
c
      real mat(3,3)
c
      r12= mat(1,2)
      r13= mat(1,3)
      r23= mat(2,3)
c
      mat(1,2)=mat(2,1)
      mat(1,3)=mat(3,1)
      mat(2,3)=mat(3,2)
c
      mat(2,1)=r12
      mat(3,1)=r13
      mat(3,2)=r23
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_change_coord(mat,Vx,Vy,Vz,n)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : change coordinate of a vector serie with a given matrix
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c *   Comment: compute V(n)= mat*V(n)
c
c *   Input  : ifc (unit for writing results), real mat(3,3)
c *   Output : writing result on unit ifc
c     ---------------------------------------------------------------+--
c
      real Vx(n),Vy(n),Vz(n)
      real mat(3,3)
c
      do i=1,n
c
         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)
c
         Vx(i)=Vxr
         Vy(i)=Vyr
         Vz(i)=Vzr
c
      enddo
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_write(ifc,com,mat)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : print on ifc unit mat(3,3) with a comment
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : ifc (unit for writing results), com, real mat(3,3)
c *   Output : writing result on unit ifc
c     ---------------------------------------------------------------+--
c
      real mat(3,3)
      character*(*) com
c
      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)
c
  100 format(1x,a)
  200 format(11x,3(1Pe15.5))
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
      subroutine mat_write_eigen_vec(ifc,lambda,mat)
c
c     ---------------------------------------------------------------+--
c *   Class  : matrix operation of Rocotlib Software
c *   Object : print on ifc unit eigen values & vectors of mat(3,3)
c *   Author : P. Robert, CETP, 2001, rev. PR 2016
c
c *   Input  : ifc (unit for writing results), lambda(3), real mat(3,3)
c *   Output : writing result on unit ifc
c     ---------------------------------------------------------------+--
c
      real lambda(3),mat(3,3),r(3),theta(3),phi(3)
c
      dpi=180./3.14159
c
      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)
c
  100 format(1x,a,3f15.4)
  110 format(1x,a,3a15  )
  120 format(1x,a,3(1Pe15.4))
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
c
c BEGIN V2.2 compatibility 
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c     P. Robert, ScientiDev, Janvier 2019
c     subroutine for compatibility with previous V2.2 versions
c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cdatdoy(idoy,iyear,imonth,iday)
      call cv_doty_to_date(idoy,iyear,imonth,iday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cdatj00(jd00,iyear,imonth,iday)
      call cv_jul2000_to_date(jd00,iyear,imonth,iday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cdatj50(jd50,iyear,imonth,iday)
      call cv_jul1950_to_date(jd50,iyear,imonth,iday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cdatwee(iweek,iyear,imonth,iday)
      call cv_weekn_to_date(iweek,iyear,imonth,iday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cdipdir(iyear,idoy,d1,d2,d3)
      call cp_geo_dipole_dir(iyear,idoy,d1,d2,d3)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cdoweek(iyear,imonth,iday,idow)
      call cv_date_to_dotw(iyear,imonth,iday,idow)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cdoyear(iyear,imonth,iday,idoy)
      call cv_date_to_doty(iyear,imonth,iday,idoy)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cfrdayn(iday,cday,nbcha)
      character*(*) cday
      call cp_fr_day_name(iday,cday,nbcha)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx


      subroutine cfrmonn(imonth,cmonth,nchar)
      character*(*) cmonth
      call cp_fr_month_name(imonth,cmonth,nchar)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine chouday(ih,im,is,houday)
      call cv_hms_to_dech(ih,im,is,houday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cjd1950(iyear,imonth,iday,jd50)
      call cv_date_to_jul1950(iyear,imonth,iday,jd50)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cjd2000(iyear,imonth,iday,jd00)
      call cv_date_to_jul2000(iyear,imonth,iday,jd00)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cmilday(ih,im,is,ims,milday)
      call cv_dhms_to_msotd(ih,im,is,ims,milday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cnbdmon(iyear,imonth,nbday)
      call cp_nbday_in_month(iyear,imonth,nbday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine coleapy(iyear,ileap)
      call cp_leap_year(iyear,ileap)
      return
      end

c     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

c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine ctimhou(houday,ih,im,is)
      call cv_dech_to_hms(houday,ih,im,is)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine ctimmil(milday,ih,im,is,ims)
      call cv_msotd_to_hmsms(milday,ih,im,is,ims)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine ctimpa2(jd1950,houday)
      call cp_time_param2(jd1950,houday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine ctimpa3(jd2000,houday)
      call cp_time_param3(jd2000,houday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine ctimpar(iyear,imonth,iday,ih,im,is)
      call cp_time_param(iyear,imonth,iday,ih,im,is)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cusdayn(iday,cday,nbcha)
      character*(*) cday
      call cp_en_day_name(iday,cday,nbcha)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cusmonn(imonth,cmonth,nchar)
      character*(*) cmonth
      call cp_en_month_name(imonth,cmonth,nchar)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine cweedoy(iyear,imonth,iday,iweek)
      call cv_date_to_weekn(iyear,imonth,iday,iweek)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine gdipdir(dxgei,dygei,dzgei,dxgeo,dygeo,dzgeo)
      call g_gei_geo_dipole_dir(dxgei,dygei,dzgei,dxgeo,dygeo,dzgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine gdiptan(diptan)
      call g_gsm_dipole_tilt_angle(diptan)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine gecldir(exgei,eygei,ezgei,exgeo,eygeo,ezgeo)
      call g_gei_geo_ecliptic_dir(exgei,eygei,ezgei,exgeo,eygeo,ezgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine gsrodir(rxgei,rygei,rzgei,rxgeo,rygeo,rzgeo)
      call g_gei_geo_sun_rot(rxgei,rygei,rzgei,rxgeo,rygeo,rzgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine gsunpar(gmst,slon,sras,sdec,obli)
      call g_gei_sun_param(gmst,slon,sras,sdec,obli)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine gvernum(vernum,verdat)
      character*14 verdat
      call g_rocot_version_number(vernum,verdat)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine plibinf
      call print_rocot_info
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine recoor(x,y,z,cs)
      character*1 cs
      call r_coordinate_values(x,y,z,cs)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine recsys(csys)
      character*(*) csys
      call r_coordinate_system(csys)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine redate(iyear,imonth,iday)
      call r_date(iyear,imonth,iday)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine retime(ih,im,is)
      call r_time(ih,im,is)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tcarsph(x,y,z,r,teta,phi)
      call t_car_to_sph(x,y,z,r,teta,phi)
      return
      end

c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeigeo(xgei,ygei,zgei,xgeo,ygeo,zgeo)
      call t_gei_to_geo(xgei,ygei,zgei,xgeo,ygeo,zgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeigse(xgei,ygei,zgei,xgse,ygse,zgse)
      call t_gei_to_gse(xgei,ygei,zgei,xgse,ygse,zgse)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeigsm(xgei,ygei,zgei,xgsm,ygsm,zgsm)
      call t_gei_to_gsm(xgei,ygei,zgei,xgsm,ygsm,zgsm)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeigsq(xgei,ygei,zgei,xgsq,ygsq,zgsq)
      call t_gei_to_gseq(xgei,ygei,zgei,xgsq,ygsq,zgsq)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeimag(xgei,ygei,zgei,xmag,ymag,zmag)
      call t_gei_to_mag(xgei,ygei,zgei,xmag,ymag,zmag)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeisma(xgei,ygei,zgei,xsma,ysma,zsma)
      call t_gei_to_sm(xgei,ygei,zgei,xsma,ysma,zsma)
      return
      end

c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeogei(xgeo,ygeo,zgeo,xgei,ygei,zgei)
      call t_geo_to_gei(xgeo,ygeo,zgeo,xgei,ygei,zgei)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeogse(xgeo,ygeo,zgeo,xgse,ygse,zgse)
      call t_geo_to_gse(xgeo,ygeo,zgeo,xgse,ygse,zgse)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeogsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm)
      call t_geo_to_gsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeogsq(xgeo,ygeo,zgeo,xgsq,ygsq,zgsq)
      call t_geo_to_gseq(xgeo,ygeo,zgeo,xgsq,ygsq,zgsq)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeomag(xgeo,ygeo,zgeo,xmag,ymag,zmag)
      call t_geo_to_mag(xgeo,ygeo,zgeo,xmag,ymag,zmag)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgeosma(xgeo,ygeo,zgeo,xsma,ysma,zsma)
      call t_geo_to_sm(xgeo,ygeo,zgeo,xsma,ysma,zsma)
      return
      end

c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsegei(xgse,ygse,zgse,xgei,ygei,zgei)
      call t_gse_to_gei(xgse,ygse,zgse,xgei,ygei,zgei)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsegeo(xgse,ygse,zgse,xgeo,ygeo,zgeo)
      call t_gse_to_geo(xgse,ygse,zgse,xgeo,ygeo,zgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsegsm(xgse,ygse,zgse,xgsm,ygsm,zgsm)
      call t_gse_to_gsm(xgse,ygse,zgse,xgsm,ygsm,zgsm)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsegsq(xgse,ygse,zgse,xgsq,ygsq,zgsq)
      call t_gse_to_gseq(xgse,ygse,zgse,xgsq,ygsq,zgsq)
      return
      end

c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsmgei(xgsm,ygsm,zgsm,xgei,ygei,zgei)
      call t_gsm_to_gei(xgsm,ygsm,zgsm,xgei,ygei,zgei)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsmgeo(xgsm,ygsm,zgsm,xgeo,ygeo,zgeo)
      call t_gsm_to_geo(xgsm,ygsm,zgsm,xgeo,ygeo,zgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsmgse(xgsm,ygsm,zgsm,xgse,ygse,zgse)
      call t_gsm_to_gse(xgsm,ygsm,zgsm,xgse,ygse,zgse)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsmgsq(xgsm,ygsm,zgsm,xgsq,ygsq,zgsq)
      call t_gsm_to_gseq(xgsm,ygsm,zgsm,xgsq,ygsq,zgsq)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsmmag(xgsm,ygsm,zgsm,xmag,ymag,zmag)
      call t_gsm_to_mag(xgsm,ygsm,zgsm,xmag,ymag,zmag)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsmsma(xgsm,ygsm,zgsm,xsma,ysma,zsma)
      call t_gsm_to_sm(xgsm,ygsm,zgsm,xsma,ysma,zsma)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsqgei(xgsq,ygsq,zgsq,xgei,ygei,zgei)
      call t_gseq_to_gei(xgsq,ygsq,zgsq,xgei,ygei,zgei)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsqgeo(xgsq,ygsq,zgsq,xgeo,ygeo,zgeo)
      call t_gseq_to_geo(xgsq,ygsq,zgsq,xgeo,ygeo,zgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsqgse(xgsq,ygsq,zgsq,xgse,ygse,zgse)
      call t_gseq_to_gse(xgsq,ygsq,zgsq,xgse,ygse,zgse)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tgsqgsm(xgsq,ygsq,zgsq,xgsm,ygsm,zgsm)
      call t_gseq_to_gsm(xgsq,ygsq,zgsq,xgsm,ygsm,zgsm)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tmaggei(xmag,ymag,zmag,xgei,ygei,zgei)
      call t_mag_to_gei(xmag,ymag,zmag,xgei,ygei,zgei)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tmaggeo(xmag,ymag,zmag,xgeo,ygeo,zgeo)
      call t_mag_to_geo(xmag,ymag,zmag,xgeo,ygeo,zgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tmaggsm(xmag,ymag,zmag,xgsm,ygsm,zgsm)
      call t_mag_to_gsm(xmag,ymag,zmag,xgsm,ygsm,zgsm)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tmagsma(xmag,ymag,zmag,xsma,ysma,zsma)
      call t_mag_to_sm(xmag,ymag,zmag,xsma,ysma,zsma)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tsmagei(xsma,ysma,zsma,xgei,ygei,zgei)
      call t_sm_to_gei(xsma,ysma,zsma,xgei,ygei,zgei)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tsmageo(xsma,ysma,zsma,xgeo,ygeo,zgeo)
      call t_sm_to_geo(xsma,ysma,zsma,xgeo,ygeo,zgeo)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tsmagsm(xsma,ysma,zsma,xgsm,ygsm,zgsm)
      call t_sm_to_gsm(xsma,ysma,zsma,xgsm,ygsm,zgsm)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tsmamag(xsma,ysma,zsma,xmag,ymag,zmag)
      call t_sm_to_mag(xsma,ysma,zsma,xmag,ymag,zmag)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tsphcar(r,teta,phi,x,y,z)
      call t_sph_to_car(r,teta,phi,x,y,z)
      return
      end

c     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

c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tsr2sre(xsr2,ysr2,spifre,spipha,deltaT,xsre,ysre)
      call t_sr2_to_sr(xsr2,ysr2,spifre,spipha,deltaT,xsre,ysre)
      return
      end

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      subroutine tsresr2(xsre,ysre,spifre,spipha,deltaT,xsr2,ysr2)
      call t_sr_to_sr2(xsre,ysre,spifre,spipha,deltaT,xsr2,ysr2)
      return
      end

c     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

c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
c
c     modulo for f77
c                                                        
      function Rmodulo(a,p) 
                                                                        
      if(abs(p).gt.0.) then 
                       if(a*p.gt.0.) then 
                                   ifloor=int(a/p) 
                                   else 
                                   ifloor=int(a/p) -1 
                       endif 
                       Rmodulo= a -float(ifloor)*p 

                       else 
                       Rmodulo=a 
      endif 

      return 
      end
c
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx

      function Imodulo(ia,ip)

      if(abs(ip).gt.0) then
                       if(ia*ip.gt.0) then 
                                   ifloor=ia/ip
                                   else
                                   ifloor=ia/ip -1
                       endif
                       Imodulo= ia -ifloor*ip

                       else
                       Imodulo=ia
      endif

      return
      end 
c
c     xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx0xx
c
c END V2.2 compatibility
c END ROCOTLIB 3.2
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
