PRO sunpos,yr,mon,day,utsecs,gmst,eclong,sra,sdec
;
;     SUBROUTINE SUN (IYR, IDAY, SECS, GST, SLONG, SRASN, SDEC)
;     Program to calculate sidereal time and position of the Sun.
;     Good for years 1901 through 2099, accuracy 0.006 degree.
;     Input is IYR, IDAY (INTEGERs), and SECS, defining Un. Time.
;     Output is Greenwich Mean Sidereal Time (GST) in degrees,
;     longitude along ecliptic (SLONG), and apparent Right Ascension
;     and Declination (SRASN, SDEC) of the Sun, all in degrees.
;
;	Routine due to G. D. Mead, privately communicated to C. T. Russell.
;	Published in Russell, C. T., Geophysical Coordinate Transformations,
;	Cosmical Electrodynamics, vol. 2, pp. 184-196, 1971.
;
; Translated to IDL by D P Steele, 3 March 1994.
;
;     REAL*8 DJ, FDAY
;
;     DATA RAD /57.29578/
;
; Treat any vector arguments properly.
ny=N_ELEMENTS(yr)
nm=N_ELEMENTS(mon)
nd=N_ELEMENTS(day)
nt=N_ELEMENTS(utsecs)
nn=ny>nm>nd>nt		; number of elements in largest argument
IF ny LT nn THEN yy=REPLICATE(yr(0),nn) ELSE yy=yr
IF nm LT nn THEN mm=REPLICATE(mon(0),nn) ELSE mm=mon
IF nd LT nn THEN dd=REPLICATE(day(0),nn) ELSE dd=day
IF nt LT nn THEN tt=REPLICATE(utsecs(0),nn) ELSE tt=utsecs

;     FDAY = SECS / 86400.
fday=tt/86400.

; DJ is the number of days since January 0, 1900, 12:00 UT (aka December
; 31, 1899, 12:00 UT, which isn't nearly so neat-looking).
;     DJ = 365 * (IYR - 1900) + (IYR - 1901) / 4 + IDAY + FDAY - 0.5D0
dj0=julday(1,1,1900)-1L

; The following lines are borrowed from $IDL_IDLUSR/ymd2date.pro
; by Ray Sterner.
year=yy
w = WHERE(year LT 50, cnt)                 ; year < 50 assumed 20xx.
IF cnt GT 0 THEN year(w) = year(w) + 2000
w = WHERE(year LT 100, cnt)                ; 50 < year < 100 assumed 19xx.
IF cnt GT 0 THEN year(w) = year(w) + 1900

;     IF ((IYR .LT. 1901) .OR. (IYR .GT. 2099)) RETURN
early=WHERE(year LT 1901,nearly)
late=WHERE(year GT 2099,nlate)
IF (nearly+nlate GE 1) THEN BEGIN
    MESSAGE,'Results accurate only for 20th, 21st centuries!'
    RETURN
    ENDIF

IF nn EQ 1 $
    THEN dj=julday(mm,dd,year)+fday-0.5D0-dj0 $
    ELSE BEGIN
	dj=DBLARR(nn)
	FOR ii=0,nn-1 DO dj(ii)=DOUBLE(julday(mm(ii),dd(ii),year(ii))) $
			       +fday(ii)-0.5D0-dj0
    ENDELSE

; T is the number of Gregorian centuries since January 0, 1900, 12:00 UT
;     T = DJ / 36525.0D0
t=dj/36525.

;     VL = DMOD (279.696678 + 0.9856473354 * DJ, 360.0D0)
vl=(279.696678D0+0.9856473354D0*dj) MOD 360.0D0

;     GST = DMOD (279.690983 + 0.9856473354 * DJ + 360. * FDAY + 180.,
;    &            360.D0)
gmst=(279.690983D0+0.9856473354D0*dj+360.0D0*fday+180.0D0) MOD 360.0D0

;     G = DMOD (358.475845 + 0.985600267 * DJ, 360.D0) / RAD
g=((358.475845D0+0.985600267D0*dj) MOD 360.D0)*!DTOR

;     SLONG = VL + (1.91946 - 0.004789 * T) * SIN (G) + 0.020094 * 
;    &        SIN (2. * G)
eclong=vl+(1.91946D0-0.004789D0*t)*SIN(g)+0.020094D0*SIN(2.0D0*g)

;     OBLIQ = (23.45229 - 0.005686 * T) / RAD
obliq=(23.45229D0-0.005686D0*t)*!DTOR

;     SLP = (SLONG - 0.005686) / RAD
slp=(eclong-0.005686D0)*!DTOR

;     SIN_D = SIN (OBLIQ) * SIN (SLP)	! SIN_D, COS_D renamed because of
;     COS_D = SQRT (1. - SIN_D ** 2)	! SIND, COSD functions in FORTRAN 77
sin_d=SIN(obliq)*SIN(slp)
cos_d=SQRT(1.0D0-sin_d*sin_d)

;     SDEC = RAD * ATAN(SIN_D / COS_D)
sdec=!RADEG*ATAN(sin_d,cos_d)

;     SRASN = 180. - RAD * ATAN2 ((1.0 / TAN (OBLIQ)) * SIN_D / COS_D,
;    &                            - COS (SLP) / COS_D)
sra=180.0D0-!RADEG*ATAN((1.0D0/TAN(obliq))*(sin_d/cos_d),(-COS(slp)/cos_d))

;     RETURN
RETURN

;     END
END