pro polcoords,glon,glat,x,y,perimlat=plat,position=pos,toplon=tlon,$
  todata=todata
;
; Given glat,glon on a polar stereographic projection, return equivalent 
;   cartesian x,y
; perimlat = perimeter latitude (default = 0.) (used only for polar stereog)
; position = [xll,yll,xur,yur] normalized lower left and upper right of plot
;            (default [0.,0.,1.,1.])
; toplon = longitude at top of projection
;
if keyword_set(todata) le 0 then begin
  if glon lt -180. or glon gt 180. then begin
    print,'>>> polcoords: bad glon=',glon
    return
  endif
  if glat lt -90. or glat gt 90. then begin
    print,'>>> polcoords: bad glat=',glat
    return
  endif
endif
if n_elements(tlon) eq 0 then tlon = 0.
if tlon lt -180. or tlon gt 180. then begin
  print,'>>> polcoords: bad toplon=',tlon
  return
endif
if n_elements(pos) eq 0 then pos = [0.,0.,1.,1.]
if n_elements(pos) ne 4 then begin
  print,'>>> polcoords: position must have 4 elements <<<'
  return
endif
if pos(0) lt 0. or pos(1) lt 0. or pos(2) lt 0. or pos(3) lt 0. or $
   pos(0) gt 1. or pos(1) gt 1. or pos(2) gt 1. or pos(3) gt 1. then begin
  print,'>>> polcoords: bad pos=',pos,' (must be in normalized coords)'
  return
endif
if pos(0) ge pos(2) or pos(1) ge pos(3) then begin
  print,'>>> polcoords: bad position=',pos,' (pos=xll,yll,xur,yur, and ',$
    'xll < xur and yll < yur)'
  return
endif
if n_elements(plat) eq 0 then plat = 0.
if plat le -90. or plat ge 90. then begin
  print,'>>> polcoords: bad perimeter latitude=',plat
  return 
endif
;
xll = pos(0) & yll = pos(1) & xur = pos(2) & yur = pos(3)
scale = 0.5*(xur-xll) / tan((90.-abs(plat))/2.*!dtor)
if keyword_set(todata) le 0 then begin
  r = scale*tan((90.-abs(glat))/2.*!dtor)
  rlon = (glon-tlon)*!dtor
  if plat ge 0. then x = .5*(xll+xur)-r*sin(rlon) else $
                     x = .5*(xll+xur)+r*sin(rlon)
  y = .5*(yll+yur)+r*cos(rlon)
endif else begin
  rtod = 1./!dtor
  xo = x-.5*(xll+xur) & yo = y-.5*(yll+yur)
  r = sqrt(xo^2+yo^2)
  if plat ge 0. then alph = atan(-xo,yo) else alph = atan(xo,yo)
  glon = alph*rtod+tlon
  if glon ge 180. then glon = glon - 360.
  if glon lt -180. then glon = glon + 360.
  glat = 90.-atan(r/scale)*2.*rtod
  if plat lt 0. then glat = -glat
; print,'polcoords: glon,glat=',glon,glat
endelse
return
end