!  
!   Programme de conversion des coordonnees geographiques
!   (rlat/rlon) vers des coordonnees digitales (line/pixel).
!
!   ( Convert pixel position from geographical co-ordinates (lat/long)
!   to digital co-ordinates (line/pixel) )
!
!
!   ATTENTION :
!   Les variables nlines et nsamps dependent de la resolution de l'image
!
!  Input :
!     rlat : latitude du pixel (+90 a -90 du Nord vers le Sud)
!     rlon : longitude du pixel (+180 a -180 de Est vers West)
!
!  Output :
!     Coordonnees du point en numeros de ligne/colonne : line pixel
!  avec
!     line  : numero de ligne mesure du Sud au Nord
!     pixel : numero du pixel, meusre de Est vers West
!     visible : flag True si pixel sur le disque
!                    False si pixel sur l'espace 
!
!  Version EUMETSAT
!  Ref: EUM TD 06 "Archive User Handbook" rev 2.6 November 2001
!  http://www.eumetsat.de/
!-----------------------------------------------------------------


      Program georef

      implicit none
! nlines/nsamps = 2500 for IR or WV image
! nlines/nsamps = 5000 for VIS image
      integer, parameter :: nlines=2500,nsamps=2500

      real(kind=4) :: rlat,rlon,lat,lon,geolat
      real(kind=4) :: altitude,req,rpol,oblate,pi,deg_to_rad,rad_to_deg
      real(kind=4) :: x,y,z,rtheta,aline,asamp,dotprod
      integer(kind=4) :: line,pixel
      logical(kind=1) :: visible

      print*,"lat long co-ordinates :"
      read*,rlat,rlon

!     Initialisation
!     altitude = distance du centre de la terre au satellite
!     req = rayon de la terre a l'equateur
!     rplo = rayon de la terre au pole
!     oblate = applatissement dela terre
!     deg_to_rad et rad_to_deg : facteurs de conversion
      
      altitude=42164.0
      req=6378.140
      rpol=6356.755
      oblate=1.0/298.257
      pi=acos(-1.)
      deg_to_rad=pi/180.0
      rad_to_deg=180.0/pi

!     Conversion en radians
      geolat=rlat*deg_to_rad
      lon=rlon*deg_to_rad

!     Conversion des latitudes geodetiques en geocentriques
      lat=atan(((1.0-oblate)**2)*tan(geolat))

!     Calcul de rtheta : distance du centre de la terre au point
!     de surface de latitude 'lat'
      rtheta=(req*rpol)/sqrt(rpol**2*cos(lat)**2+req**2*sin(lat)**2)

!     Calcul des coordonnees cartesiennes. Le systeme de coordonnes est
!     geocentrique avec l'axe des X vers le satellite, l'axe des Y vers
!     l'Est et l'axe des X vers le pole Nord
      x=rtheta*cos(lat)*cos(lon)
      y=rtheta*cos(lat)*sin(lon)
      z=rtheta*sin(lat)

!     Determination si le point est sur la terre ou dans l'espace.
!     Principe : le produit de 2 vectuers A et B est egal a :
!     |A||B| cos(theta)    theta etant l'angle forme par A et B.
!     L'horizon est defini comme le 'locus' de points pour lesquels la
!     normale est perpendiculaire a la ligne de visee du satellite. 
!     Tous les points visibles ont une valeur de theta inferieure a 90
!     et inversement.
!     Le test se resume a regarder le signe du produit de vecteur : si
!     le produit est negatif le point est invisible.
      
!     Les composantes du vecteur du point vers le satellite sont
!     Rs-x,-y,-z avec Rs la distance de l'origine au satellite. Les
!     composantes du vecteur normal sont x y z(Re/Rp)^2

      dotprod=(altitude-x)*x-y*y-z*z*((req/rpol)**2)
      if(dotprod.le.0) then
       visible=.false.
       line=0
       pixel=0
      else
       visible=.true.
!     Dans ce systeme, le satellite (S) a pour coordonnees (altitude,0,0)
!     le centre de la terre (O) (0,0,0) et le point point (P) (x,y,z).
!     Soit O'=(x,y,0) la projection de P a lavertical du plan equatorial
!     et O""=(x,0,z) la projection de P a la vertical du plan de Greenwich.
!     Les angles SO'P et SO"P permettent de determiner les coordonnes
!     angulaires (aline,asamp) de P dans le FOV
       asamp=atan(y/(altitude-x))
       aline=atan(z/sqrt(y**2+(altitude-x)**2))

!      conversion en degres
       asamp=asamp*rad_to_deg
       aline=aline*rad_to_deg

!      calcul de line/pixel. Les pixels sont mesures a partir de la
!      droite de l'image et la conversion angulaire etait mesuree dans la
!      direction x (est),il faut inclure une correction de signe pour les
!      pixels. L'image represente un champ 18°x18° divise sur une base
!      equi-angulaire
       asamp=asamp/(18.0/float(nsamps))
       aline=aline/(18.0/float(nlines))

       if(asamp.ge.0.0)then
        pixel=nsamps/2-int(asamp)
       else
        pixel=nsamps/2+1-int(asamp)
       endif
       if(aline.ge.0.0)then
        line=nlines/2+1+int(aline)
       else
        line=nlines/2+int(aline)
       endif

      endif

      print*,"Digital co-ordinates (line pixel): ",line,pixel
      stop
      end