;
;*  Copyright (c) 1995-2004 University Corporation for Atmospheric Research
;* All rights reserved
;
; /**----------------------------------------------------------------------    
; @file       xyz2g.f
;
; Convert XYZ position coordinates to lat/lon/height
;
; @author     S. V. Sokolovskiy, S. Syndergaard
; @version    $URL: svn://lynx.cosmic.ucar.edu/trunk/src/gmrion/xyz2g.f $ $Id: xyz2g.f 13339 2010-11-03 17:34:37Z jasonlin $
;
; -----*/
; /**----
; @sub        	xyz2g
; @parameter	r_ecef=[x,y,z] -- Input position
; @return       h              -- height above the Earth surface
; @             glat,glon      -- Latitude and longitude
; -----------------------------------------------------------------------*/

pro SUB_xyz2g,r_ecef,h,glat,glon

      remax=6378.1370d0
      remin=6356.7523d0

      a=remax^2
      b=remin^2
      c=remax*remin
      d=(remax/remin)^2
      e=180/acos(-1d0)

      ro2=r_ecef[*,0]^2+r_ecef[*,1]^2
      z2=r_ecef[*,2]^2
      r=sqrt(ro2+z2)
      ro=sqrt(ro2)

      glon=e*atan(r_ecef[*,1],r_ecef[*,0])
      glat=e*atan(d*r_ecef[*,2]/ro)
      h=r*(1.d0-c/sqrt(a*z2+b*ro2))

end



;
;*  Copyright (c) 1995-2004 University Corporation for Atmospheric Research
;* All rights reserved
;
; /**----------------------------------------------------------------------    
; @file       impact.f
;
; Compute impact parameters and corresponding heights, latitudes,
; longitudes, and index for identifying occultation flag (rising/settting).
;
; @author     S.V.Sokolovskiy
; @version    $URL: svn://lynx.cosmic.ucar.edu/trunk/src/gmrion/impact.f $ $Id: impact.f 13339 2010-11-03 17:34:37Z jasonlin $
;
; -----*/
; /**----
; @sub        	impact
; @parameter r1=[x1,y1,z1]  -- GPS position
; @		     r2=[x2,y2,z2]  -- LEO position
; @return	p         -- straight-line impact parameter (km)
; @         h         -- approximate height above Earth surface (km)
; @         glat,glon -- latitude and longitude of tangent point (deg)
; @         ind       -- integer index used for identifying occultation flag
; -----------------------------------------------------------------------*/
; /** @call xyz2g.f 
;*/

pro SUB_impact_parameter, r1, r2, rp;,glat=glat,glon=glon;,ind

      dim=size(r1,/dim) ; [time, x/y/z]

      rn=r2-r1
      ;xn=x2-x1
      ;yn=y2-y1
      ;zn=z2-z1
      rabs=sqrt(total(rn*rn, 2))
      rabs=transpose(congrid(transpose(rabs), dim[1],dim[0]))
      ;rabs=sqrt(xn^2+yn^2+zn^2)
      rn=rn/rabs
      ;xn=xn/rabs
      ;yn=yn/rabs
      ;zn=zn/rabs
      sp1=total(r1*rn,2)
      sp1=transpose(congrid(transpose(sp1), dim[1],dim[0]))
      ;sp1=x1*xn+y1*yn+z1*zn
      ;;sp2=x2*xn+y2*yn+z2*zn
      rp=r1-rn*sp1
      ;px=x1-xn*sp1
      ;py=y1-yn*sp1
      ;pz=z1-zn*sp1


      ;SUB_xyz2g,px,py,pz,h,glat,glon

      ;p=sqrt(px^2+py^2+pz^2) ; p = h + earth_radius(in tangent point)

      ;prd=sp1*sp2
      ;print, sp1 ; sp1 is negative array 
      ;print, sp2 ; sp2 is possitive array
      ;print, prd ; prd is negative array
      ;stop
      ;if (prd gt 0.) then ind=1
      ;if (prd eq 0.) then ind=0
      ;if (prd lt 0.) then ind=-1

end