PRO sun_crab_offset, time, RADEC2000=radec2000, WNOFFSET=wnoffset, ECLIPTIC_MODE = ecliptic_mode, SKYOFF=skyoff
;
; Informal procedure to calculate angular offset between the sun and Crab, and express result in terms of
; 		RHESSI offpointing command parameters.
;
; time = date,time at which offset is to be calculated (eg '2003-jun-15 15:00:00')
; RADEC2000 = optional 2-element vector = RA, DEC coordinates of the source (J2000 degrees).
;					Default = CRAB coordinates.
; WNOFFSET  = optional 2-element vector = target offset relative to source (degrees +ve, equatorial west, north)
;					Default = [0,0].
; ECLIPTIC_MODE = if set, output is expressed as ecliptic offsets.  New default is equatorial.
; SKYOFF		= output variable to hold target offset on sky relative to Sun (+W, +N, mag)
;
; 15-Jun-02 	Initial version (ghurford@ssl.berkeley.edu)
;  7-May-03 gh 	Add RADEC2000 keyword
;  8-May-03 gh 	Corrected labels for ecliptic offsets
;			   	Add /ECLIPTIC_MODE switch.
;				Change default from ecliptic to equatorial.
; 23-May-03 gh	Minor tweak to output format
;  9-Jun-03 gh	Minor code simplification
; 13-Jun-03 gh	Add WNOFFSET keyword and documentation.
; 16-Apr-04 gh	Minor change to output format.
; 24-Jun-04 gh	Add SKYOFF keyword
;
; Precess and convert Crab coordinates
;
jdate		= 1979 + ANYTIM(time, /UTIME) / (86400*365.25)	; Julian year
racrab = 5.57556*15			; = 05h34m32s (Chandra web site, J2000)
decrab = 22.01444			; = 22d00m52s
IF N_ELEMENTS(radec2000) NE 0 THEN BEGIN
	racrab = radec2000[0]
	decrab = radec2000[1]
ENDIF
; racrab = 270.294		; Mean of G5-1 and GRS 1758, J2000, per DS 4-Nov-02
; decrab = -25.411
fssunit = 0.9674		; = degree equivalent of 1 'FSS degree'
;
; Intermediate outputs
;
PRINT
PRINT, time, 								FORMAT="('TIME:', 35X, A)"
PRINT, fssunit,								FORMAT="('FSSunit (deg):',      22X,  F10.4)"
PRINT,   racrab, decrab,					FORMAT="('CRAB(input):',        25X, 2F10.4)"
PRECESS, racrab, decrab, 2000.0, jdate
PRINT,   racrab, decrab,					FORMAT="('CRAB(precessed):',    21X, 2F10.4)"
IF N_ELEMENTS( WNOFFSET) GT 0 THEN BEGIN
	PRINT, wnoffset,						FORMAT="('OFFSET(eq.deg.+W+N onsky)',12X, 2F10.4)"
	racrab = racrab - wnoffset[0]/COS(decrab*!DTOR)
	decrab = decrab + wnoffset[1]
	PRINT, racrab, decrab,					FORMAT="('TARGET(ra,dec):', 	22X, 2F10.4)"
ENDIF
EULER, racrab, decrab, eclongcrab, eclatcrab, 3			; convert from ra,dec to ecliptic longitude,latitude
IF KEYWORD_SET(ecliptic_mode) THEN	$
	PRINT, eclongcrab, eclatcrab,			FORMAT="('TARGET ecl long,lat:',17X, 2F10.4)"
skyoff = FLTARR(3)
solradec 	= solephut(time)							; Output is consistent with FASAR,FASR epherides to 3 arcsec
PRINT, solradec[0], solradec[1],			FORMAT="('SOLAR RA,DEC:', 	    24X, 2F10.4)"
;
; Ecliptic-based calculations
;
IF KEYWORD_SET(ecliptic_mode) THEN BEGIN
	EULER, solradec[0], solradec[1], eclongsol, eclatsol, 3		; Convert solar coords to ecliptic longitude,latitude
	PRINT, eclongsol, eclatsol,				FORMAT="('SOLAR ECL long,lat:', 18X, 2F10.4)"
	avlat = (eclatsol + eclatcrab) / 2
	skyoff[0] = (eclongsol-eclongcrab)*COS(avlat*!DTOR)
	skyoff[1] =  eclatcrab - eclatsol
	skyoff[2] = SQRT(skyoff[0]^2 + skyoff[1]^2)
	PRINT, skyoff,							FORMAT="('TARGET rel to Sun(ecl.deg:+W,+N,mag):', 3F10.4)"
	PRINT
	xycmd = skyoff/fssunit
	PRINT, xycmd[1]*!DTOR, xycmd[0]*!DTOR, 	FORMAT="('RHESSI thetaX,thetaY (ecl.rad:+N,+W):', 2F10.3)" ; X=+ve N,Y=+veW
;
; Equatorial-based calculations
;
ENDIF ELSE BEGIN
	avlat = (solradec[1]+decrab)/2.
	skyoff[0] = (solradec[0] - racrab)*COS(avlat*!DTOR)
	skyoff[1] = decrab - solradec[1]
	skyoff[2] = SQRT(skyoff[0]^2 + skyoff[1]^2)
	PRINT, skyoff,								FORMAT="('TARGET rel to Sun(eq.deg:+W,+N,mag): ', 3F10.4)"
	PRINT
	xycmd = skyoff/fssunit
	PRINT, xycmd[1]*!DTOR, xycmd[0]*!DTOR, 		$
			FORMAT="('RHESSI thetaX,thetaY (eq.rad:+N,+W):', 2F10.3, 10X, '<=========')" ; X=+ve N,Y=+veW
ENDELSE
RETURN
END