!
! XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! XX                                                                  XX
! XX   Robert's Coordinates Transformation Library  ROCOTLIB          XX
! XX                                                                  XX
! XX   Revised  version 1.8, November 2003 (Version for CDPP)         XX
! XX                    ___________________                           XX
! XX                                                                  XX
! XX             Version 1.0 initially supported by                   XX
! XX                                                                  XX
! XX                  EUROPEAN  SPACE  AGENCY                         XX
! XX                                                                  XX
! XX               Study of the Cluster Mission                       XX
! XX                 Planning Related Aspects                         XX
! XX          within the Numerical Simulations Network                XX
! XX                                                                  XX
! XX            Patrick ROBERT, CRPE, November 1992                   XX
! XX                    ___________________                           XX
! XX                                                                  XX
! XX lowing version developped by the Author:                         XX
! XX                                                                  XX
! XX               version 1.0, November 1992                         XX
! XX               version 1.1, July     1993                         XX
! XX               version 1.2, January  1995                         XX
! XX               version 1.3, July     2000                         XX
! XX               version 1.4, June     2001                         XX
! XX               version 1.5, December 2001                         XX
! XX               version 1.6, Juin     2002 (upgrade IGRF -> 2005)  XX
! XX               version 1.7, December 2002 (Version for CDPP)      XX
! XX               version 1.8, November 2003 (adding tsr2sre)        XX
! XX               version 1.9, March    2004 (compatibility with IDL)XX
! XX               version 2.0, June     2010 (version f90           )XX
! XX                                                                  XX
! XX Copyright 1992, Patrick ROBERT, CNRS-ESA, All Rights reserved    XX
! XX                    ___________________                           XX
! XX                                                                  XX
! XX            For details, see the orginal document untitled:       XX
! XX "CLUSTER Software Tools, part I: Coordinate Transformation LibrarXX
! XX                  Document de travail DT/CRPE/1231                XX
! XX               Patrick Robert, CRPE/TID, Juillet 1993             XX
! XX                                                                  XX
! XX         and recently, for CDPP:                                  XX
! XX         ROCOTLIB: a coordinate Transformation Library            XX
! XX         for Solar-Terrestrial studies                            XX
! XX         Patrick ROBERT, version 1.7 - January 2003,              XX
! XX         Rapport Interne no RI-CETP/02/2003                       XX
! XX                                                                  XX
! XX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!
      subroutine cangrat(ux,uy,uz,vx,vy,vz,angle,ratio)
!
! ----------------------------------------------------------------------
!
! *   Class  : basic compute modules of Rocotlib Software
! *   Object : compute_angle_and_ratio beetween U and V vectors
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : ux,uy,uz
!              vx,vy,vz
!
! *   output : sp=U.V
!              angle=angle beetween U and V (radians)
!              ratio= mod(U)/mod(V)
!
! ----------------------------------------------------------------------
!
!
      real (kind=8) :: u1,u2,u3,v1,v2,v3,dp,ru,rv,cot
!
      u1=dble(ux)
      u2=dble(uy)
      u3=dble(uz)
      v1=dble(vx)
      v2=dble(vy)
      v3=dble(vz)
!
      dp= u1*v1 + u2*v2 + u3*v3
      ru= dsqrt(u1*u1 + u2*u2 + u3*u3)
      rv= dsqrt(v1*v1 + v2*v2 + v3*v3)
      cot=dp/(ru*rv)
      cot=cot -sign(1.d-13,cot)
      ratio=real(ru/rv)
      if(cot.ge.1.d0)  then
                 print*, ' *** Rocotlib/cangrat: cos greater than 1 !!!'
                 print*, '                       angle set to 0.'
                 angle=0.
                 return
                 endif
      if(cot.lt.-1.d0) then
                 print*, ' *** Rocotlib/cangrat: cos less than 1 !!!'
                 print*, '                       angle set to 0.'
                 angle=0.
                 return
                 endif
!
      angle=real(dacos(cot))
!
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdatdoy(idoy,iyear,imonth,iday)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_date_from_day_of_year and for a given year
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : iyear,idoy (idoy=1 for january 1)
! *   output : imonth,iday
!
! ----------------------------------------------------------------------
!
!
      dimension month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
      save month,ily
!
!
      call coleapy(iyear,ily)
!
      if(idoy.gt.365.and.ily.eq.0) then
                print*,'*** Rocotlib/cdatdoy: idoy= ',idoy
                print*,'*** Rocotlib/cdatdoy: iyear is not a leap year'
                print*,'                      no more than 365 day'
      stop 'CDATDOY        *** ABORTED ! IYEAR IS NOT A LEAP YEAR   ***'
                                    endif
!
      if(idoy.lt.1) then
             print*,'*** Rocotlib/cdatdoy: idoy= ',idoy
             print*,'*** Rocotlib/cdatdoy: day of the year must be GT 0'
      stop 'CDATDOY        *** ABORTED ! DAY OF THE YEAR MUST BE >0 ***'
                     endif
!
      if(ily.eq.1) then
                   month(2)=29
                   else
                   month(2)=28
                   endif
!
      m=0
!
      do 10 im=1,12
      mp=m
      m=m+month(im)
      if(idoy.le.m) go to 20
   10 continue
   20 continue
!
      imonth=im
      iday=idoy-mp
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdatj00(jd00,iyear,imonth,iday)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_date_from_julian_day_2000 with jd00=0 for jan. 1
! *   Author : P. Robert, CRPE, 1992
! *   Comment: compute date as year, month, day from julian day 2000
!
! *   input  : jd00 julian day 2000 (0= 1/1/2000)
! *   output : iyear,imonth,iday
!
! ----------------------------------------------------------------------
!
!
      jd50 = jd00 +18262
!
      call cdatj50(jd50,iyear,imonth,iday)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdatj50(jd50,iyear,imonth,iday)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_date_from_julian_day_1950 with jd50=0 for jan. 1
! *   Author : P. Robert, CRPE, 1992
! *   Comment: compute date as year, month, day from julian day 1950
!
! *   input  : jd50  julian day 1950 (0= 1/1/1950)
! *   output : iyear,imonth,iday
!
! ----------------------------------------------------------------------
!
      data ily /-1/
      save ily
!
!
      jd1= -1
!
      do 10 iy=1950,3000
      call coleapy(iy,ily)
      jdp=jd1
      if(ily.eq.1) then
                   jd1=jd1+366
                   else
                   jd1=jd1+365
                   endif
      if(jd1.ge.jd50) then
                       iyear=iy
                       go to 20
                       endif
!
   10 continue
      print*,'*** Rocotlib/cdatj50: jd50= ',jd50
      print*,'*** Rocotlib/cdatj50: jul. day correspond to year GT 3000'
      stop 'CDATJ50        *** ABORTED ! JUL. DAY C. TO YEAR >3000 /***'
!
   20 continue
      jd=jd50-jdp
      call cdatdoy(jd,iy,imonth,iday)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdatwee(iweek,iyear,imonth,iday)
!
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_date_for_first_day_of_week_number
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : iweek
! *   output : iyear,imonth,iday
!
! ----------------------------------------------------------------------
!
!
!  calcul de la date correspondant au premier jour de la semaine
!
      if(iweek.eq.1) then
                     imonth=1
                     iday=1
                     return
                     endif
!
! *** premier lundi de l'annee
!
      do 10 ipl=1,7
      call cdoweek(iyear,1,ipl,idow)
      if(idow.eq.1) go to 20
   10 continue
   20 continue
!
! *** jour de l'annee correspondant au premier jour de la semaine
!
      idoy= (iweek-2)*7 +ipl
      if(ipl.eq.1) idoy=idoy+7
!
! *** date correspondante
!
      call cdatdoy(idoy,iyear,imonth,iday)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdipdir(iyear,idoy,d1,d2,d3)
!
! ----------------------------------------------------------------------
!
! *   Class  : basic compute modules of Rocotlib Software
! *   Object : compute_dipole_direction in GEO system
! *   Author : P. Robert, CRPE, 1992, revised 2002
!
! *   Comment: Compute geodipole axis direction from International
!              Geomagnetic Reference Field (IGRF) models for time
!              interval 1965 to 2005. For time out of interval,
!              computation is made for nearest boundary.
!              Code extracted from geopack, N.A. Tsyganenko, Jan. 5 2001
!
! *   input  :  iyear (1965 - 2005), idoy= day of year (1/1=1)
! *   output :  d1,d2,d3  cartesian dipole components in GEO system
!
! ----------------------------------------------------------------------
!
!
      save iy,id,ipr
!
      data iy,id,ipr /-1,-1,-1/
!
!
      if(iyear.eq.iy.and.idoy.eq.id.and.ipr.eq.1) return
!
      iy=iyear
      id=idoy
!
      if (iy.lt.1965) iy=1965
      if (iy.gt.2005) iy=2005
      if (iy.ne.iyear) then
                       if(ipr.eq.0) print 10, iy
                       ipr=1
                       endif
!
! ***  linear interpolation of the geodipole moment components between
!      the values for the nearest epochs
!      (here g10 has opposite sign to that in IGRF tables)
!
! *   1965-1970
!
      if (iy.lt.1970) then
                      f2=(float(iy)+float(idoy)/365.-1965.)/5.
                      f1=1.-f2
                      g10=30334.*f1+30220.*f2
                      g11=-2119.*f1-2068.*f2
                      h11=5776.*f1+5737.*f2
!
! *   1970-1975
!
      elseif (iy.lt.1975) then
                          f2=(float(iy)+float(idoy)/365.-1970.)/5.
                          f1=1.-f2
                          g10=30220.*f1+30100.*f2
                          g11=-2068.*f1-2013.*f2
                          h11=5737.*f1+5675.*f2
!
! *   1975-1980
!
      elseif (iy.lt.1980) then
                          f2=(float(iy)+float(idoy)/365.-1975.)/5.
                          f1=1.-f2
                          g10=30100.*f1+29992.*f2
                          g11=-2013.*f1-1956.*f2
                          h11=5675.*f1+5604.*f2
!
! *   1980-1985
!
      elseif (iy.lt.1985) then
                          f2=(float(iy)+float(idoy)/365.-1980.)/5.
                          f1=1.-f2
                          g10=29992.*f1+29873.*f2
                          g11=-1956.*f1-1905.*f2
                          h11=5604.*f1+5500.*f2
!
! *   1985-1990
!
      elseif (iy.lt.1990) then
                          f2=(float(iy)+float(idoy)/365.-1985.)/5.
                          f1=1.-f2
                          g10=29873.*f1+29775.*f2
                          g11=-1905.*f1-1848.*f2
                          h11=5500.*f1+5406.*f2
!
! *   1990-1995
!
      elseif (iy.lt.1995) then
                          f2=(float(iy)+float(idoy)/365.-1990.)/5.
                          f1=1.-f2
                          g10=29775.*f1+29682.*f2
                          g11=-1848.*f1-1789.*f2
                          h11=5406.*f1+5318.*f2
