FUNCTION geom_postel, phi0, lat0, nx, ny, mapscale

; Function which calculates the solar latitude and 
; longitude of a Postel projected space.
; INPUT:
; phi0 - the central longitude of the space.
; lat0 - the central latitude of the space.
; nx - the x-dimension of the space in pixels.
; ny - the y-dimension of the space in pixels.
; mapscale - the spatial scaling of the pixels in deg.
;
; This works for even nx, ny 
;
; OUTPUT:
; A a 3D datacube of 2 x 2D images. The first 2D image is the
; longitude and the second is the  latitude.
;
; EXAMPLE
; IDL> coord_map=geom_postel(133.28d0,-8.3d0,512,512,1.4)
;




phi0=phi0*!dpi/180.d0
lat0=lat0*!dpi/180.d0

dx=mapscale*!dpi/180.d0
x=(findgen(nx)-(nx-1)/2.d0)*dx
y=(findgen(ny)-(ny-1)/2.d0)*dx

x=x#replicate(1.d0,ny)
y=y##replicate(1.d0,nx)

c=sqrt(x^2+y^2)

lat = asin( cos(c)*sin(lat0) + y*sin(c)*cos(lat0)/c )  *180.d0/!dpi

phi = phi0*180.d0/!dpi + atan( x*sin(c) / (c*cos(lat0)*cos(c)-y*sin(lat0)*sin(c)) ) *180.d0/!dpi


;!p.multi=[0,2,1]
;plot, phi(*,100) - geom(*,100,0), ytitle='longitude difference'
;plot, lat(100,*) - geom(100,*,1), xtitle='latitude difference'

ofile=fltarr(512,512,2)
ofile(*,*,0)=phi
ofile(*,*,1)=lat

RETURN,ofile


end