1 files changed, 139 insertions, 0 deletions
diff --git a/noao/imred/vtel/dephem.x b/noao/imred/vtel/dephem.x
new file mode 100644
index 00000000..6c8c315d
--- /dev/null
+++ b/noao/imred/vtel/dephem.x
@@ -0,0 +1,139 @@
+# EPHEM -- Calculate ephemeris data for the sun, return latitude and
+# longitude of sub-earth point.
+
+procedure ephem (month, day, year, hour, minute, second, image_r,
+	bn_degrees, cldc_degrees, verbose)
+ 
+int	month			# time of observation
+int	day			#
+int	year			#
+int	hour			#
+int	minute			#
+int	second			#
+real	image_r			# image radius
+real	bn_degrees		# solar latitude of sub-earth point (degrees)
+real	cldc_degrees		# Carrington longitude of disk center
+bool	verbose			# verbose flag
+
+double	radians_per_degree, pi, two_pi, st, d, dd
+double	ma, sin_ma, sin_two_ma, ml, e, e_squared, e_cubed
+double	ep, ea, r, image_r_squared, tl
+double	lan, bn, p, p_degrees
+double	sl1, sl2, cldc, cos_bn, x, cl1
+double	sin_three_ma, sec_bn, y
+double	dd_squared, dd_cubed, c, s, cl2, sln
+int	mac[12]
+
+data	mac/0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334/
+
+begin
+	# This version ignores lunar and planetary perturbations.
+	radians_per_degree = .017453292519943d+0
+	pi = 3.1415926536d+0
+	two_pi = pi + pi
+
+	d = double(365 * year + (year - 1)/4 + mac[month] + day)
+	if (month >= 3 && mod(year, 4) == 0)
+	    d = d + 1.d+0
+	st = double(second / 3600. + minute / 60. + hour)
+	d = d + st/24.d+0 -.5d+0
+	dd = d / 10000.d+0
+	dd_squared = dd * dd
+	dd_cubed = dd * dd * dd
+
+	# Mean anomaly.
+	ma = radians_per_degree * (358.475845d+0 + .985600267d+0 *
+	    d - 1.12d-5 * dd_squared - 7.d-8 * dd_cubed)
+	ma = mod(ma, two_pi)
+	sin_ma = sin(ma)
+	sin_two_ma = sin(2.d+0 * ma)
+	sin_three_ma = sin(3.d+0 * ma)
+
+	# Mean longitude.
+	ml = radians_per_degree *
+	    (279.696678d+0 + .9856473354d+0 * d + 2.267d-5 *  dd_squared)
+	ml = mod(ml, two_pi)
+
+	# Ecentricity.
+	e = 0.01675104d+0 - 1.1444d-5 * dd - 9.4d-9 * dd_squared
+	e_squared = e * e
+	e_cubed = e_squared * e
+
+	# Obliquity.
+	ep = radians_per_degree * (23.452294d+0 -
+	    3.5626d-3 * dd - 1.23d-7 * dd_squared + 1.03d-8 * dd_cubed)
+
+	# Eccentric anomaly.
+	ea = ma + (e - e_cubed/8.d+0) * sin_ma + e_squared * sin_two_ma/2.d+0 +
+	    3.d+0 * e_cubed * sin_three_ma/8.d+0
+
+	# Radius vector.
+	r = 1.00000003d+0 * (1.d+0 - e * cos(ea))
+
+	# Image radius.
+	image_r = real(961.18d+0 / r)
+	image_r_squared = double(image_r * image_r)
+
+	# True longitude.
+	tl = ml + (2.d+0 * e - e_cubed/4.d+0) * sin_ma + 5.d+0 * e_squared *
+	    sin_two_ma/4.d+0 + 13.d+0 * e_cubed * sin_three_ma/12.d+0
+	tl = mod(tl, two_pi)
+
+	# Longitude of ascending node of solar equator.
+	lan = radians_per_degree * (73.666667d+0 + 0.0139583d+0 *
+	    (year + 50.d+0))
+
+	# Solar latitude of sub-earth point.
+	bn = asin(sin(tl - lan) * .12620d+0)
+	bn_degrees = real(bn / radians_per_degree)
+	if (verbose) {
+	    call printf("B0 (degrees) = %10.5f\n")
+	        call pargr(bn_degrees)
+	}
+
+	# Position angle of rotation axis.
+	p = atan(-cos(tl) * tan(ep)) + atan(-cos(tl - lan) * .12722d+0)
+	p_degrees = p/radians_per_degree
+	if (verbose) {
+	    call printf("P-angle (degrees) = %10.5f\n")
+	        call pargr(real(p_degrees))
+	}
+
+	# Carrington longitude of disk center.
+	sl1 = (d + 16800.d+0) * 360.d+0/25.38d+0
+	sl2 = mod(sl1, 360.d+0)
+	sln = 360.d+0 - sl2
+	sln = radians_per_degree * sln
+
+	cos_bn = cos(bn)
+	sec_bn = 1.d+0/cos_bn
+	c = +1.d+0
+	s = +1.d+0
+	x = -sec_bn * cos(tl - lan)
+	if (x < 0.)
+	    c = -1.d+0
+	y = -sec_bn * sin(tl - lan) * .99200495d+0
+	if (y < 0.)
+	    s = -1.d+0
+
+	cl1 = tan(tl - lan) * 0.99200495d+0
+	cl2 = atan(cl1)
+	if (s == 1.d+0 && c == 1.d+0)
+	    cldc = sln + cl2
+	if (s == -1.d+0 && c == -1.d+0)
+	    cldc = sln + cl2 + pi
+	if (s == 1.d+0 && c == -1.d+0)
+	    cldc = sln + cl2 + pi
+	if (s == -1.d+0 && c == 1.d+0)
+	    cldc = sln + cl2
+	if (cldc < 0.d+0)
+	    cldc = cldc + two_pi
+	if (cldc > two_pi)
+	    cldc = mod(cldc, two_pi)
+
+	cldc_degrees = real(cldc / radians_per_degree)
+	if (verbose) {
+	    call printf ("L0 (degrees) = %10.5f\n")
+	        call pargr (cldc_degrees)
+	}
+end