      program test_dir_sun
c
c     ***************************************************************0**
c *   Test_Dir_Sun: Sun direction in various coordinate systems
c     P. Robert, CETP, June 2001 ; use Rocotlib library                                  
c     ***************************************************************0**
c
c
      print*
      print 200, '*****************************************************'
      print*
      print 200, 'Sun direction in various coordinate systems'
      print*
c 
c *   reading date and time & compute general matrix
c
      call r_date(iyear,imonth,iday)
      call r_time(ih,im,is)
      call cp_time_param(iyear,imonth,iday,ih,im,is)
c
c *   reading position of point of observation in GEO
c
      print 200, 'latitude, longitude, distance in GEO ? (deg & any u.)'
      read *, rlageo, rlogeo,rdis
      print 200, ' ',rlageo, rlogeo,rdis
c
      rlagr=rlageo*3.1415927/180.
      rlogr=rlogeo*3.1415927/180.
      teta=(90.-rlageo)*3.1415927/180.
c
c *   geographic cartesian coordinates
c
      call t_sph_to_car(rdis,teta,rlogr,xgeo,ygeo,zgeo)
c
c *   same in other systems
c
      call t_geo_to_gei(xgeo,ygeo,zgeo,xgei,ygei,zgei)
      call t_geo_to_gse(xgeo,ygeo,zgeo,xgse,ygse,zgse)
      call t_geo_to_gsm(xgeo,ygeo,zgeo,xgsm,ygsm,zgsm)
      call t_geo_to_vdh(xgeo,ygeo,zgeo,rlagr,rlogr,xvdh,yvdh,zvdh)
c
c *   latitude, longitude in other systems
c
      call t_car_to_sph(xgeo,ygeo,zgeo,r,tetgeo,phigeo)
      call t_car_to_sph(xgei,ygei,zgei,r,tetgei,phigei)
      call t_car_to_sph(xgse,ygse,zgse,r,tetgse,phigse)
      call t_car_to_sph(xgsm,ygsm,zgsm,r,tetgsm,phigsm)
      call t_car_to_sph(xvdh,yvdh,zvdh,r,tetvdh,phivdh)
c
      call colalon(tetgeo,phigeo,rlageo,rlogeo)
      call colalon(tetgei,phigei,rlagei,rlogei)
      call colalon(tetgse,phigse,rlagse,rlogse)
      call colalon(tetgsm,phigsm,rlagsm,rlogsm)
      call colalon(tetvdh,phivdh,rlavdh,rlovdh)
c
      print*,'_________________________________________________________'
      print*
      print 300, 'date: YY/MM/JJ=', iyear,imonth,iday
      print 300, 'time: HH/MM/SS=', ih,im,is
      print*
c
      print 100, '            ','x','y','z','R','Lat','Long'
      print 100
c
      print 200, 'Pos. in GEO :',xgeo,ygeo,zgeo, r,rlageo,rlogeo
      print 200, 'Pos. in GEI :',xgei,ygei,zgei, r,rlagei,rlogei
      print 200, 'Pos. in GSE :',xgse,ygse,zgse, r,rlagse,rlogse
      print 200, 'Pos. in GSM :',xgsm,ygsm,zgsm, r,rlagsm,rlogsm
      print 200, 'Pos. in VDH :',xvdh,yvdh,zvdh, r,rlavdh,rlovdh
c 
c *   Sun direction
c
      call g_gei_geo_sun_dir(sxgei,sygei,szgei,sxgeo,sygeo,szgeo)
      call t_gei_to_gse(sxgei,sygei,szgei,sxgse,sygse,szgse)
      call t_gei_to_gsm(sxgei,sygei,szgei,sxgsm,sygsm,szgsm)
      call t_geo_to_vdh(sxgeo,sygeo,szgeo,rlagr,rlogr,sxvdh,syvdh,szvdh)
c
      call t_car_to_sph(sxgeo,sygeo,szgeo,rgeo,tetgeo,phigeo)
      call t_car_to_sph(sxgei,sygei,szgei,rgei,tetgei,phigei)
      call t_car_to_sph(sxgse,sygse,szgse,rgse,tetgse,phigse)
      call t_car_to_sph(sxgsm,sygsm,szgsm,rgsm,tetgsm,phigsm)
      call t_car_to_sph(sxvdh,syvdh,szvdh,rvdh,tetvdh,phivdh)
c
      call convdeg(tetgeo,phigeo)
      call convdeg(tetgei,phigei)
      call convdeg(tetgse,phigse)
      call convdeg(tetgsm,phigsm)
      call convdeg(tetvdh,phivdh)
c
      print*
      print 200, 'Sun  in GEO :',sxgeo,sygeo,szgeo,  rgeo,tetgeo,phigeo
      print 200, 'Sun  in GEI :',sxgei,sygei,szgei,  rgei,tetgei,phigei
      print 200, 'Sun  in GSE :',sxgse,sygse,szgse,  rgse,tetgse,phigse
      print 200, 'Sun  in GSM :',sxgsm,sygsm,szgsm,  rgsm,tetgsm,phigsm
      print 200, 'Sun  in VDH :',sxvdh,syvdh,szvdh,  rvdh,tetvdh,phivdh
      print*
      print*
      print*, '********************************************************'
      print*
c
  100 format(a,4a10  ,2a10  )
  200 format(a,4f10.3,2f10.3)
  300 format(a,3I5)
       
      stop 'test_dir_sun.exe: normal termination'
      end
c
c     ---------------------------------------------------------------0--
c
      subroutine convdeg(x,y)
c
      pisd=3.1415927/180.
c
      x=x/pisd
      y=y/pisd
c
      return
      end
c
c     ---------------------------------------------------------------0--
c
      subroutine colalon(teta,phi,rlat,rlon)
c
c *   compute lat and long in degres from teta, phi in radians
c
      pisd=3.1415927/180.
c
      rlat= 90. -teta/pisd
      rlon= phi/pisd
c
      return
      end
c
c     ---------------------------------------------------------------0--
c
