<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_XYLL_PS'><A href='../../html_code/tools/da_xyll_ps.inc.html#DA_XYLL_PS' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_xyll_ps(x, y, proj, lat, lon) 2,2
! This is the inverse subroutine da_of llij_ps. It returns the
! latitude and longitude of an x/y point given the projection info
! structure.
implicit none
real, intent(in) :: x ! Column
real, intent(in) :: y ! Row
type (proj_info), intent(in) :: proj
real, intent(out) :: lat ! -90 -> 90 North
real, intent(out) :: lon ! -180 -> 180 East
real :: reflon
real :: scale_top
real :: xx,yy
real :: gi2, r2
real :: arccos
if (trace_use_frequent) call da_trace_entry
("da_xyll_ps")
! Compute the reference longitude by rotating 90 degrees to the east
! to find the longitude line parallel to the positive x-axis.
reflon = proj%stdlon + 90.0
! Compute numerator term of map scale factor
scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
! Compute radius to point of interest
xx = x - proj%polei
yy = (y - proj%polej) * proj%hemi
r2 = xx**2 + yy**2
! Now the magic code
if (r2 == 0.0) then
lat = proj%hemi * 90.0
lon = reflon
else
gi2 = (proj%rebydx * scale_top)**2.0
lat = deg_per_rad * proj%hemi * ASin((gi2-r2)/(gi2+r2))
arccos = ACOS(xx/sqrt(r2))
if (yy > 0) then
lon = reflon + deg_per_rad * arccos
else
lon = reflon - deg_per_rad * arccos
end if
end if
! Convert to a -180 -> 180 East convention
if (lon > 180.0) lon = lon - 360.0
if (lon < -180.0) lon = lon + 360.0
if (trace_use_frequent) call da_trace_exit("da_xyll_ps")
end subroutine da_xyll_ps