!
! *   1995-2000
!
      elseif (iy.lt.2000) then
                          f2=(float(iy)+float(idoy)/365.-1995.)/5.
                          f1=1.-f2
                          g10=29682.*f1+29615.*f2
                          g11=-1789.*f1-1728.*f2
                          h11=5318.*f1+5186.*f2
!
! *   2000-2005
!
!     linear extrapolation beyond 2000 by using secular velocity coeff.
!
                          else
                          dt=float(iy)+float(idoy)/366.-2000.
                          g10=29615.-14.6*dt
                          g11=-1728.+10.7*dt
                          h11=5186.-22.5*dt
                          endif
!
!     now calculate the unit vector ezmag in geo coord. system:
!
!     sin(teta0)*cos(lambda0), sin(teta0)*sin(lambda0), and cos(teta0)
!         st0 * cl0                st0 * sl0                ct0
!
!     sqr is the geodipole moment
!
      sq=g11**2+h11**2
      sqq=sqrt(sq)
      sqr=sqrt(g10**2+sq)
!
      sl0=-h11/sqq
      cl0=-g11/sqq
      st0=sqq/sqr
      ct0=g10/sqr
      stcl=st0*cl0
      stsl=st0*sl0
!
      d1=stcl
      d2=stsl
      d3=ct0
!
   10 format(' * ROCOTLIB/cdipdir: Warning &! dipole direction ', &
             'is computed for year ',i4)
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdoweek(iyear,imonth,iday,idow)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_day_of_the_week
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : iyear,imonth,iday
! *   output : idow (1-7)
!
! ----------------------------------------------------------------------
!
      data jdref,julday /-100000,-100000/
      save jdref,julday
!
!
! *** jour julien du lundi de reference (le 1/1/2001 etait un lundi)
!
      call cjd2000(2001,1,1,jdref)
!
! *** jour julien de la date demandee
!
      call cjd2000(iyear,imonth,iday,julday)
!
! *** calcul du jour de la semaine
!
      idiff= julday-jdref
!
      iweek= idiff/7
      idow= idiff -7*iweek +1
      if(idow.lt.1) idow=idow+7
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdoyear(iyear,imonth,iday,idoy)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_day_of_year with idoy=1 for january 1
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : iyear,imonth,iday   ex: 1990,10,17
! *   output : idoy                ex: 290
!
! ----------------------------------------------------------------------
!
!
      dimension month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
      save month,ily
!
      call coleapy(iyear,ily)
!
      if(ily.eq.1) then
                   month(2)=29
                   else
                   month(2)=28
                   endif
!
      if(imonth.lt.1.or.imonth.gt.12) then
          print*, '*** Rocotlib/cdoyear: imonth = ',imonth
          print*, '*** Rocotlib/cdoyear: imonth must be between 1 or 12'
      stop 'CDOYEAR        *** ABORTED ! IMONTH MUST BE >1 AND < 12 ***'
                                      endif
!
      if(iday.gt.month(imonth)) then
           print*, '*** Rocotlib/cdoyear: iday= ',iday
           print*, '*** Rocotlib/cdoyear: this month has only', &
                                          month(imonth),'days'
      stop 'CDOYEAR        *** ABORTED ! ERROR ON IDAY              ***'
                                endif
!
      idoy=iday
      do 10 i=1,imonth-1
      idoy=idoy+month(i)
   10 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cfrdayn(iday,cday,nbcha)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_french_day_name, ex: 'Lundi' for iday=1
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : iday (1-7, otherwise modulo is done)
! *   output : cday,nchar
!
! ----------------------------------------------------------------------
!
      character*(*) cday
      character*8 days (7)
      integer nbca(7)
      save days,nbca
!
      data days /'Lundi','Mardi','Mercredi','Jeudi','Vendredi', &
                 'Samedi','Dimanche'/
      data nbca /5,5,8,5,8,6,8/
!
      iweek= iday/7
      nday= iday -7*iweek
      if(nday.lt.1) nday=nday+7
!
      cday= days(nday)
      nbcha=nbca(nday)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cfrmonn(imonth,cmonth,nchar)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_French_month_name
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : imonth (1-12)
! *   output : cmonth,nchar
!
! ----------------------------------------------------------------------
!
!
      character*(*) cmonth
      integer nbca(12)
      character*9 cara(12)
      save nbca,cara
!
      data cara /'Janvier','Février','Mars','Avril','Mai','Juin', &
                 'Juillet','Août','Septembre','Octobre','Novembre', &
                 'Décembre'/
!
      data nbca /7,7,4,5,3,4,7,4,9,7,8,8/
!
!
      if(imonth.lt.1.or.imonth.gt.12) then
      print*,'***/ IMONTH=',imonth
             print*,'*** Rocotlib/cfrmonn: month= ',imonth
             print*,'*** Rocotlib/cfrmonn: month must be LE 12 and GE 1'
      stop 'CFRMONN        *** ABORTED ! MONTH MUST BE <1 AND >12   ***'
                                      endif
!
      cmonth=cara(imonth)
      nchar =nbca(imonth)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine chouday(ih,im,is,houday)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_hour_of_day from hours, minutes, seconds
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : ih,im,is
! *   output : houday  decimal hour of the day
!
! ----------------------------------------------------------------------
!
      houday=float(ih)+float(im)/60.+float(is)/3600.
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cjd1950(iyear,imonth,iday,jd50)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_julian_day_1950 with jd50=0 for january 1, 1950
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : iyear,imonth,iday   ex: 1990,10,17
! *   output : jd50
!
! ----------------------------------------------------------------------
!
      data ily,idoy /-1,-1/
      save ily,idoy
!
!
      if(iyear.lt.1950) then
      print*,'***/ IYEAR=',iyear
         print*, '*** Rocotlib/cjd1950: iyear= ',iyear
         print*, '*** Rocotlib/cjd1950: iyear must be grather than 1950'
      stop 'CDJ1950        *** ABORTED ! IYEAR MUST BE > 1950       ***'
                        endif
!
      call cdoyear(iyear,imonth,iday,idoy)
!
      jd50=idoy-1
      do 10 i=1950,iyear-1
      call coleapy(i,ily)
      if(ily.eq.1) then
                   jd50=jd50+366
                   else
                   jd50=jd50+365
                   endif
   10 continue
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cjd2000(iyear,imonth,iday,jd00)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_julian_day_2000 with jd00=0 for january 1, 2000
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : iyear,imonth,iday   ex: 2001,10,17
! *   output : jd00 (may be negative)
!
! ----------------------------------------------------------------------
!
      data jd50 /-1/
      save jd50
!
!
      if(iyear.lt.1950) then
      print*,'***/ IYEAR=',iyear
         print*, '*** Rocotlib/cjd2000: iyear= ',iyear
         print*, '*** Rocotlib/cjd2000: iyear must be grather than 1950'
      stop 'CDJ2000        *** ABORTED ! IYEAR MUST BE > 1950       ***'
                        endif
!
      call cjd1950(iyear,imonth,iday,jd50)
!
      jd00= jd50 - 18262
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cmilday(ih,im,is,ims,milday)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_millisec_of_day from hours, minutes, seconds, ms
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : ih,im,is,ims
! *   output : milday  millisecond of the day
!
! ----------------------------------------------------------------------
!
      milday= ih*3600000 + im*60000 + is*1000 +ims
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cnbdmon(iyear,imonth,nbday)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_number_of_day_of_the_month
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : iyear,imonth (1-12)
! *   output : nbday
!
! ----------------------------------------------------------------------
!
!
      dimension nday(12)
      data nday/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
      save nday,ily
!
      if(imonth.lt.1.or.imonth.gt.12) then
      print*,'***/ IMONTH=',imonth
             print*,'*** Rocotlib/cnbdmon: month= ',imonth
             print*,'*** Rocotlib/cnbdmon: month must be LE 12 and GE 1'
      stop 'CNBDMON        *** ABORTED ! MONTH MUST BE >1 AND <12   ***'
                                      endif
!
      call coleapy(iyear,ily)
!
      if(ily.eq.1) then
                   nday(2)=29
                   else
                   nday(2)=28
                   endif
!
      nbday=nday(imonth)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine coleapy(iyear,ileap)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_leap_year with ileap=1 for leap year, 0 if not
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : iyear (ex: 1980)
! *   output : ileap (1 or 0 if iyear is or not a leap year)
!
! ----------------------------------------------------------------------
!
      if(iyear.lt.1900) then
      print*,'***/ IYEAR=',iyear
                    print*,'*** Rocotlib/coleapy: iyear= ',iyear
                    print*,'*** Rocotlib/coleapy: iyear must be > 1900'
      stop 'COLEAPY        *** ABORTED ! IYEAR MUST BE > 1900       ***'
                        endif
!
      ir=iyear-(iyear/4)*4
      if(ir.eq.0) then
                  ileap=1
                  else
                  ileap=0
                  endif
!
      is=iyear-(iyear/100)*100
      if(is.eq.0) then
                  ir=iyear-(iyear/400)*400
                  if(ir.eq.0) then
                              ileap=1
                              else
                              ileap=0
                              endif
                  else
                  return
                  endif
!
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine csundir(iyear,idoy,ih,im,is,gst,slong,sra,sdec,obliq)
!
! ----------------------------------------------------------------------
!
! *   Class  : basic compute modules of Rocotlib Software
! *   Object : compute_sun_direction in GEI system
! *   Author : C.T. Russel, 1971, rev. P.Robert,1992,2001,2002
!
! *   Comment: calculates four quantities in gei system necessary for
!              coordinate transformations dependent on sun position
!              (and, hence, on universal time and season)
!              Initial code from C.T. Russel, cosmic electro-dynamics,
!              v.2, 184-196, 1971.
!              Adaptation P.Robert, November 1992.
!              Revised and F90 compatibility, P. Robert June 2001.
!              Optimisation of DBLE computations and comments,
!              P. Robert, December 2002
!
! *   input  : iyear : year (1901-2099)
!              idoy : day of the year (1 for january 1)
!              ih,im,is : hours, minutes, seconds U.T.
!
! *   output : gst      greenwich mean sideral time (radians)
!              slong    longitude along ecliptic (radians)
!              sra      right ascension (radians)
!              sdec     declination of the sun (radians)
!              obliq    inclination of Earth's axis (radians)
!
! ----------------------------------------------------------------------
!
      double precision dj,fday
!
!
      if(iyear.lt.1901.or.iyear.gt.2099) then
      print*,'***/ IYEAR=',iyear
      print*,'*** Rocotlib/csundir: year = ',iyear
      print*,'*** Rocotlib/csundir: year must be > 1901 and < 2099'
      stop 'CSUNDIR        *** ABORTED ! YEAR MUST BE >1901 & <2099 ***'
                                         endif
