function ephmrs,hr,min,doy,nyr

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	EPHMRS
; PURPOSE:
;	Calculate local ephemerides of the SUN.
; CATAGORY:
;	Coronal Solar Library
; CALLING SEQUENCE:
;	eph( hr, min, doy, nyr)
; INPUTS:
;	hr = Hour(UT)   ... integer
;	min = Minute(UT)   ... integer
;	doy = Day of the year   ... integer
;	nyr = Year   ... integer
; KEYWORD PARAMETERS:
; OUTPUTS:
;	Nine element structure named "angle" containing
;	  el = True anomaly + longitude of perigee
;	  rldst = 1.92236d3 / diam
;	  ra = Right Ascension
;	  dec = Declination
;	  dt = Mean sidereal time (?)
;	  p = P Angle
;	  b = B Angle
;	  diam = Solar Diameter in Arcseconds
;	  xlo = Carrington Longitude
; COMMON BLOCKS:
; NOTES:
;	Converted from ephmrs.f with certain "errors" corrected
;	References:
;	1)  Explanatory Supplement to the Astronomical Almanac, 1961
;	2)  Textbook on Spherical Astronomy,W.M.Smart,1977
;	3)  American Ephemeris & Nautical Almanac,1980,Current Year
; MODIFICATION HISTORY:
;	T. Henry, 20 Dec, 1991  ---  Created
;	          14 Mar, 1991  ---  Changes to XLO & DIAM
;-

if n_params(0) eq 0 or keyword_set(hlp) then begin
    print,'  '
    print,'  Calculate local ephemerides of the Sun.'
    print,'  '
    print,'  ephmrs, hr, min, doy, yr'
    print,'  '
    print,'       hr = Hour (UT)'
    print,'       min = Minute (UT)'
    print,'       doy = Day of the Year'
    print,'       yr = Year (YY)'
    print,'  '
    print,'  OUTPUTS:'
    print,'         Nine element structure named "angle" containing'
    print,'           el = True anomaly + longitude of perigee'
    print,'           rldst = 1.92236d3 / diam'
    print,'           ra = Right Ascension'
    print,'           dec = Declination'
    print,'           dt = Mean sidereal time (?)'
    print,'           p = P Angle'
    print,'           b = B Angle'
    print,'           diam = Solar Diameter in Arcseconds'
    print,'           xlo = Carrington Longitude'
    print,'  '
    return, {angle,el: 0.,rldst: 0.,ra: 0.,dec: 0.,dt: 0.,p: 0.,b: 0., diam: 0.,xlo: 0.}
endif

pi = 3.141592653589793d0
pi2 = 2.d0 * pi
pi5 = pi * 5.0d-1
rpd = pi / 1.8d2

nhr = fix(hr)
smin = size(min)
nsmin = n_elements(smin) - 2
if smin(nsmin) ne 5 then nmin = fix(min)
ndoy = fix(doy)
nnyr = fix(nyr)
dhr = double(nhr)
if smin(nsmin) eq 5 then dmin = min else dmin = double(nmin)
ddoy = double(ndoy)
dyr = double(nnyr)
if dyr gt 1.9d3 then begin
    dyr = dyr - 1.9d3
    yr = nnyr - 1900
endif else	$
    yr = nnyr

;Calculate the number of Leap days since the Epoch.
n = (yr - 1) / 4 - (yr-1) / 100 + (yr + 299) / 400

ut = (dhr + dmin/6.0d1) / 2.4d1
da = 3.65d2 * dyr + double(n) - 5.d-1 + ddoy + ut
delt = (da * 5.80566821d-8 - 2.75978511d-4) * da + 9.36861754d0
d = da + delt / 8.64d4
t = d / 3.6525d4
tsq = t * t
tcb = t * tsq

;   Calculate mean longitude of perigee,mean longitude of date
;Explanatory Supplement,pp 98
g1 = 2.812208444444444d2			;281d13'15.04"
g2 = 1.719175d0					;6189.03"
g3 = 4.527777777777778d-4			;1.63"
g4 = 3.333333333333333d-6			;0.012"
;g = (g1 + g2 * t + g3 * tsq + g4 * tcb) * rpd
g = (((g4 * t + g3) * t + g2) * t + g1) * rpd
;   Calculate mean anomaly
;Explanatory Supplement,pp 98
am1 = 3.584758333333334d2			;358d28'33.00"
am2 = 3.599904975000000d4			;129596579.1"
am3 = 1.5d-4					;0.54"
am4 = g4					;0.012"
;am = (am1 + am2 * t - am3 * tsq - am4 * tcb) * rpd
am = (am1 + t * (am2 - t * (am3 + t * am4))) * rpd
am = am mod pi2

;   Calculate eccentricity
;Explanatory Supplement,pp 98
e = 1.675104d-2 - t * (4.18d-5 + t * 1.26d-7)
esq = e * e
ecb = e * esq

;   Calculate true anomaly
;Smart,1977,pp 120 eqn 87   "Equation of the Center"
at = am + e * (sin(am) * (2.d0 - 2.5d-1 * esq) + $
	sin(2.d0 * am) * (1.25d0 * e) + sin(3.d0 * am) * (esq * 1.3d1 / 1.2d1))
;Smart,1977,pp121 eqn 89  Longitude = true anomaly + longitude of perihelion/perigee
el = (g+at)
el = el mod pi2
sel = sin(el)
cel = cos(el)
tel = tan(el)

