!  
!   Programme de conversion des coordonnees digitales
!   d'un point de l'image (line/pixel) vers des coordonnees
!   geographiques (rlat/rlon) 
!
!   ( Convert pixel position from digital co-ordinates (line/pixel)
!   to geographical co-ordinates (rlat/rlon) )
!
!   ATTENTION :
!   La variable step depend de la resolution de l'image
!
!  Input :
!     line  : numero de ligne mesure du Sud au Nord
!     pixel : numero du pixel, meusre de Est vers West
!
!  Output :
!     Coordonnees geographiques (lat/lon): rlat rlon
!  ou Point hors planete
!  avec
!     rlat : latitude du pixel (+90 a -90 du Nord vers le Sud)
!     rlon : longitude du pixel (+180 a -180 de Est vers West)
!
!  Version EUMETSAT 
!  Ref: EUM TD 06 "Archive User Handbook" rev 2.6 November 2001
!  http://www.eumetsat.de/
!----------------------------------------------------------------

      Program georef

      implicit none
!     step : depends on resolution image
!            IR or WV image : 2500*2500 pixels
!            VSI image : 5000*5000 pixels
!     step : pas du radiometre vu par le satellite, en degree. L'image
!     represente un champ de 18°x18°. Ici image IR 2500*2500 
      real, parameter :: step=18.0/2500

      real(kind=4) :: rlat,rlon,line,pixel,cenlat,h,rflon
      real(kind=4) :: altitude,req,rpol,oblate,pi,deg_to_rad,rad_to_deg
      real(kind=4) :: tanx,tany,val1,val2,val,vmu
      real(kind=4) :: x,y,xt,yt,zt,teta,yk,cosrf,sinrf
      logical(kind=1) :: visible

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

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


!     Conversion de line/pixel en offset d'angle du point central
      x=-(pixel-1250.5)*step
      y=(line-1250.5)*step

      x=x*deg_to_rad
      y=y*deg_to_rad

!     Calcul des tangentes d'angle
      tanx=tan(x)
      tany=tan(y)

      val1 = 1. + (tanx**2)
      val2 = 1. + (tany**2)*((1.+oblate)**2)
      val=val1*val2
      yk=altitude/req
      yk=yk**2

      if(val.gt.(yk/(yk-1.))) then
       visible=.false.
       rlat=999.
       rlon=999.
       print*,"Point hors planete"
       goto 999
      else
       visible=.true.
      endif

      vmu=(altitude-(req*sqrt(yk-(yk-1)*val)))/val
      cosrf=cos(rflon)
      sinrf=sin(rflon)

      xt=(altitude*cosrf)+(vmu*(tanx*sinrf-cosrf))
      yt=(altitude*sinrf)+(vmu*(tanx*cosrf+sinrf))
      zt=vmu*tany/cos(x)

      teta=asin(zt/rpol)
      rlat=(atan(((tan(teta))*req)/rpol))*rad_to_deg 
      rlon=atan((yt/xt))*rad_to_deg


      print*,"Geographical co-ordinates (lat lon): ",rlat,rlon
 999  continue
      end