!
      pi=acos(-1.)
      pisd= pi/180.
!
! *** Julian day and greenwich mean sideral time
!
      fday=dble(ih*3600+im*60+is)/86400.d0
      ileap=(iyear-1901)/4
      dj=dble(365*(iyear-1900) +ileap +idoy) -0.5d0 +fday
      gst=real(dmod(279.690983d0 +0.9856473354d0*dj+360.d0*fday +180.d0, &
                     360.d0) )*pisd
!
! *** longitude along ecliptic
!
      vl= real( dmod(279.696678d0 +0.9856473354d0*dj,360.d0) )
      t=real(dj/36525.d0)
      g=real( dmod(358.475845d0 +0.985600267d0*dj, 360.d0) )*pisd
      slong=(vl+(1.91946 -0.004789*t)*sin(g) +0.020094*sin(2.*g))*pisd
!
! *** inclination of Earth's axis
!
      obliq=(23.45229 -0.0130125*t)*pisd
      sob=sin(obliq)
      cob=cos(obliq)
!
!     precession of declination (about 0.0056 deg.)
!
      pre= (0.0055686 - 0.025e-4*t)*pisd
!
! *** declination of the sun
!
      slp=slong -pre
      sind=sob*sin(slp)
      cosd=sqrt(1. -sind**2)
      sc=sind/cosd
      sdec=atan(sc)
!
! *** right ascension of the sun
!
      sra=pi -atan2((cob/sob)*sc, -cos(slp)/cosd)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine csunset(iyear,imon,iday,rlat,rlon,tmer,tris,tset,durd)
!
! ----------------------------------------------------------------------
!
! *   Class  : basic compute modules of Rocotlib Software
! *   Object : compute_sunset_time and others
! *   Author : P. Robert, CRPE, 2001
! *   Comment: Spherical Earth assumed
!
! *   input  : iyear,imon,iday, rlat, rlong (rad)
! *   output : tmer,tris,tset,durd on char*8 format, as '23:42:37'
!              tmer: Sun Meridian time (local noon)
!              tris: Sunrise time
!              tset: Sunset  time
!              durd: Duration of the day
!
!              Accuracy: 10 sec.
!
! ----------------------------------------------------------------------
!
!
      character*(*) tmer,tris,tset,durd
!
      lmer=len(tmer)
      lris=len(tris)
      lset=len(tset)
      ldur=len(durd)
!
      if(lmer.lt.8) then
           print*, '*** Rocotlib/csunset: tmer musr be char**8 variable'
      stop 'CSUNSET        *** ABORTED ! TMER MUST BE CHARACTER**8  ***'
                    endif
!
      if(lris.lt.8) then
           print*, '*** Rocotlib/csunset: tris musr be char**8 variable'
      stop 'CSUNSET        *** ABORTED ! TRIS MUST BE CHARACTER**8  ***'
                    endif
!
      if(lset.lt.8) then
           print*, '*** Rocotlib/csunset: tset musr be char**8 variable'
      stop 'CSUNSET        *** ABORTED ! TSET MUST BE CHARACTER**8  ***'
                    endif
!
      if(ldur.lt.8) then
           print*, '*** Rocotlib/csunset: durd musr be char**8 variable'
      stop 'CSUNSET        *** ABORTED ! DURD MUST BE CHARACTER**8  ***'
                    endif
!
      milris= -1
      milset= -1
      milmer= -1
!
      sdp=0.
      svp=0.
!
! *** recherche des temps par rapport au centre du soleil
!     precision 10 sec
!
      do 10 isec=0,86400-1,10
!
      call ctimmil(isec*1000,ih,im,is,ims)
      call ctimpar(iyear,imon,iday,ih,im,is)
      call gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      call tgeovdh(sxgeo,sygeo,szgeo,rlat,rlon,sv,sd,sh)
!
      if(isec.eq.0) go to 20
