PRO spotpos, lat, lon, jd
;+
; NAME:
;       Spotpos
; PURPOSE:
;       Print statistics on the position of a spot on the solar disk
;       at a given date
; CALLING SEQUENCE:
;       Spotpos, Lat, Lon, JD
; INPUTS:
;       Lon:  Latitude of spot, North is positive
;       Lat:  Longitude of spot, west is positive
;       JD:   Julian day number
; OUTPUTS:
;       None, prints to screen
; RESTRICTIONS:
;       
; PROCEDURE:
;       
; MODIFICATION HISTORY:
;       25-Apr-2001  P.Suetterlin, SIU
;-

IF n_params() LT 3 THEN BEGIN
    message, 'Usage: Spotpos, Lat, Lon, JD', /inf
    message, '       North, West positive', /inf
    return
ENDIF

  ;;; compute angles P0 and B0
sun_pb, p0, b0, jd

  ;;; polar coordinates

tet = (90-lat)*!dtor
phi = (90-lon)*!dtor

  ;;; convert to cartesian

v = [ sin(tet)*cos(phi), sin(tet)*sin(phi), cos(tet)]

  ;;; rotate it by B0

v1 = rot3d(v, b0, /X)

  ;;; the required positions are the projection to x-z

r = sqrt(v1(0)^2+v1(2)^2)    ;;; distance from center
an = atan(-v1(0), -v1(2))    ;;; angle to center *without* p0
an1 = an - p0                ;;; angle to north direction

print
writeu, -1, 'Location at '+strtrim(abs(lat), 2)
IF lat LT 0 THEN writeu, -1, 'S 'ELSE writeu, -1, 'N '
writeu, -1, strtrim(abs(lon), 2)
IF lon LT 0 THEN writeu, -1, 'E' ELSE writeu, -1, 'W'
print
print, r, $
       form="('  Distance from center: ', f8.3)"
print, cos(asin(r)), $
       form="('  Cosine mu:            ', f8.3)"
print, an1*!radeg, $
       form="('  Angle Center<->North: ', f6.1)"

END