#! /usr/bin/env python3 # def analemma ( ecc, lon, obliq ): #*****************************************************************************80 # ## analemma() computes the analemma. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2017 # # Author: # # Original C version by Brian Tung. # Python version by John Burkardt. # # Input: # # double ECC, the orbital eccentricity. # # double LON, the longitude of the perihelion in degrees. # # double OBLIQ, the obliquity in degrees. # import matplotlib.pyplot as plt import numpy as np days = 365.242 # # Internally, longitude and obliquity are in radians. # lon = lon * np.pi / 180.0 obliq = obliq * np.pi / 180.0 # # Compute the data using vector operations. # n = 1001 t01 = np.zeros ( n ) eot = np.zeros ( n ) dec = np.zeros ( n ) f = np.linspace ( 0.0, 1.0, n ) tau = 2.0 * np.pi * f # # Set theta to the current longitude. # theta = np.arctan2 ( np.sqrt ( 1.0 - ecc * ecc ) * np.sin ( tau ), np.cos ( tau ) - ecc ) # # Rotate clockwise in XY plane by theta, corrected by lon. # x1 = np.cos ( theta - ( lon - np.pi / 2.0 ) ) y1 = np.sin ( theta - ( lon - np.pi / 2.0 ) ) z1 = 0.0 # # Rotate counter-clockwise in XZ plane by obliq. # x2 = np.cos ( obliq ) * x1 + np.sin ( obliq ) * z1 y2 = y1 z2 = - np.sin ( obliq ) * x1 + np.cos ( obliq ) * z1 # # Set t equal to real time from tau and # rotate counter-clockwise by t, corrected by lon # t = tau - ecc * np.sin ( tau ) x3 = np.cos ( t - ( lon - np.pi / 2.0 ) ) * x2 + np.sin ( t - ( lon - np.pi / 2.0 ) ) * y2 y3 = - np.sin ( t - ( lon - np.pi / 2.0 ) ) * x2 + np.cos ( t - ( lon - np.pi / 2.0 ) ) * y2 z3 = z2 t01 = t / 2.0 / np.pi eot = - np.arctan2 ( y3, x3 ) * 4.0 * 180 / np.pi * days / ( days + 1.0 ) dec = np.arcsin ( z3 ) * 180.0 / np.pi # # Plot the equation of time. # plt.plot ( t01, eot, 'b-' ) plt.grid ( True ) plt.xlabel ( '<---Normalized Date--->' ) plt.ylabel ( '<---Minutes Early/Late--->' ) plt.title ( 'Equation of Time' ) plt.savefig ( 'eot.png' ) plt.show ( block = False ) plt.close ( ) # # Plot the declination. # plt.plot ( t01, dec, 'b-' ) plt.grid ( True ) plt.xlabel ( '<---Normalized Date--->' ) plt.ylabel ( '<---Degrees North/South--->' ) plt.title ( 'Declination' ) plt.savefig ( 'declination.png' ) plt.show ( block = False ) plt.close ( ) # # Plot the analemma. # plt.plot ( eot, dec, 'b-' ) plt.grid ( True ) plt.xlabel ( '<---Minutes Early/Late--->' ) plt.ylabel ( '<---Degrees North/South--->' ) plt.title ( 'Analemma' ) plt.savefig ( 'analemma.png' ) plt.show ( block = False ) plt.close ( ) return def analemma_test ( ): #*****************************************************************************80 # ## analemma_test() tests analemma(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2017 # # Author: # # Original C version by Brian Tung. # Python version by John Burkardt. # import numpy as np import platform print ( '' ) print ( 'analemma_test():' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' analemma() computes the analemma curve for given values of' ) print ( ' eccentricity, longitude, and obliquity.' ) ecc = 0.01671 lon = 1.347 * 180.0 / np.pi obliq = 0.4091 * 180.0 / np.pi analemma ( ecc, lon, obliq ) # # Terminate. # print ( '' ) print ( 'analemma_test():' ) print ( ' Normal end of execution.' ) return def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 06 April 2013 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) analemma_test ( ) timestamp ( )