!
! *** calcul du midi local, quand D (vers l'est) change de signe
!
      if(sv.gt.0..and.sd*sdp.lt.0.) milmer=isec*1000
!
!
! *** calcul des heures de lever et coucher quand V (vertical sortante)
!     change de signe (unite: millisec du jour)
!
! *   cas ou sv est nul
!
      if(sv.eq.0.) then
                   if(svp.lt.0.) then
                                 milris= isec*1000
                                 else
                                 milset= isec*1000
                                 endif
                    go to 10
                    endif
!
! *   sv et svp sont non nuls
!
      if(sv*svp.lt.0.) then
                       if(sv.gt.0.) then
                                    milris= isec*1000
                                    else
                                    milset= isec*1000
                                    endif
                       endif
   20 continue
      sdp=sd
      svp=sv
   10 continue
!
! *** correction du rapport apparent du soleil et de la refraction
!     a l'horizon. Les valeurs prises sont extraites de l'ouvrage
!     "introduction aux ephemerides astronomiques", publie par le
!     Bureau Des Longitudes.
!     Diametre apparent du soleil: 16'
!     Refraction a l'horizon     : 34'
!
! *   calcul de la vitesse apparente du soleil a l'horizon pour le lever
!
      idt= 60000
!
      call ctimmil(milris-idt,ih,im,is,ims)
      call ctimpar(iyear,imon,iday,ih,im,is)
      call gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      call tgeovdh(sxgeo,sygeo,szgeo,rlat,rlon,sv1,sd1,sh1)
!
      call ctimmil(milris,ih,im,is,ims)
      call ctimpar(iyear,imon,iday,ih,im,is)
      call gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      call tgeovdh(sxgeo,sygeo,szgeo,rlat,rlon,sv2,sd2,sh2)
!
      call cangrat(sv1,sd1,sh1,sv2,sd2,sh2,angle,ratio)
!
! *   calcul de la correction a apporter au lever
!
      angle= angle*(180./acos(-1.))
!
      corr= float(16+34)/60.
      icor= int((corr/angle)*float(idt) +0.5)
      milris= milris -icor
!
! *   idem pour le coucher:
!     calcul de la vitesse apparente du soleil a l'horizon
!
      call ctimmil(milset,ih,im,is,ims)
      call ctimpar(iyear,imon,iday,ih,im,is)
      call gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      call tgeovdh(sxgeo,sygeo,szgeo,rlat,rlon,sv1,sd1,sh1)
!
      call ctimmil(milset+idt,ih,im,is,ims)
      call ctimpar(iyear,imon,iday,ih,im,is)
      call gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      call tgeovdh(sxgeo,sygeo,szgeo,rlat,rlon,sv2,sd2,sh2)
!
      call cangrat(sv1,sd1,sh1,sv2,sd2,sh2,angle,ratio)
!
! *   calcul de la correction a apporter au coucher
!
      angle= angle*(180./acos(-1.))
      corr= float(16+34)/60.
      icor= int((corr/angle)*float(idt) +0.5)
      milset= milset +icor
!
! *** encodage des resultats sous la forme HH:MM:SS
!
      tmer='undefin.'
      tris='undefin.'
      tset='undefin.'
      durd='undefin.'
!
      if(milmer.gt.0) then
                      call ctimmil(milmer,ih,im,is,ims)
                      write(tmer,100) ih,im,is
                      endif
!
      if(milris.gt.0) then
                      call ctimmil(milris,ih,im,is,ims)
                      write(tris,100) ih,im,is
                      endif
!
      if(milset.gt.0) then
                      call ctimmil(milset,ih,im,is,ims)
                      write(tset,100) ih,im,is
                      endif
!
      idiff= milset-milris
      if(idiff.lt.0) idiff=idiff +86400000
!
      call ctimmil(idiff,ih,im,is,ims)
!
      write(durd,200) ih,im,is
!
  100 format(i2.2,':',i2.2,':',i2.2)
  200 format(i2  ,':',i2.2,':',i2.2)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ctimhou(houday,ih,im,is)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_time from decimal hour of the day
! *   Author : P. Robert, CRPE, 1992
!
! *   input  : houday decimal hour of the day
! *   output : ih,im,is
!
! ----------------------------------------------------------------------
!
!
      nj=int(abs(houday)/24.)
      rdech=abs(houday)-float(nj*24)
!
      ih=int(rdech)
      im=int((rdech-float(ih))*60.)
      is=int((rdech-float(ih)-float(im)/60.)*3600. +0.001)
!
      if(is.eq.60) then
                   is=0
                   im=im +1
                   endif
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ctimmil(milday,ih,im,is,ims)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_time from millisec. of the day
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : milday millisecond of the day
! *   output : ih,im,is,ims
!
! ----------------------------------------------------------------------
!
!
      nj= milday/86400000
      mil2= milday-nj*86400000
!
      if(nj.ne.0) then
          print*, '*** Rocotlib/ctimmil: milday=',milday,' > 86 400 000'
          print*, '                      assumed:',mil2
                  endif
!
      ih=  milday/3600000
      im= (milday-ih*3600000)/60000
      is= (milday-ih*3600000-im*60000)/1000
      ims= milday-ih*3600000 -im*60000 -is*1000
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ctimpa2(jd1950,houday)
!
! ----------------------------------------------------------------------
!
! *   Class  : basic compute modules of Rocotlib Software
! *   Object : compute_time_parameters and time-dependent matrix
! *   Author : P. Robert, CRPE, 2001
!
! *   Comment: Prepare all time varying quantities for computations of
!              coordinate transforms of the library.
!              Same as ctimpar, only input arguments are changed.
!
! *   input  : jd1950: Julian day 1950 (0= 1/1/1950)
!              houday: decimal hour of the day (U.T.)
! *   output : in common statements
!
! ----------------------------------------------------------------------
!
      data iyear,imonth,iday /-1,-1,-1/
      data ih,im,is /-1,-1,-1/
      save iyear,imonth,iday, ih,im,is
!
      call cdatj50(jd1950,iyear,imonth,iday)
      call ctimhou(houday,ih,im,is)
!
      call ctimpar(iyear,imonth,iday,ih,im,is)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ctimpa3(jd2000,houday)
!
! ----------------------------------------------------------------------
!
! *   Class  : basic compute modules of Rocotlib Software
! *   Object : compute_time_parameters and time-dependent matrix
! *   Author : P. Robert, CRPE, 2001
!
! *   Comment: Prepare all time varying quantities for computations of
!              coordinate transforms of the library.
!              Same as ctimpar, only input arguments are changed.
!
! *   input  : jd2000: Julian day 2000 (0= 1/1/2000)
!              houday: decimal hour of the day (U.T.)
! *   output : in common statements
!
! ----------------------------------------------------------------------
!
      data iyear,imonth,iday /-1,-1,-1/
      data ih,im,is /-1,-1,-1/
      save iyear,imonth,iday, ih,im,is
!
      call cdatj50(jd2000,iyear,imonth,iday)
      call ctimhou(houday,ih,im,is)
!
      call ctimpar(iyear,imonth,iday,ih,im,is)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ctimpar(iyear,imonth,iday,ih,im,is)
!
! ----------------------------------------------------------------------
!
! *   Class  : basic compute modules of Rocotlib Software
! *   Object : compute_time_parameters and time-dependent matrix
! *   Author : P. Robert, CRPE, 1992
!
! *   Comment: Prepare all time varying quantities for computations of
!              coordinate transforms of the library.
!
! *   input  : iyear,imonth,iday (1901<year<2099)
!              ih,im,is : hours, minutes, seconds U.T.
! *   output : in common statements
!
! ----------------------------------------------------------------------
!
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
!
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3
      common /timp04/ grgs1,grgs2,grgs3,  grge1,grge2,grge3
      common /timp05/ gdps1,gdps2,gdps3
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3
!
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3
!
      data idoy /-1/
      save idoy
!
! *** day of year
!
      call cdoyear(iyear,imonth,iday,idoy)
!
! *** Dipole direction in GEO system
!
      call cdipdir(iyear,idoy,gd1,gd2,gd3)
!
      qd12=sqrt(gd1**2 + gd2**2)
!
! *** Sun direction in GEI (from right ascension and declination),
!     greenwich mean sideral time and longitude along ecliptic
!
      call csundir(iyear,idoy,ih,im,is,gst,slong,srasn,sdecl,obliq)
!
      gs1=cos(srasn)*cos(sdecl)
      gs2=sin(srasn)*cos(sdecl)
      gs3=sin(sdecl)
!
! *** sin and cos of GMST
!
      sgst=sin(gst)
      cgst=cos(gst)
!
! *** ecliptic pole in GEI system
!
      ge1=  0.
      ge2= -sin(obliq)
      ge3=  cos(obliq)
!
! *** direction of the rotation axis of the sun in GEI system
!     (from C.T. Russell, Cosmic Electro-dynamics, V.2, 184-196, 1971.)
!
      rad=  57.29578
!
      ras= -74.0/rad
      dec=  63.8/rad
!
      gr1=cos(ras)*cos(dec)
      gr2=sin(ras)*cos(dec)
      gr3=sin(dec)
!
! *** dipole direction in GEI system
!
      gm1= gd1*cgst - gd2*sgst
      gm2= gd1*sgst + gd2*cgst
      gm3= gd3
!
      qm12=sqrt(gm1**2 + gm2**2)
!
! *** direction of the sun in GEO system
!
      ps1=  gs1*cgst + gs2*sgst
      ps2= -gs1*sgst + gs2*cgst
      ps3=  gs3
!
! *** direction of the ecliptic in GEO system
!
      pe1=  ge1*cgst + ge2*sgst
      pe2= -ge1*sgst + ge2*cgst
      pe3=  ge3
!
! *** direction of the rotation axis of the sun in GEO system
!
      pr1=  gr1*cgst + gr2*sgst
      pr2= -gr1*sgst + gr2*cgst
      pr3=  gr3
!
! *** cross product MxS in GEI system
!
      gmgs1= gm2*gs3 - gm3*gs2
      gmgs2= gm3*gs1 - gm1*gs3
      gmgs3= gm1*gs2 - gm2*gs1
!
      rgmgs=sqrt(gmgs1**2 + gmgs2**2 + gmgs3**2)
!
! *** cross product ExS in GEI system
!
      gegs1= ge2*gs3 - ge3*gs2
      gegs2= ge3*gs1 - ge1*gs3
      gegs3= ge1*gs2 - ge2*gs1
!
! *** cross product RxS in GEI system
!
      grgs1= gr2*gs3 - gr3*gs2
      grgs2= gr3*gs1 - gr1*gs3
      grgs3= gr1*gs2 - gr2*gs1
!
      rgrgs=sqrt(grgs1**2 + grgs2**2 + grgs3**2)
!
! *** cross product RxE in GEI system
!
      grge1= gr2*ge3 - gr3*ge2
      grge2= gr3*ge1 - gr1*ge3
      grge3= gr1*ge2 - gr2*ge1
!
! *** cross product DxS in GEO system
!
      gdps1= gd2*ps3 - gd3*ps2
      gdps2= gd3*ps1 - gd1*ps3
      gdps3= gd1*ps2 - gd2*ps1
!
      rgdps=sqrt(gdps1**2 + gdps2**2 + gdps3**2)
!
! *** cross product ExS in GEO system
!
      peps1= pe2*ps3 - pe3*ps2
      peps2= pe3*ps1 - pe1*ps3
      peps3= pe1*ps2 - pe2*ps1
!
! *** cross product RxS in GEO system
!
      prps1= pr2*ps3 - pr3*ps2
      prps2= pr3*ps1 - pr1*ps3
      prps3= pr1*ps2 - pr2*ps1
!
      rprps=sqrt(prps1**2 + prps2**2 + prps3**2)
!
! *** computation of gei to mag vectors
!
      xeima1=  gm1*gm3/qm12
      xeima2=  gm2*gm3/qm12
      xeima3= -qm12
!
      yeima1= -gm2/qm12
      yeima2=  gm1/qm12
      yeima3=  0.
!
! *** computation of gei to sm vectors
!
      yeism1=gmgs1/rgmgs
      yeism2=gmgs2/rgmgs
      yeism3=gmgs3/rgmgs
!
      xeism1= yeism2*gm3 - yeism3*gm2
      xeism2= yeism3*gm1 - yeism1*gm3
      xeism3= yeism1*gm2 - yeism2*gm1
!
! *** computation of gei to gsm vectors
!
      yeigm1= gmgs1/rgmgs
      yeigm2= gmgs2/rgmgs
      yeigm3= gmgs3/rgmgs
!
      zeigm1= gs2*yeigm3 - gs3*yeigm2
      zeigm2= gs3*yeigm1 - gs1*yeigm3
      zeigm3= gs1*yeigm2 - gs2*yeigm1
!
! *** computation of gei to gseq vectors
!
      yeigq1= grgs1/rgrgs
      yeigq2= grgs2/rgrgs
      yeigq3= grgs3/rgrgs
!
      zeigq1= gs2*yeigq3 - gs3*yeigq2
      zeigq2= gs3*yeigq1 - gs1*yeigq3
      zeigq3= gs1*yeigq2 - gs2*yeigq1
!
! *** computation of tetq angle
!
      stetq= (grge1*gs1 + grge2*gs2 + grge3*gs3)/rgrgs
      ctetq= sqrt(1.-stetq**2)
!
! *** computation of mu angle
!
      smu= ps1*gd1 + ps2*gd2 + ps3*gd3
      cmu= sqrt(1.-smu*smu)
!
! *** computation of dzeta angle
!
      cdze= (ge1*gm1   + ge2*gm2   + ge3*gm3)/rgmgs
      sdze= (ge1*gmgs1 + ge2*gmgs2 + ge3*gmgs3)/rgmgs
      epsi=1.e-6
      if(abs(sdze**2 +cdze**2-1.).gt.epsi) then
      stop 'CTIMPAR        *** ABORTED ! ERROR 3                    ***'
      endif
!
! *** computation of phi angle
!
      q=qd12*rgdps
!
      cphi=(gd1*gd3*ps1 + gd2*gd3*ps2 - (gd1**2+gd2**2)*ps3)/q
      sphi=(gd2*ps1 - gd1*ps2)/q
      if(abs(sphi**2 +cphi**2-1.).gt.epsi) then
      stop 'CTIMPAR        *** ABORTED ! ERROR 4                    ***'
      endif
!
! *** computation of geo to mag vectors
!
      yeoma1= -gd2/qd12
      yeoma2=  gd1/qd12
      yeoma3=  0.
!
      xeoma1=  yeoma2*gd3
      xeoma2= -yeoma1*gd3
      xeoma3=  yeoma1*gd2 - yeoma2*gd1
!
! *** computation of geo to sm vectors
!
      yeosm1= gdps1/rgdps
      yeosm2= gdps2/rgdps
      yeosm3= gdps3/rgdps
!
      xeosm1= yeosm2*gd3 - yeosm3*gd2
      xeosm2= yeosm3*gd1 - yeosm1*gd3
      xeosm3= yeosm1*gd2 - yeosm2*gd1
!
! *** computation of geo to gsm vectors
!
      yeogm1=yeosm1
      yeogm2=yeosm2
      yeogm3=yeosm3
!
      zeogm1= ps2*yeogm3 - ps3*yeogm2
      zeogm2= ps3*yeogm1 - ps1*yeogm3
      zeogm3= ps1*yeogm2 - ps2*yeogm1
!
! *** computation of geo to gsq vectors
!
      yeogq1= prps1/rprps
      yeogq2= prps2/rprps
      yeogq3= prps3/rprps
!
      zeogq1= ps2*yeogq3 - ps3*yeogq2
      zeogq2= ps3*yeogq1 - ps1*yeogq3
      zeogq3= ps1*yeogq2 - ps2*yeogq1
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cusdayn(iday,cday,nbcha)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_US_day_name, ex: 'Monday' for iday=1
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : iday (1-7, otherwise modulo is done)
! *   output : cday,nchar
!
! ----------------------------------------------------------------------
!
      character*(*) cday
      character*9 days (7)
      integer nbca(7)
      save days,nbca
!
      data days /'Monday','Tuesday','Wednesday','Thurday','Friday', &
                 'Saturday','Sunday'/
      data nbca /6,7,9,7,6,8,6/
!
      iweek= iday/7
      nday= iday -7*iweek
      if(nday.lt.1) nday=nday+7
!
      cday= days(nday)
      nbcha=nbca(nday)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cusmonn(imonth,cmonth,nchar)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_US_month_name
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : imonth (1-12)
! *   output : cmonth,nchar
!
! ----------------------------------------------------------------------
!
!
      character*(*) cmonth
      integer nbca(12)
      character*9 cara(12)
      save nbca,cara
!
      data cara/'January','February','Mars','April','May','June','July', &
                'August','September','October','November','December'/
      data nbca /7,8,4,5,3,4,4,6,9,7,8,8/
!
!
      if(imonth.lt.1.or.imonth.gt.12) then
      print*,'***/ IMONTH=',imonth
             print*,'*** Rocotlib/cusmonn: month= ',imonth
             print*,'*** Rocotlib/cusmonn: month must be LE 12 and GE 1'
      stop 'CUSMONN        *** ABORTED ! MONTH MUST BE >1 AND <12   ***'
                                      endif
!
      cmonth=cara(imonth)
      nchar =nbca(imonth)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cweedoy(iyear,imonth,iday,iweek)
!
! ----------------------------------------------------------------------
!
! *   Class  : calendar modules of Rocotlib Software
! *   Object : compute_week_of_the_year
! *   Author : P. Robert, CRPE, 2001
!
! *   input  : iyear,imonth,iday
! *   output : iweek
!
! ----------------------------------------------------------------------
!
      data idoty1,idoty /-1,-1/
      save idoty1,idoty
!
!
! *** calcul du 1er lundi de l'annee
!
      do 10 id=1,7
      call cdoweek(iyear,1,id,idow)
      if(idow.eq.1) then
                    idl=id
                    go to 20
                    endif
   10 continue
      print*, '*** Rocotlib/cweedoy: error in id1 computation'
      stop 'CWEEDOY        *** ABORTED ! ERROR ON ID1 COMPUTATION   ***'
   20 continue
!
      if(idl.eq.1) then
                   iw1=1
                   else
                   iw1=2
                   endif
!
! *** jour de l'annee du 1er lundi
!
      call cdoyear(iyear,1,idl,idoty1)
!
! *** jour de l'annee de la date demandee
!
      call cdoyear(iyear,imonth,iday,idoty)
!
! *** calcul de la semaine de l'annee
!
      idiff= idoty-idoty1
!
      if(idiff.ge.0) then
                     iweek= idiff/7 +iw1
                     else
                     iweek=1
                     endif
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gdipdir(dxgei,dygei,dzgei,dxgeo,dygeo,dzgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : give modules of Rocotlib Software
! *   Object : give_dipole_direction in GEI and GEO system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: values are extracted from common
!
! *   input  : none
! *   output : dxgei,dygei,dzgei  cartesian dipole GEI coordinates
!              dxgeo,dygeo,dzgeo  cartesian dipole GEO coordinates
!
! ----------------------------------------------------------------------
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
!
      dxgei=gm1
      dygei=gm2
      dzgei=gm3
!
      dxgeo=gd1
      dygeo=gd2
      dzgeo=gd3
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gdiptan(diptan)
!
! ----------------------------------------------------------------------
!
! *   Class  : give modules of Rocotlib Software
! *   Object : give_dipole_tilt_angle in radians
! *   Author : P. Robert, CRPE, 1992
! *   Comment: values are extracted from common
!
! *   input  : none
! *   output : diptan (radians)
!
! ----------------------------------------------------------------------
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      diptan=atan2(smu,cmu)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gecldir(exgei,eygei,ezgei,exgeo,eygeo,ezgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : give modules of Rocotlib Software
! *   Object : give_ecliptic_direction in GEI and GEO system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: values are extracted from common
!
! *   input  : none
! *   output : exgei,eygei,ezgei  cartesian ecliptic GEI coordinates
!              exgeo,eygeo,ezgeo  cartesian ecliptic GEO coordinates
!
! ----------------------------------------------------------------------
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
!
      exgei=ge1
      eygei=ge2
      ezgei=ge3
!
      exgeo=pe1
      eygeo=pe2
      ezgeo=pe3
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gsrodir(rxgei,rygei,rzgei,rxgeo,rygeo,rzgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : give modules of Rocotlib Software
! *   Object : give_sun_rotation_direction in GEI and GEO system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: values are extracted from common
!
! *   input  : none
! *   output : rxgei,rygei,rzgei cartesian sun rotation GEI coordinates
!              rxgeo,rygeo,rzgeo cartesian sun rotation GEO coordinates
!
! ----------------------------------------------------------------------
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
!
      rxgei=gr1
      rygei=gr2
      rzgei=gr3
!
      rxgeo=pr1
      rygeo=pr2
      rzgeo=pr3
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gsundir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : give modules of Rocotlib Software
! *   Object : give_sun_direction in GEI and GEO system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: values are extracted from common
!
! *   input  : none
! *   output : sxgei,sygei,szgei  cartesian sun GEI coordinates
!              sxgeo,sygeo,szgeo  cartesian sun GEO coordinates
!
! ----------------------------------------------------------------------
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
!
      sxgei=gs1
      sygei=gs2
      szgei=gs3
!
      sxgeo=ps1
      sygeo=ps2
      szgeo=ps3
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gsunpar(gmst,slon,sras,sdec,obli)
!
! ----------------------------------------------------------------------
!
! *   Class  : give modules of Rocotlib Software
! *   Object : give_sun_parameter dependant of time in GEI system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: values are extracted from common
!
! *   input  : none
! *   output : gmst   greenwich mean sideral time (radians)
!              slon   longitude along ecliptic (radians)
!              sras   right ascension (radians)
!              sdec   declination of the sun (radians)
!
! ----------------------------------------------------------------------
!
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq
!
      gmst= gst
      slon= slong
      sras= srasn
      sdec= sdecl
      obli= obliq
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine gvernum(vernum,verdat)
!
      character*14 verdat
!
! ----------------------------------------------------------------------
!
! *   Class  : give modules of Rocotlib Software
! *   Object : give_version_number and modification date of the library
! *   Author : P. Robert, CRPE, 1992
! *   Comment: values are extracted from common
!
! *   input  : none
! *   output : vernum (ex 1.0) and verdat (ex: 'January 1995')
!
! ----------------------------------------------------------------------
!
!     version 1.0, November 1992
!     version 1.1, July     1993
!     version 1.2, January  1995
!     version 1.3, July     2000
!     version 1.4, June     2001
!     version 1.5, December 2001
!     version 1.6, Juin     2002 (upgrade IGRF -> 2005)
!     version 1.7, December 2002 (Version for CDPP)
!     version 1.8, November 2003 (adding tsr_sr2)
!     version 1.9, March    2004 (compatibility with V 1.9 IDL)
!     version 2.0, June     2010 (version f90           )
!
      vernum=2.0
      verdat='June 2010'
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine plibinf
!
      character*14 verdat
!
! ----------------------------------------------------------------------
!
! *   Class  : print modules of Rocotlib Software
! *   Object : print_library_informations
! *   Author : P. Robert, CRPE, 1992
! *   Comment: could be a main program
!
! *   input  : none
! *   output : none; print infos on output
!
! ----------------------------------------------------------------------
!
      call gvernum(vernum,verdat)
!
!
      print 100, ' '
      print 100, '*****************************************************'
      print 100, ' '
      print 100, '        Coordinates Transformation Library  ROCOTLIB'
      print 100, ' '
      print 200, '              Revised Version ',vernum, ' - ',verdat
      print 100, ' '
      print 100, '                  ___________________'
      print 100, ' '
      print 100, '                initially supported by '
      print 100, ' '
      print 100, '                EUROPEAN  SPACE  AGENCY'
      print 100, ' '
      print 100, '             Study of the Cluster Mission'
      print 100, '               Planning Related Aspects'
      print 100, '        within the Numerical Simulations Network'
      print 100, ' '
      print 100, '          Patrick ROBERT, CRPE, November 1992'
      print 100, '                  ___________________'
      print 100, ' '
      print 100, '         First  version 1.0, PR, CRPE, November 1992'
      print 100, '        revised version 1.1, PR, CRPE, July     1993'
      print 100, '        revised version 1.2, PR, CETP, January  1995'
      print 100, '        revised version 1.3, PR, CETP, July     2000'
      print 100, '        revised version 1.4, PR, CETP, June     2001'
      print 100, '        revised version 1.5, PR, CETP, December 2001'
      print 100, '        revised version 1.6, PR, CETP, June     2002'
      print 100, '        revised version 1.7, PR, CETP, December 2002'
      print 100, '        revised version 1.8, PR, CETP, November 2003'
      print 100, '        revised version 1.9, PR, CETP, March    2004'
      print 100, ' '
      print 100, ' '
      print 100, 'Copyright 1992, Patrick ROBERT, CNRS-ESA, All Rights', &
                 ' reserved'
      print 100, '                  ___________________'
      print 100, ' '
      print 100, '      For details, see the orginal document untitled:'
      print 100, '"CLUSTER Software Tools, 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, '                  ___________________'
      print 100, ' '
      print 100, '*****************************************************'
      print 100, ' '
!
  100 format(3a)
  200 format(a,f4.1,2a)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine recoor(x,y,z,cs)
!
! ----------------------------------------------------------------------
!
! *   Class  : read modules of Rocotlib Software
! *   Object : read coordinate values from input
! *   Author : P. Robert, CRPE, 2002
!
! *   Comment: read cs and x,y,z cartesian  or spherical coordinates
! *            print error if cs is not valid, and ask again
!
! *   input  : cs (c or s) and x,y,z on standard input
! *   output : cs and x,y,z always in cartesian coordinates
!
! ----------------------------------------------------------------------
!
      character*1 cs
!
!
      pisd=acos(-1.)/180.
!
   10 continue
!
      print 100, 'input coordinates: cartesian (c) or spherical (s)'
      read *, cs
      print 100, cs
!
      if(cs.ne.'s'.and.cs.ne.'c') then
                                  print*, 'only c or s'
                                  print*, 'again...'
                                  go to 10
                                  endif
      if(cs.eq.'s') then
                    print 100, 'R,tet, phi (RE,deg.) ?'
                    read *,  r,tet, phi
                    print 110,  r,tet, phi
!
                    tet=tet*pisd
                    phi=phi*pisd
                    call tsphcar(r,tet,phi,x,y,z)
                    print 120, 'then x,y,z = ',x,y,z,'(RE)'
!
                    else
                    print*, 'x,y,z ? (RE) '
                    read *,  x,y,z
                    print*,  x,y,z
!
                    call tcarsph(x,y,z,r,tet,phi)
                    tet=tet/pisd
                    phi=phi/pisd
                    print 200, 'then R,tet,phi= ',r,tet,phi,' (RE,deg.)'
                    endif
!
  100 format(a)
  110 format(e14.6,1x,f8.3,1x,f8.3)
  120 format(1x,a,3e11.3,a)
  200 format(1x,a,e11.3,f8.3,f8.3,a)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine recsys(csys)
!
! ----------------------------------------------------------------------
!
! *   Class  : read modules of Rocotlib Software
! *   Object : read coordinate system from input and check validity
! *   Author : P. Robert, CRPE, 2002
!
! *   Comment: read csys string variable and check validity
! *            (only gei, geo, mag, sma, gsm, gse, gsq)
! *            print error if csys is not valid, and ask again
!
! *   input  : csys on standard input
! *   output : csys
!
! ----------------------------------------------------------------------
!
      character*3 csys
!
   10 continue
!
      print 100, 'repere ? (gei, geo, mag, sma, gsm, gse, gsq)'
      read *,csys
      print 100,csys
  100 format(a)
!
      if(csys.eq.'gei') return
      if(csys.eq.'geo') return
      if(csys.eq.'mag') return
      if(csys.eq.'sma') return
      if(csys.eq.'gsm') return
      if(csys.eq.'gse') return
      if(csys.eq.'gsq') return
!
      print*, 'Only gei, geo, mag, sma, gsm, gse or gsq please...'
      print*, 'again...'
      go to 10
!
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine redate(iyear,imonth,iday)
!
! ----------------------------------------------------------------------
!
! *   Class  : read modules of Rocotlib Software
! *   Object : read_date from input and check validity
! *   Author : P. Robert, CRPE, 1992
!
! *   Comment: test if year is gt 1900
! *            test if imonth is not greater than 12
!              test if iday is not greather than lengh of the month,
!                        takink into account the leap years.
!               print error if date is not valid, and ask again
!
! *   input  : iyear,imonth,iday given in standard input
! *   output : iyear,imonth,iday
!
! ----------------------------------------------------------------------
!
!
      dimension month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      data ily /-1/
      save month,ily
!
   10 continue
      print 100, ' iyear,imonth,iday ? (ex: 1990,10,17)'
      read *,  iyear,imonth,iday
      print 110,  iyear,imonth,iday
!
  100 format(a)
  110 format(1x,i4,1x,i2,1x,i2)
!
      if(iyear.lt.1900) then
         print*, '*** Rocotlib/redate: iyear must be greater than 1900'
         print*, '                     again ...'
                        go to 10
                        endif
!
      call coleapy(iyear,ily)
!
!
      if(ily.eq.1) then
                   month(2)=29
                   else
                   month(2)=28
                   endif
!
      if(imonth.lt.1.or.imonth.gt.12) then
          print*, '*** Rocotlib/redate: imonth must be between 1 or 12 '
          print*, '                     again...'
                     go to 10
                                      endif
!
      if(iday.gt.month(imonth)) then
                     print*, '*** Rocotlib/redate: this month has only', &
                                  month(imonth),'days'
                     print*, '                     again...'
                go to 10
                                endif
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine retime(ih,im,is)
!
! ----------------------------------------------------------------------
!
! *   Class  : read modules of Rocotlib Software
! *   Object : read_time from input and check validity
! *   Author : P. Robert, CRPE, 1992
! *   Comment: read hour, minute and second and verifie validity
!              ih must be between 1 and 23, im and is between 1 and 59
!               print error if time is not valid, and ask again
!
! *   input  : ih,im,is on standard input
! *   output : ih,im,is
!
! ----------------------------------------------------------------------
!
   10 continue
      print 100, ' hour, minute, second ? (ex: 10,45,50)'
      read *, ih,im,is
      print 110, ih,im,is
!
  100 format(a)
  110 format(1x,i2,1x,i2,1x,i2)
!
      if(ih.lt.0.or.ih.gt.23) then
                              print*, 'hour between 0 and 23 please'
                              print*, 'again...'
                              go to 10
                              endif
!
      if(im.lt.0.or.im.gt.59) then
                              print*, 'minute between 0 and 59 please'
                              print*, 'again...'
                              go to 10
                              endif
!
      if(is.lt.0.or.is.gt.59) then
                              print*, 'second between 0 and 59 please'
                              print*, 'again...'
                              go to 10
                              endif
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tcarsph(x,y,z,r,teta,phi)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_car_to_sph: CAR -> SPH  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: none
!
! *   input  :   x,y,z      cartesian coordinates
! *   output :   r,teta,phi spherical coordinates (angles in radians)
!
! ----------------------------------------------------------------------
!
!
      teta=0.
      phi=0.
      sq=x**2+y**2
      r=sqrt(sq+z**2)
!
      if(r.lt.1.e-30) return
!
! *** en dessous de 1/10000 degres, on considere que teta =0 ou 180
!     et phi indefini, mis a zero
!
      if(sq/r.gt.1.7e-6) then
                   phi=atan2(y,x)
                   if(phi.lt.-179.999.and.phi.gt.-180.) phi=180.
                   if(phi.lt. 0.0001.and.phi.gt.-0.0001) phi=0.
                   teta=acos(z/r)
                   return
                   endif
!
      if (z.lt.0.) then
                   teta=acos(-1.)
                   else
                   teta=0.
                   endif
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tdmegeo(xdme,ydme,zdme,rlat,rlong,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_dme_to_geo: DM  -> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xdme,ydme,zdme   cartesian dme coordinates
!              rlat,rlong       latitude and longitude of the point
!                               of observation (radians)
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
!
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
!
      y1= real(dprod(gd2,r3) - dprod(gd3,r2))
      y2= real(dprod(gd3,r1) - dprod(gd1,r3))
      y3= real(dprod(gd1,r2) - dprod(gd2,r1))
!
      q=sqrt(y1**2 + y2**2 + y3**2)
!
      y1=y1/q
      y2=y2/q
      y3=y3/q
!
      x1= real(dprod(y2,gd3) - dprod(y3,gd2))
      x2= real(dprod(y3,gd1) - dprod(y1,gd3))
      x3= real(dprod(y1,gd2) - dprod(y2,gd1))
!
!
      xgeo= x1*xdme + y1*ydme + gd1*zdme
      ygeo= x2*xdme + y2*ydme + gd2*zdme
      zgeo= x3*xdme + y3*ydme + gd3*zdme
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeigeo(xgei,ygei,zgei,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gei_to_geo: GEI -> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgei,ygei,zgei   cartesian gei coordinates
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq
!
      xgeo=  cgst*xgei + sgst*ygei
      ygeo= -sgst*xgei + cgst*ygei
      zgeo=  zgei
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeigse(xgei,ygei,zgei,xgse,ygse,zgse)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gei_to_gse: GEI -> GSE  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgei,ygei,zgei cartesian gei coordinates
! *   output : xgse,ygse,zgse cartesian gse coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3
!
      xgse=   gs1*xgei +   gs2*ygei +   gs3*zgei
      ygse= gegs1*xgei + gegs2*ygei + gegs3*zgei
      zgse=   ge1*xgei +   ge2*ygei +   ge3*zgei
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeigsm(xgei,ygei,zgei,xgsm,ygsm,zgsm)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gei_to_gsm: GEI -> GSM  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgei,ygei,zgei   cartesian gei coordinates
! *   output : xgsm,ygsm,zgsm   cartesian gsm coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3
!
      xgsm=    gs1*xgei +    gs2*ygei +    gs3*zgei
      ygsm= yeigm1*xgei + yeigm2*ygei + yeigm3*zgei
      zgsm= zeigm1*xgei + zeigm2*ygei + zeigm3*zgei
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeigsq(xgei,ygei,zgei,xgsq,ygsq,zgsq)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gei_to_gsq: GEI -> GSEQ system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgei,ygei,zgei   cartesian gei coordinates
! *   output : xgsq,ygsq,zgsq   cartesian gsq coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3
!
      xgsq=    gs1*xgei +    gs2*ygei +    gs3*zgei
      ygsq= yeigq1*xgei + yeigq2*ygei + yeigq3*zgei
      zgsq= zeigq1*xgei + zeigq2*ygei + zeigq3*zgei
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeimag(xgei,ygei,zgei,xmag,ymag,zmag)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gei_to_mag: GEI -> MAG  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgei,ygei,zgei   cartesian gei coordinates
! *   output : xmag,ymag,zmag   cartesian mag coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3
!
      xmag= xeima1*xgei + xeima2*ygei + xeima3*zgei
      ymag= yeima1*xgei + yeima2*ygei + yeima3*zgei
      zmag=    gm1*xgei +    gm2*ygei +    gm3*zgei
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeisma(xgei,ygei,zgei,xsma,ysma,zsma)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gei_to_sma: GEI -> SM   system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgei,ygei,zgei   cartesian gei coordinates
! *   output : xsma,ysma,zsma   cartesian sma coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3
!
      xsma= xeism1*xgei + xeism2*ygei + xeism3*zgei
      ysma= yeism1*xgei + yeism2*ygei + yeism3*zgei
      zsma=    gm1*xgei +    gm2*ygei +    gm3*zgei
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeodme(xgeo,ygeo,zgeo,rlat,rlong,xdme,ydme,zdme)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_dme: GEO -> DM   system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgeo,ygeo,zgeo   cartesian geo coordinates
!              rlat,rlong       latitude and longitude of the point
!                               of observation (radians)
! *   output : xdme,ydme,zdme   cartesian dme coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
!
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
!
      y1= real(dprod(gd2,r3) - dprod(gd3,r2))
      y2= real(dprod(gd3,r1) - dprod(gd1,r3))
      y3= real(dprod(gd1,r2) - dprod(gd2,r1))
!
      q=sqrt(y1**2 + y2**2 + y3**2)
!
      y1=y1/q
      y2=y2/q
      y3=y3/q
!
      x1= real(dprod(y2,gd3) - dprod(y3,gd2))
      x2= real(dprod(y3,gd1) - dprod(y1,gd3))
      x3= real(dprod(y1,gd2) - dprod(y2,gd1))
!
      xdme=  x1*xgeo +  x2*ygeo +  x3*zgeo
      ydme=  y1*xgeo +  y2*ygeo +  y3*zgeo
      zdme= gd1*xgeo + gd2*ygeo + gd3*zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeogei(xgeo,ygeo,zgeo,xgei,ygei,zgei)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_gei: GEO -> GEI  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgeo,ygeo,zgeo cartesian geo coordinates
! *   output : xgei,ygei,zgei cartesian gei coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp00/ sgst,cgst, gst,slong,srasn,sdecl,obliq
!
      xgei= cgst*xgeo - sgst*ygeo
      ygei= sgst*xgeo + cgst*ygeo
      zgei= zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeogse(xgeo,ygeo,zgeo,xgse,ygse,zgse)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_gse: GEO -> GSE  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgeo,ygeo,zgeo   cartesian geo coordinates
! *   output : xgse,ygse,zgse   cartesian gse coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3
!
      xgse=   ps1*xgeo +   ps2*ygeo +   ps3*zgeo
      ygse= peps1*xgeo + peps2*ygeo + peps3*zgeo
      zgse=   pe1*xgeo +   pe2*ygeo +   pe3*zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeogsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_gsm: GEO -> GSM  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgeo,ygeo,zgeo   cartesian geo coordinates
! *   output : xgsm,ygsm,zgsm   cartesian gsm coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3
!
      xgsm=    ps1*xgeo +    ps2*ygeo +    ps3*zgeo
      ygsm= yeogm1*xgeo + yeogm2*ygeo + yeogm3*zgeo
      zgsm= zeogm1*xgeo + zeogm2*ygeo + zeogm3*zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeogsq(xgeo,ygeo,zgeo,xgsq,ygsq,zgsq)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_gsq: GEO -> GSEQ system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgeo,ygeo,zgeo   cartesian geo coordinates
! *   output : xgsq,ygsq,zgsq   cartesian gsq coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3
!
      xgsq=    ps1*xgeo +    ps2*ygeo +    ps3*zgeo
      ygsq= yeogq1*xgeo + yeogq2*ygeo + yeogq3*zgeo
      zgsq= zeogq1*xgeo + zeogq2*ygeo + zeogq3*zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeomag(xgeo,ygeo,zgeo,xmag,ymag,zmag)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_mag: GEO -> MAG  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgeo,ygeo,zgeo   cartesian geo coordinates
! *   output : xmag,ymag,zmag   cartesian mag coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3
!
      xmag= xeoma1*xgeo + xeoma2*ygeo + xeoma3*zgeo
      ymag= yeoma1*xgeo + yeoma2*ygeo + yeoma3*zgeo
      zmag=    gd1*xgeo +    gd2*ygeo +    gd3*zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeosma(xgeo,ygeo,zgeo,xsma,ysma,zsma)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_sma: GEO -> SM   system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgeo,ygeo,zgeo   cartesian geo coordinates
! *   output : xsma,ysma,zsma   cartesian sma coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3
!
      xsma= xeosm1*xgeo + xeosm2*ygeo + xeosm3*zgeo
      ysma= yeosm1*xgeo + yeosm2*ygeo + yeosm3*zgeo
      zsma=    gd1*xgeo +    gd2*ygeo +    gd3*zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgeovdh(xgeo,ygeo,zgeo,rlat,rlong,xvdh,yvdh,zvdh)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_geo_to_vdh: GEO -> VDH  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: local system, non time dependent
!
! *   input  : xgeo,ygeo,zgeo   cartesian geo coordinates
!              rlat,rlong       latitude and longitude of the point
!                               of observation (radians)
! *   output : xvdh,yvdh,zvdh   cartesian vdh coordinates
!
! ----------------------------------------------------------------------
!
!
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
!
      v1=r1
      v2=r2
      v3=r3
!
      q12=sqrt(r1**2 + r2**2)
!
      d1= -r2/q12
      d2=  r1/q12
      d3=  0.
!
      h1= -r1*r3/q12
      h2= -r2*r3/q12
      h3=  q12
!
      xvdh= v1*xgeo + v2*ygeo + v3*zgeo
      yvdh= d1*xgeo + d2*ygeo + d3*zgeo
      zvdh= h1*xgeo + h2*ygeo + h3*zgeo
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsegei(xgse,ygse,zgse,xgei,ygei,zgei)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gse_to_gei: GSE -> GEI  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgse,ygse,zgse    cartesian gse coordinates
! *   output : xgei,ygei,zgei    cartesian gei coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp03/ gmgs1,gmgs2,gmgs3,  gegs1,gegs2,gegs3
!
      xgei= gs1*xgse + gegs1*ygse + ge1*zgse
      ygei= gs2*xgse + gegs2*ygse + ge2*zgse
      zgei= gs3*xgse + gegs3*ygse + ge3*zgse
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsegeo(xgse,ygse,zgse,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gse_to_geo: GSE -> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgse,ygse,zgse   cartesian gse coordinates
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp06/ peps1,peps2,peps3,  prps1,prps2,prps3
!
      xgeo= ps1*xgse + peps1*ygse + pe1*zgse
      ygeo= ps2*xgse + peps2*ygse + pe2*zgse
      zgeo= ps3*xgse + peps3*ygse + pe3*zgse
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsegsm(xgse,ygse,zgse,xgsm,ygsm,zgsm)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gse_to_gsm: GSE -> GSM  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgse,ygse,zgse   cartesian gse coordinates
! *   output : xgsm,ygsm,zgsm   cartesian gsm coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xgsm= xgse
      ygsm=  cdze*ygse + sdze*zgse
      zgsm= -sdze*ygse + cdze*zgse
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsegsq(xgse,ygse,zgse,xgsq,ygsq,zgsq)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gse_to_gsq: GSE -> GSEQ system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgse,ygse,zgse   cartesian gse coordinates
! *   output : xgsq,ygsq,zgsq   cartesian gsq coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xgsq= xgse
      ygsq= ctetq*ygse - stetq*zgse
      zgsq= stetq*ygse + ctetq*zgse
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsesr2(xgse,ygse,zgse,rotx,roty,rotz,xsr2,ysr2,zsr2)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gse_to_sr2: GSE -> SR2  system
! *   Author : P. Robert, CETP, 2001
! *   Comment: local system, non time dependent
!
! *   input  : xgse,ygse,zgse cartesian gse coordinates
!              rotx,roty,rotz cartesian gse coordinates of rotation axis
! *   output : xsr2,ysr2,zsr2 cartesian sr2 coordinates
!
! ----------------------------------------------------------------------
!
!
! *** set transform matrix with Spin axis terms
!
      rmod= sqrt(rotx**2 + roty**2 + rotz**2)
!
      rx=rotx/rmod
      ry=roty/rmod
      rz=rotz/rmod
!
      a= 1./sqrt(ry*ry + rz*rz)
!
      x1= (ry*ry + rz*rz)*a
      x2= -rx*ry*a
      x3= -rx*rz*a
!
      y1=  0.
      y2=  rz*a
      y3= -ry*a
!
      z1= rx
      z2= ry
      z3= rz
!
! *** tranform the input vector from gse to sr2
!
      xsr2= x1*xgse + x2*ygse + x3*zgse
      ysr2= y1*xgse + y2*ygse + y3*zgse
      zsr2= z1*xgse + z2*ygse + z3*zgse
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsmgei(xgsm,ygsm,zgsm,xgei,ygei,zgei)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsm_to_gei: GSM -> GEI  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
! *   output : xgei,ygei,zgei   cartesian gei coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp09/ yeigm1,yeigm2,yeigm3,  zeigm1,zeigm2,zeigm3
!
      xgei= gs1*xgsm + yeigm1*ygsm + zeigm1*zgsm
      ygei= gs2*xgsm + yeigm2*ygsm + zeigm2*zgsm
      zgei= gs3*xgsm + yeigm3*ygsm + zeigm3*zgsm
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsmgeo(xgsm,ygsm,zgsm,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsm_to_geo: GSM -> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp14/ yeogm1,yeogm2,yeogm3,  zeogm1,zeogm2,zeogm3
!
      xgeo= ps1*xgsm + yeogm1*ygsm + zeogm1*zgsm
      ygeo= ps2*xgsm + yeogm2*ygsm + zeogm2*zgsm
      zgeo= ps3*xgsm + yeogm3*ygsm + zeogm3*zgsm
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsmgse(xgsm,ygsm,zgsm,xgse,ygse,zgse)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsm_to_gse: GSM -> GSE  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
! *   output : xgse,ygse,zgse   cartesian gse coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xgse= xgsm
      ygse= cdze*ygsm - sdze*zgsm
      zgse= sdze*ygsm + cdze*zgsm
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsmgsq(xgsm,ygsm,zgsm,xgsq,ygsq,zgsq)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsm_to_gsq: GSM -> GSQ  system
! *   Author : P. Robert, CRPE, 2002
!
! *   input  : xgsm,ygsm,zgsm cartesian gsm coordinates
! *   output : xgsq,ygsq,zgsq cartesian gsq coordinates
!
! ----------------------------------------------------------------------
!
!
      call tgsmgse(xgsm,ygsm,zgsm,ax,ay,az)
      call tgsegsq(ax,ay,az,xgsq,ygsq,zgsq)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsmmag(xgsm,ygsm,zgsm,xmag,ymag,zmag)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsm_to_mag: GSM -> MAG  system
! *   Author : P. Robert, CRPE, 2002
!
! *   input  : xgsm,ygsm,zgsm cartesian gsm coordinates
! *   output : xmag,ymag,zmag cartesian mag coordinates
!
! ----------------------------------------------------------------------
!
!
      call tgsmsma(xgsm,ygsm,zgsm,ax,ay,az)
      call tsmamag(ax,ay,az,xmag,ymag,zmag)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsmsma(xgsm,ygsm,zgsm,xsma,ysma,zsma)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsm_to_sma: GSM -> SM   system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgsm,ygsm,zgsm   cartesian gsm coordinates
! *   output : xsma,ysma,zsma   cartesian sma coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xsma= cmu*xgsm - smu*zgsm
      ysma= ygsm
      zsma= smu*xgsm + cmu*zgsm
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsqgei(xgsq,ygsq,zgsq,xgei,ygei,zgei)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsq_to_gei: GSEQ-> GEI  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgsq,ygsq,zgsq   cartesian gsq coordinates
! *   output : xgei,ygei,zgei   cartesian gei coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp10/ yeigq1,yeigq2,yeigq3,  zeigq1,zeigq2,zeigq3
!
      xgei= gs1*xgsq + yeigq1*ygsq + zeigq1*zgsq
      ygei= gs2*xgsq + yeigq2*ygsq + zeigq2*zgsq
      zgei= gs3*xgsq + yeigq3*ygsq + zeigq3*zgsq
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsqgeo(xgsq,ygsq,zgsq,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsq_to_geo: GSEQ-> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgsq,ygsq,zgsq   cartesian gsq coordinates
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp15/ yeogq1,yeogq2,yeogq3,  zeogq1,zeogq2,zeogq3
!
      xgeo= ps1*xgsq + yeogq1*ygsq + zeogq1*zgsq
      ygeo= ps2*xgsq + yeogq2*ygsq + zeogq2*zgsq
      zgeo= ps3*xgsq + yeogq3*ygsq + zeogq3*zgsq
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsqgse(xgsq,ygsq,zgsq,xgse,ygse,zgse)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsq_to_gse: GSEQ-> GSE  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xgsq,ygsq,zgsq   cartesian gsq coordinates
! *   output : xgse,ygse,zgse   cartesian gse coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xgse= xgsq
      ygse= ctetq*ygsq + stetq*zgsq
      zgse=-stetq*ygsq + ctetq*zgsq
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tgsqgsm(xgsq,ygsq,zgsq,xgsm,ygsm,zgsm)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_gsq_to_gsm: GSQ -> GSM  system
! *   Author : P. Robert, CRPE, 2002
!
! *   input  : xgsq,ygsq,zgsq cartesian gsq coordinates
! *   output : xgsm,ygsm,zgsm cartesian gsm coordinates
!
! ----------------------------------------------------------------------
!
!
      call tgsqgse(xgsq,ygsq,zgsq,xx,yy,zz)
      call tgsegsm(xx,yy,zz,xgsm,ygsm,zgsm)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tmaggei(xmag,ymag,zmag,xgei,ygei,zgei)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_mag_to_gei: MAG -> GEI  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xmag,ymag,zmag cartesian mag coordinates
! *   output : xgei,ygei,zgei cartesian gei coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp07/ xeima1,xeima2,xeima3,  yeima1,yeima2,yeima3
!
      xgei= xeima1*xmag + yeima1*ymag + gm1*zmag
      ygei= xeima2*xmag + yeima2*ymag + gm2*zmag
      zgei= xeima3*xmag + yeima3*ymag + gm3*zmag
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tmaggeo(xmag,ymag,zmag,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_mag_to_geo: MAG -> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xmag,ymag,zmag   cartesian mag coordinates
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp12/ xeoma1,xeoma2,xeoma3,  yeoma1,yeoma2,yeoma3
!
      xgeo= xeoma1*xmag + yeoma1*ymag + gd1*zmag
      ygeo= xeoma2*xmag + yeoma2*ymag + gd2*zmag
      zgeo= xeoma3*xmag + yeoma3*ymag + gd3*zmag
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tmaggsm(xmag,ymag,zmag,xgsm,ygsm,zgsm)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_mag_to_gsm: MAG -> GSM  system
! *   Author : P. Robert, CRPE, 2002
!
! *   input  : xmag,ymag,zmag cartesian mag coordinates
! *   output : xgsm,ygsm,zgsm cartesian gsm coordinates
!
! ----------------------------------------------------------------------
!
!
      call tmagsma(xmag,ymag,zmag,xx,yy,zz)
      call tsmagsm(xx,yy,zz,xgsm,ygsm,zgsm)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tmagsma(xmag,ymag,zmag,xsma,ysma,zsma)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_mag_to_sma: MAG -> SM   system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xmag,ymag,zmag   cartesian mag coordinates
! *   output : xsma,ysma,zsma   cartesian sma coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xsma= cphi*xmag - sphi*ymag
      ysma= sphi*xmag + cphi*ymag
      zsma= zmag
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsmagei(xsma,ysma,zsma,xgei,ygei,zgei)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sma_to_gei: SM  -> GEI  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xsma,ysma,zsma   cartesian sma coordinates
! *   output : xgei,ygei,zgei   cartesian gei coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp01/ gs1,gs2,gs3, gm1,gm2,gm3, ge1,ge2,ge3,gr1,gr2,gr3
      common /timp08/ xeism1,xeism2,xeism3,  yeism1,yeism2,yeism3
!
      xgei=xeism1*xsma + yeism1*ysma + gm1*zsma
      ygei=xeism2*xsma + yeism2*ysma + gm2*zsma
      zgei=xeism3*xsma + yeism3*ysma + gm3*zsma
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsmageo(xsma,ysma,zsma,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sma_to_geo: SM  -> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xsma,ysma,zsma   cartesian sma coordinates
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp02/ ps1,ps2,ps3, gd1,gd2,gd3, pe1,pe2,pe3, pr1,pr2,pr3
      common /timp13/ xeosm1,xeosm2,xeosm3,  yeosm1,yeosm2,yeosm3
!
      xgeo= xeosm1*xsma + yeosm1*ysma + gd1*zsma
      ygeo= xeosm2*xsma + yeosm2*ysma + gd2*zsma
      zgeo= xeosm3*xsma + yeosm3*ysma + gd3*zsma
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsmagsm(xsma,ysma,zsma,xgsm,ygsm,zgsm)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sma_to_gsm: SM  -> GSM  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xsma,ysma,zsma   cartesian sma coordinates
! *   output : xgsm,ygsm,zgsm   cartesian gsm coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xgsm=  cmu*xsma + smu*zsma
      ygsm=  ysma
      zgsm= -smu*xsma + cmu*zsma
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsmamag(xsma,ysma,zsma,xmag,ymag,zmag)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sma_to_mag: SM  -> MAG  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: terms of transformation matrix are given in common
!
! *   input  : xsma,ysma,zsma   cartesian sma coordinates
! *   output : xmag,ymag,zmag   cartesian mag coordinates
!
! ----------------------------------------------------------------------
!
!
      common /timp11/ stetq,ctetq, sdze,cdze, smu,cmu, sphi,cphi
!
      xmag=  cphi*xsma + sphi*ysma
      ymag= -sphi*xsma + cphi*ysma
      zmag=  zsma
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsphcar(r,teta,phi,x,y,z)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sph_to_car: SPH -> CAR  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: none
!
! *   input  :   r,teta,phi  spherical coordinates (angles in radians)
! *   output :   x,y,z       cartesian coordinates
!
! ----------------------------------------------------------------------
!
!
      sq=r*sin(teta)
      x=sq*cos(phi)
      y=sq*sin(phi)
      z=r*cos(teta)
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsr2gse(xsr2,ysr2,zsr2,rotx,roty,rotz,xgse,ygse,zgse)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sr2_to_gse: SR2 -> GSE  system
! *   Author : P. Robert, CETP, 2001
! *   Comment: local system, non time dependent
!
! *   input  : xsr2,ysr2,sr2z cartesian sr2 coordinates
!              rotx,roty,rotz cartesian gse coordinates of rotation axis
! *   output : xgse,ygse,zgse cartesian gse coordinates
!
! ----------------------------------------------------------------------
!
! *** set transform matrix with Spin axis terms
!
      rmod= sqrt(rotx**2 + roty**2 + rotz**2)
!
      rx=rotx/rmod
      ry=roty/rmod
      rz=rotz/rmod
!
      a= 1./sqrt(ry*ry + rz*rz)
!
      x1= (ry*ry + rz*rz)*a
      x2= -rx*ry*a
      x3= -rx*rz*a
!
      y1=  0.
      y2=  rz*a
      y3= -ry*a
!
      z1= rx
      z2= ry
      z3= rz
!
! *** tranform the input vector from sr2 to gse
!
      xgse= x1*xsr2 + y1*ysr2 + z1*zsr2
      ygse= x2*xsr2 + y2*ysr2 + z2*zsr2
      zgse= x3*xsr2 + y3*ysr2 + z3*zsr2
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsr2mfa(xsr2,ysr2,zsr2,bx,by,bz,rox,roy,roz,xm,ym,zm)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sr2_to_mfa: SR2 -> MFA  system
! *   Author : P. Robert, CETP, 2001
! *   Comment: local system, non time dependent
!
! *   input  : xsr2,ysr2,sr2z cartesian sr2 coordinates
!              bx,  by,  bz   cartesian sr2 coordinates of DC mag field
!              rox, roy, roz  cartesian gse coordinates of rotation axis
! *   output : xm,  ym,  zm   cartesian mfa coordinates
!
! ----------------------------------------------------------------------
!
!
! *** normalise R supposed undependant of time
!
      rmod= sqrt(rox**2 + roy**2 + roz**2)
!
      rx=rox/rmod
      ry=roy/rmod
      rz=roz/rmod
!
! *** tranform the vector from sr2 to mfa
!
      bperp= sqrt(bx*bx + by*by)
      b0= sqrt(bx*bx + by*by + bz*bz)
!
! **  first rotation
!
      sinphi=by/bperp
      cosphi=bx/bperp
!
      xp=  cosphi*xsr2 + sinphi*ysr2
      yp= -sinphi*xsr2 + cosphi*ysr2
      zp=  zsr2
!
! **  second rotation
!
      sintet= bperp/b0
      costet= bz/b0
!
      xmag= costet*xp - sintet*zp
      ymag= yp
      zmag= sintet*xp + costet*zp
!
! **  third rotation
!
      b= sqrt(ry*ry + rz*rz)
!
      sxm=  b*costet*cosphi - rx*sintet
      sym= -b*sinphi
!
      smperp= sqrt(sxm*sxm + sym*sym)
!
      sinpsi= sym/smperp
      cospsi= sxm/smperp
!
      xm=  cospsi*xmag + sinpsi*ymag
      ym= -sinpsi*xmag + cospsi*ymag
      zm=  zmag
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsr2sre(xsr2,ysr2,spifre,spipha,deltaT,xsre,ysre)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sr2_to_sre: SR2 -> SRef system
! *   Author : P. Robert, CRPE, 2001
! *   Comment: Z component is unchanged (spin axis)
!
! *   input  : xsr2, ysr2 cartesian sr2 coordinates
!              spifre     spin frequency in Hz
!              spipha     spin phase in radians, growing with time
!                         spipha= positive angle between the xsr axis
!                         and the component of the direction of the Sun
!                         in the xsr-ysr plane.
!              deltaT     (T -To) time (sec.), between the current time
!                         and the time where is measured the spin phase.
! *   output : xsre,ysre  cartesian sr coordinates
!
! ----------------------------------------------------------------------
!
!
      pi2= 2.*3.1415927
!
      depift= mod(pi2*spifre*deltaT,pi2)
      phicr= spipha -depift
      sinphi=sin(phicr)
      cosphi=cos(phicr)
!
      xsre=  cosphi*xsr2 -sinphi*ysr2
      ysre=  sinphi*xsr2 +cosphi*ysr2
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tsresr2(xsre,ysre,spifre,spipha,deltaT,xsr2,ysr2)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_sre_to_sr2: SRef-> SR2 system
! *   Author : P. Robert, CRPE, 2001
! *   Comment: Z component is unchanged (spin axis)
!
! *   input  : xsre, ysre cartesian sr coordinates
!              spifre     spin frequency in Hz
!              spipha     spin phase in radians, growing with time
!                         spipha= positive angle between the xsr axis
!                         and the component of the direction of the Sun
!                         in the xsr-ysr plane.
!              deltaT     (T -To) time (sec.), between the current time
!                         and the time where is measured the spin phase.
! *   output : xsr2,ysr2  cartesian sr coordinates
!
! ----------------------------------------------------------------------
!
!
      pi2= 2.*3.1415927
!
      depift= mod(pi2*spifre*deltaT,pi2)
      phicr= spipha -depift
      sinphi=sin(phicr)
      cosphi=cos(phicr)
!
      xsr2=  cosphi*xsre +sinphi*ysre
      ysr2= -sinphi*xsre +cosphi*ysre
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine tvdhgeo(xvdh,yvdh,zvdh,rlat,rlong,xgeo,ygeo,zgeo)
!
! ----------------------------------------------------------------------
!
! *   Class  : transform modules of Rocotlib Software
! *   Object : transforms_vdh_to_geo: VDH -> GEO  system
! *   Author : P. Robert, CRPE, 1992
! *   Comment: local system, non time dependent
!
! *   input  : xvdh,yvdh,zvdh   cartesian vdh coordinates
!              rlat,rlong       latitude and longitude of the point
!                               of observation (radians)
! *   output : xgeo,ygeo,zgeo   cartesian geo coordinates
!
! ----------------------------------------------------------------------
!
!
      q=cos(rlat)
      r1=q*cos(rlong)
      r2=q*sin(rlong)
      r3=sin(rlat)
!
      v1=r1
      v2=r2
      v3=r3
!
      q12=sqrt(r1**2 + r2**2)
!
      d1= -r2/q12
      d2=  r1/q12
      d3=  0.
!
      h1= -r1*r3/q12
      h2= -r2*r3/q12
      h3=  q12
!
      xgeo= v1*xvdh + d1*yvdh + h1*zvdh
      ygeo= v2*xvdh + d2*yvdh + h2*zvdh
      zgeo= v3*xvdh + d3*yvdh + h3*zvdh
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
!