;Smart,1977,pp 116 eqn 74i   WHAT IS THE VALUE OF "a"?
;setting a = 1.00000023 astronomical units reduces error by ~2.5d-4
rldst = 1.00000023d0 * (1.d0 - esq) / (1.d0 + e*cos(at))
;     diam = 1.92236d3 / rldst
;Semi-diameter at unit radius = 959'.63   Explanatory Supplement,pp489
;Semi-daimeter = 206264.8062ds / 2.<re>   Explanatory Supplement,1992, pp 398
;   where ds = diameter of the sun = 696000km                       eqn 7.2-4
;         re = earth-sun distance
diam = 1.91926d3 / rldst
;diam = 1.92231d3 / rldst
;diam = 1.91924d3 / rldst
;If I use diam=1.91931d3/rldst for 1988-1995, the difference between computed and
;AA book values are range from +.06 to -.09.
;If I use diam=1.92231d3/rldst, then I have 2 sine waves superimposed over the
;difference between computed values and the AA values for the year 1981!  The
;short period is eliminated with a 30 day offset (doy+30) given in the call to
;this program.  The long period is 384.61 days.
;longitude of ascending node of solar equator on the ecliptic   E. S. 1992,pp397
;omega = 75.76degrees + 0.01397degrees * T
;


;   Convert from Ecliptic to Equatorial coordinates
;   Mean obliquity of the Ecliptic
;Explanatory Supplement,pp 98
eom1 = 2.345229444444444d1		;23d27'8.26"
eom2 = 1.30125d-2	        		;46.845"
eom3 = 1.638888888888889d-6		;0.0059"
eom4 = 5.027777777777778d-7		;0.00181"
eom = (((eom4 * t - eom3) * t - eom2) * t + eom1) * rpd
seom = sin(eom)

;Explanatory Supplement, pp 26
;sin(delta) = sin(lamda)*sin(eps) iff |beta| << 1.0  => cos(beta)~1.0 & sin(beta) ~ 0.0
;maximum beta for 1991 is 0.97"  sin(beta) ~ 4e-6 cos(beta) ~ .9999999999
dec = asin(sel * seom)
p2 = cel / cos(dec)
ra = acos(p2)
i = fix(1.d0 + el / pi5)
case 1 of
(i eq 1): ra = ra
(i eq 3): ra = pi - ra
(i eq 4): ra = pi2 - ra
else: ra = pi + ra
endcase

;   Calculate mean sidereal time
;Astronomical Almanac 1980, pp 532
h1 = 6.646065555555556d0			;6h38m45.836s
h2 = 2.400051261666667d3			;8640185.542s
h3 = 2.580555555555556d-5			;0.0929s
;h = (h1 + h2 * t + h3 * tsq) * (3.6d2/2.4d1) * rpd
h = ((h3 * t + h2) * t + h1) * (3.6d2 / 2.4d1) * rpd
h = h mod pi2
dt = ra - h - pi
if dt lt -4. then dt = dt + pi2

;Longitude of ascending node of solar equator on ecliptic
;Explanatory Supplement, pp 307
;1854 Jan 1 at Greenwich mean noon = JD 2398220.0
;Sidereal Period of Rotation = 25.38 Mean Solar Days
;Elements: I = 7d15m & N = 73d40m for year 1850.0 from
;"Observations of Spots on the SUN" by R. C. Carrington pp 221,244
z1 = 7.366666666666667d1			;73d40'
z2 = 1.395833333333333d-2			;50.25"
;tyme = 50 years from 1850 to 1900 in which there were 12 leap years
;       minus 1 day to the epoch
daze = 1.8261d4
tyme = (daze + da)/3.6525d2
omega = (z1 + z2 * tyme) * rpd
elmz = el - omega
selmz = sin(elmz)
celmz = cos(elmz)
;Explanatory Supplement,pp 309-310
sinb = 1.261989691358298d-1*selmz		;sin(I)*sin(lam-omeg)
cosb = sqrt(1.d0 - sinb*sinb)
;cos(B) * cos(L - M) = -cos(lam - omega)
;cos(B) * sin(L - M) = -sin(lam - omega)* cos(I)
;=> tan(L - M) = tan(lam - omega) * cos(I)
clm = -celmz / cosb
slm = -selmz*9.920049496797150d-1 / cosb
xla = atan(slm,clm)
;xlb3 = 14.18439716d = 0.2475644328525329rad
;xlb2 = 2430000.5
;xlb1 = 112.766d = 1.968137964063901rad
;M - 180d = 112.766d + (2430000.5 -JD) * 14.18439716d
xlb = pi2 - (1.68d4 + d) * 2.475644328525329d-1
xlb = xlb mod pi2
;Carrington Longitude
;Rotation 1 commenced on 1853 Nov 9
;The synodic period of rotation is the interval of time during which
;L(lam) decreases by 360 degrees.  The mean synodic period is 27.2753 days.
;(L - 180d) = (L - M) + (M - 180d) ... + N * pi2 => positive L
xlo = xla + xlb + 5.d0 * pi2
xlo = xlo mod pi2
b = atan(sinb,cosb)
tanx = - cel * tan(eom)			;-cos(lam)tan(eps)
tany = - celmz * 1.272160680010469d-1	;-cos(lam-omeg)tan(I)
p = atan(tanx) + atan(tany)
a = {angle,el: el,rldst: rldst,ra: ra,dec: dec,dt: dt,p: p,b: b, diam: diam,xlo: xlo}
;openw,six,'eph.d',/get_lun,/append
;printf,six,nhr,nmin,ndoy,nnyr,da,d,delt,t,g,am,e,at,el,eom,h,dt,tyme,omega,elmz,xla,xlb,xlo
;free_lun,six
return,a
end