Compute and plot times scale differences

deltat.png

leaps.png

log_yearly.txt

WARNING: iauUtctai: dubious year for UTC (before 1960.0 or 5 years beyond last known leap second) [astropy.time.core]
WARNING: iauTaiutc: dubious year for UTC (before 1960.0 or 5 years beyond last known leap second) [astropy.time.core]
*** Yearly ***
1959 0.0 -0.944130239077
1960 -0.944130239077 -0.473363581114
1961 -1.41749382019 -0.428364076652
1962 -1.84585789684 -0.410249107517
1963 -2.25610700436 -0.51033468917
1964 -2.76644169353 -0.673364382237
1965 -3.43980607577 -0.873363926075
1966 -4.31317000184 -0.94672779087
1967 -5.25989779271 -0.946728419513
1968 -6.20662621222 -0.846727704629
1969 -7.05335391685 -0.94672779087
1970 -8.00008170772 -0.946728419513
1971 -8.94681012724 -1.05318969581
1972 -9.99999982305 -1.7499999376
1973 -11.7499997607 -1.25000013504
1974 -12.9999998957 -1.00000023376
1975 -14.0000001295 -0.999999605119
1976 -14.9999997346 -0.750000332482
1977 -15.7500000671 -1.25000013504
1978 -17.0000002021 -0.999999605119
1979 -17.9999998072 -1.00000023376
1981 -19.000000041 -1.00000023376
1982 -20.0000002747 -0.999999605119
1983 -20.9999998799 -1.00000023376
1985 -22.0000001136 -0.999999605119
1987 -22.9999997187 -1.00000023376
1989 -23.9999999525 -1.00000023376
1990 -25.0000001863 -0.999999605119
1992 -25.9999997914 -1.00000023376
1993 -27.0000000251 -1.00000023376
1994 -28.0000002589 -0.999999605119
1995 -28.999999864 -1.00000023376
1997 -30.0000000978 -0.999999605119
1998 -30.9999997029 -1.00000023376
2005 -31.9999999367 -1.00000023376
2008 -33.0000001704 -0.749999703839
2009 -33.7499998743 -0.24999990128
2012 -33.9999997756 -1.00000023376

plot_deltat.py

"""
Check out astropy.time.
https://github.com/astropy/astropy/issues/357
"""
import numpy as np
from astropy.time import Time

def time_scale_offset(times, format='jyear', scale1='utc', scale2='tai'):
	"""Compute offset between two time scales in seconds
	for a given list of times."""
	t1 = Time(times, format=format, scale=scale1)
	t2 = Time(times, format=format, scale=scale2)
	dt = t2 - t1
	return dt.sec
	
def plot_deltat():
	"""Reproduce time scale difference plot with astropy.time
	http://www.ucolick.org/~sla/leapsecs/deltatpng.html"""
	import matplotlib.pyplot as plt
	# Create an array of times at which we will compute the time scale difference
	times = np.arange(1950, 2020, 1./365)
	format = 'jyear'
	plt.clf()
	# Once delta_ut1_utc is set automatically we can add 'ut1':
	# https://github.com/astropy/astropy/issues/351
	for scale in ['tai', 'utc', 'tt', 'tcg', 'tdb', 'tcb']:
		diff = time_scale_offset(times, format, scale)
		plt.plot(times, diff, label='%s' % scale)
	plt.xlim(1950, 2020)
	plt.ylim(-60, 60)
	plt.title('Time scale differenct TAI - (OTHER SCALE)')
	plt.legend(loc='best')
	plt.savefig('deltat.png')


def check_leaps(times=np.arange(1950, 2020, 1)):
	"""TAI - UTC should always be an integer and only change when leap seconds occur
	http://en.wikipedia.org/wiki/Leap_second	
	
	This comment and the following code suggest that this might not be the case:
	https://github.com/astropy/astropy/issues/357#issuecomment-8068224
	"""
	import matplotlib.pyplot as plt
	diff = time_scale_offset(times, 'jyear', 'utc')
	jump = np.diff(diff)

	# print the numbers to illustrate jump
	for ii in np.nonzero(jump)[0]:
		print times[ii], diff[ii], jump[ii]

	# make a plot illustrating jump
	plt.clf()
	plt.plot(times[:-1], jump)
	plt.xlim(1950, 2020)
	plt.ylim(-1.2, 1.2)
	plt.title('np.diff(TAI - UTC) with one-day sampling')
	plt.savefig('leaps.png')

if __name__ == '__main__':
	plot_deltat()
	print('*** Daily ***')
	check_leaps(np.arange(1950, 2020, 1./365))
	print('*** Yearly ***')
	check_leaps(np.arange(1950, 2020, 1))