#==============================================================================
# SOLAR CALCULATIONS
# Version 3.0 of 23rd Feb 2014
# Kevin Karney - Kevin@Karney.com

# You need to put this file in the same directory as the calling file 
# and include the following three lines at the top of the calling file
#       import Astro_Subs 
#       reload(Astro_Subs) # Recompile AstroSubs, in case yoiu have made any changes
#       from Astro_Subs import*

# Just call Lets_See(0) which shows the complete range of output in number form
# Just call Lets_See(1) which shows the complete range of output in explanatory text & number form
#==============================================================================
# PLEASE FEEL FREE TO USE AND ABUSE THESE ROUTINES IN ANY FASHION
# JUST HAVE SOME FUN - email comments, complaints, ideas, bugs, enhancements to Kevin@Karney.com
from time import timezone, localtime
from math import degrees, radians, atan, atan2, tan, asin, sin, acos, cos, sqrt, pi
#==============================================================================
# YOU MUST CONFIGURE Using_Pythonista, on the next line, to True or False depending whether you are using the Pythonista program
Using_Pythonista = False # Does operating system have locational data as Apple IOS for iPad or iPhone
#==============================================================================
if Using_Pythonista : from location import get_location

# ===============================================================================================================================
# This Routine (astronomically totally boring & uninteresting) just sorts out the input, checks its type/range 
# and, if required and you are using Pythonista, accesses locational info from iPads and iPhones 
# and, depending on the first parameter,
# calls subroutines . . .
#  . . . . 'Solar_Position'       WHICH DOES THE MAIN ASTRONOMICAL CALCULATIONS. 
#                                 If you want to do a lot of calculations, you can call Solar_Position directly, not via this routine
#  . . . . 'Solar_Dawns_Dusks'    to do the iterative calculations needed to correctly compute sunrises/sets, twilight etc
#  . . . . 'Solar_Noon_Midnight'  to do the iterative calculations needed to correctly compute solar noon and midnight
# Sun('?') dumps Help
# ===============================================================================================================================
def Sun (*args):
	Output_Type = args[0].upper()
	if len(args) == 1 or args[0] == "?":
		Help()
	else: 
		if len(args) == 10 or len(args) == 12 or len(args) == 14:
			Output_Unit = args[1]
			if not(Output_Unit == 0 or Output_Unit == 1) : return "Wrong Output Unit, must be 0 or 1"
			Place       = str(args[2])
			
			Lat         = args[3]
			if Using_Pythonista  and  Lat == "*" : Lat = get_location()["latitude"]
			elif not(isinstance(Lat, int)    or isinstance(Lat, float))    or Lat<-90    or Lat > 90    : return "Latitude (entered as " + str(Lat)    + ") must be a number between -90 & 90"
			
			Long        = args[4]
			if Using_Pythonista  and  Long == "*" : Long = get_location()["longitude"]
			elif not(isinstance(Long, int)   or isinstance(Long, float))   or Long <-360 or Lat > 360   : return "Longitude (entered as " + str(Long)   + ") must be a number between -360 & 360"
			
			Zone        = args[5]
			if Zone == "*" : Zone = -timezone/60. # n.b. Python uses Timezone in seconds and +ve west
			elif not(isinstance(Zone, int)   or isinstance(Zone, float))   or Zone <-24  or Zone > 24   : return "Zone (entered as "     + str(Zone)   + ") must be a number between -24 & 24"
			
			DST         = args[6]
			if DST == "*" : DST = localtime()[8] 
			elif not(isinstance(DST, int)    or isinstance(DST, float))    or DST <0     or DST > 2     : return "DST (entered as "      + str(DST)    + ") must be a number between 0 & 2"
			
			Year        = args[7]
			if Year == "*" : Year = localtime()[0] 
			elif not(isinstance(Year, int)   or isinstance(Year, float))   or Year <200  or Year > 2050 : return "Year (entered as "     + str(Year)   + ") must be a number between 2000 & 2050"
			
			Month       = args[8]
			if Month == "*" : Month = localtime()[1] 
			elif not(isinstance(Month, int)  or isinstance(Month, float))  or Month <1   or Month > 12  : return "Month (entered as "    + str(Month)  + ") must be a number between 1 & 12"
			
			Day         = args[9]
			if Day == "*" : Day = localtime()[2] 
			elif not(isinstance(Day, int)    or isinstance(Day, float))    or Day <-1    or Day > 31    : return "Day (entered as "      + str(Day)   + ") must be a number between -1 & 31"
			
			Hour        = 12  #default value
			Minute      = 0   #default value
		
		if len(args) == 12 or len(args) == 14:
			Hour        = args[10]
			if Hour == "*" : Hour = localtime()[3] 
			elif not(isinstance(Hour, int)   or isinstance(Hour, float))   or Hour <-1   or Hour > 24   : return "Hour (entered as "     + str(Hour)   + ") must be a number between -1 & 24"
			
			Minute      = args[11]
			if Minute == "*" : Minute = localtime()[4] + localtime()[5]/60.
			elif not(isinstance(Minute, int) or isinstance(Minute, float)) or Minute <-1 or Minute  > 60: return "Minute (entered as "   + str(Minute) + ") must be a number between -1 & 60"
			
		Temperature     = 20.   #default value
		Pressure        = 1020. #default value
		
		if len(args) == 14 :
			Temperature = args[12]
			Pressure    = args[13]
		if not(len(args) == 10 or len(args) == 12 or len(args) == 14)  :
			return "Wrong Number of Arguments for " + Output_Type  + " in " + str(args) + ": must be 10,12 or 14"

		if  Output_Type[0:4]	 == "DAWN" or Output_Type[0:4] == "DUSK":
			Results = Solar_Dawns_Dusks(Output_Type,Place,Lat,Long,Zone,DST,Year,Month,Day)
		
		elif Output_Type[0:4]    == "DAYL":
			x1                   = Solar_Dawns_Dusks("DUSK" + Output_Type[4:] ,Place,Lat,Long,Zone,DST,Year,Month,Day)[0]
			x2                   = Solar_Dawns_Dusks("DAWN" + Output_Type[4:],Place,Lat,Long,Zone,DST,Year,Month,Day)[0]
			Results = x1 - x2
		
		elif Output_Type[0:4]    == "SOLA":
			if   Output_Type     == "SOLAP" : Results = Solar_Noon_Midnight(Place,Lat,Long,Zone,DST,Year,Month,Day,-1)
			elif Output_Type     == "SOLAN" : Results = Solar_Noon_Midnight(Place,Lat,Long,Zone,DST,Year,Month,Day, 0)
			elif Output_Type     == "SOLAF" : Results = Solar_Noon_Midnight(Place,Lat,Long,Zone,DST,Year,Month,Day, 1)
			elif Output_Type     == "SOLAL" : 
				x1               = Solar_Noon_Midnight(Place,Lat,Long,Zone,DST,Year,Month,Day, 0)[0]
				y1               = Solar_Noon_Midnight(Place,Lat,Long,Zone,DST,Year,Month,Day+1, 0)[0]
				Results = 24. + y1 - x1
			else : Results = "Did not understand the Type of Output required"
		else :
			Results = Solar_Position(Output_Type,Place,Lat,Long,Zone,DST,Year,Month,Day,Hour,Minute,Temperature,Pressure)
		return Output_Formatting(Output_Type, Output_Unit,Results)
			
# ===============================================================================================================================
# This Subroutune performs all the astronomical calculations concerning the Sun's position at a given LOCAL TIME at given location 
# It uses  Output_Type (e.g. RAD for Right Ascension in Degrees) to return what is actually required (see Help for full details)
#          Place (e.e. "Athens") is used just as a reference string
#          local Lat/Long/Time Zone/DST as place input : Long & Time Zone are +ve East of Greenwich: Lat is +ve North
#          local Year/Month/Day/Hour/Minute as time input
# It calls the sub-routines 'Get_Calendar_Day' & 'Get_Julian_Day' to get changing month/day if the UTC date differs from the local date
# ===============================================================================================================================
def Solar_Position(Output_Type,Place,Lat,Long,Zone,DST,Year,Month,Day,Hour,Minute,Temperature,Pressure) :
	if Output_Type == "LOC"     : return [Place,Lat,Long,Zone,DST]    # just use this to check your input
	if Output_Type == "CT"      : return [Year,Month,Day,Hour,Minute] # just use this to check your input
		
	#=========================================================================
	# CONVERT TO UTC TIME & UTC DATE
	# There is probably a smarter way to do this using Python's Time functions, but this does the job...
	Civil_Hrs			  = Hour + Minute/60.
	Standard_Hrs		  = Civil_Hrs - DST
	UTC_Hours			  = Standard_Hrs - Zone
	if UTC_Hours < 0 : # UTC is on the day before
		UTC_Hours         = UTC_Hours + 24.
		UTC_Date          = Get_Calendar_Day(Get_Julian_Day(Year,Month,Day, Hour,Minute) - 1)
		UTC_Year          = UTC_Date[0]
		UTC_Month         = UTC_Date[1]
		UTC_Day           = UTC_Date[2]
	elif UTC_Hours >= 0. and UTC_Hours < 24. : # UTC is today
		UTC_Year          = Year
		UTC_Month         = Month
		UTC_Day           = Day
	else : # UTC is on the day after
		UTC_Hours         = UTC_Hours - 24.
		UTC_Date          = Get_Calendar_Day(Get_Julian_Day(Year,Month,Day, Hour,Minute) + 1)
		UTC_Year          = UTC_Date[0]
		UTC_Month         = UTC_Date[1]
		UTC_Day           = UTC_Date[2]
	UTC_Deg				  = UTC_Hours * 15.
	if Output_Type == "UTC"    : return [UTC_Year, UTC_Month, UTC_Day,Hour,Minute]

	#=========================================================================
	# GET DAYS & JULIAN CENTURIES FROM EPOCH 2000 - which is 12 noon on 1st Jan 2000
	# I found this little routine on Wikipedia - it works during this century. No idea WHY it works. . . 
	Bbb	     			  = 367. * UTC_Year - 730531.5
	Ccc				      = int((7. * int(UTC_Year + (UTC_Month + 9.)/12.))/4.)
	Ddd					  = int(275. * UTC_Month / 9.) + UTC_Day
	J2000_Today		      = Bbb - Ccc + Ddd
	
	Days_2000_Now         = J2000_Today + UTC_Hours/24.
	T_Cent                = Days_2000_Now/36525.

	#=========================================================================
	# GET GREENWICH MEAN SIDEREAL TIME 
	#GMST_Hrs             = (6.697374558   + 0.06570982441908 * J2000_Today + 1.00273790935 * UTC_Hours + 0.000026 * T_Cent * T_Cent) % 24
	GMST_Hrs              = (6.69737483851 + 0.0657098242761  * J2000_Today + 1.00273790935 * UTC_Hours + 0.000026 * T_Cent * T_Cent) % 24
	GMST_Deg		      = 15. * GMST_Hrs	
	if Output_Type == "GMSTH"   : return GMST_Hrs
	if Output_Type == "GMSTD"   : return GMST_Deg
	if Output_Type == "MICA2"   : return  str(UTC_Year) + "-" + str(UTC_Month) + "-" + str(UTC_Day) + " " + str(Hour) + ":" + str(Minute) + "\t" + str(GMST_Hrs)		
	#=========================================================================
	# SET ECCENTRICITY, MEAN LONGITUDE AT PERIGEE & OBLIQUITY
	Mean_Long_Perigee_Deg = 248.54536 + 0.017196 * Year         # this formula adapted from Page C1 The Astronomical Almanac 2009
	Mean_Long_Perigee_Rad = radians(Mean_Long_Perigee_Deg)
	Eccentricity		  = 0.017585 - 0.438   * Year / 1000000.# this formula adapted from Page C1 The Astronomical Almanac 2009
	Obliquity_Deg         = 23.6993  - 0.00013  * Year          # this formula adapted from Page C1 The Astronomical Almanac 2009
	Obliquity_Rad		  = radians(Obliquity_Deg)
	if Output_Type == "M_LONP"  : return  Mean_Long_Perigee_Deg
	if Output_Type == "ECCENT"  : return  Eccentricity
	if Output_Type == "OBL"     : return  Obliquity_Deg

	#=========================================================================
	# GET SUN'S MEAN LONGITUDE
	Mean_Longitude_Deg    = (GMST_Deg - 180 - UTC_Deg) % 360
	Mean_Longitude_Rad    = radians(Mean_Longitude_Deg)
	if Output_Type == "M_LONS"  : return  Mean_Longitude_Deg

	#=========================================================================
	# SOLVE KEPLER'S EQUATION TO GET SUN'S TRUE LONGITUDE
	Mean_Anomaly_Rad	  = Mean_Longitude_Rad - Mean_Long_Perigee_Rad
	Eccentric_Anomaly_Rad = Mean_Anomaly_Rad - Eccentricity * sin(Mean_Anomaly_Rad) / (Eccentricity * cos(Mean_Anomaly_Rad) - 1)
	# Second Unneccesary Kepler iteration
	# Eccentric_Anomaly_Rad = Eccentric_Anomaly_Rad - (Mean_Anomaly_Rad - Eccentric_Anomaly_Rad + Eccentricity * sin(Eccentric_Anomaly_Rad)) / (-1 + Eccentricity * cos(Eccentric_Anomaly_Rad))

	True_Anomaly_Rad	  = 2 * atan(tan(Eccentric_Anomaly_Rad / 2) * sqrt((1 + Eccentricity) / (1 - Eccentricity)))
	True_Longitude_Rad	  = True_Anomaly_Rad + Mean_Long_Perigee_Rad
	True_Longitude_Deg    = degrees(True_Longitude_Rad)
	if Output_Type == "M_ANOM"  : return  Mean_Anomaly_Rad
	if Output_Type == "T_ANOM"  : return  True_Anomaly_Rad
	if Output_Type == "E_ANOM"  : return  Eccentric_Anomaly_Rad
	if Output_Type == "T_LONS"  : return  True_Longitude_Rad

	#=========================================================================
	# FIND SUN'S DISTANCE & ANGULAR SIZE
	Factor               = (1 + Eccentricity * cos(True_Anomaly_Rad)) / (1 - Eccentricity**2)
	Radius_Vector_AU     = 1 / Factor
	Angular_Size         = 0.533128 * Factor
	if Output_Type == "DIST"    : return  Radius_Vector_AU
	if Output_Type == "ANG_S"   : return  Angular_Size
	
	#=========================================================================
	# CONVERT TRUE LONGITUDE & OBLIQUITY TO RIGHT ASCENSION & DECLINATION
	Decl_Rad			  = asin(sin(Obliquity_Rad) * sin(True_Longitude_Rad))
	Decl_Deg			  = degrees(Decl_Rad)
			
	RA_Rad				  = atan2(cos(Obliquity_Rad) * sin(True_Longitude_Rad),cos(True_Longitude_Rad))
	RA_Deg				  = degrees(RA_Rad) % 360.
	RA_Hrs  			  = RA_Deg / 15.
	if Output_Type == "DECL"    : return Decl_Deg
	if Output_Type == "RAH"     : return RA_Hrs
	if Output_Type == "RAD"     : return RA_Deg
	if Output_Type == "RA_DEC"  : return [RA_Hrs,Decl_Deg]

	#=========================================================================
	# FIND LOCAL HOUR ANGLE
	LMST_Deg 			  = (GMST_Deg + Long) % 360.
	LMST_Hrs	          = LMST_Deg / 15.
	LHA_Deg	              = (LMST_Deg - RA_Deg) % 360
	if(LHA_Deg > 180): LHA_Deg = LHA_Deg - 360
	LHA_Rad	              = radians(LHA_Deg)
	if Output_Type == "LHA"     : return LHA_Deg
	if Output_Type == "MICA3"   : return str(UTC_Year) + "-" + str(UTC_Month) + "-" + str(UTC_Day) + " " + str(Hour) + ":" + str(Minute) + "\t" + str(LHA_Deg)
	
	#=========================================================================
	# FIND EQUATION OF TIME
	# EoT_Astro is the Equation of Time as formally defined by the IAU
	EoT_Astro_Deg_1       = GMST_Deg - RA_Deg - UTC_Hours * 15 + 180
	EoT_Astro_Deg         = EoT_Astro_Deg_1
	if EoT_Astro_Deg_1 < -180 :	 EoT_Astro_Deg = EoT_Astro_Deg + 360.
	if EoT_Astro_Deg_1 >  180 :	 EoT_Astro_Deg = EoT_Astro_Deg - 360.

	# EoT_Gnomon is the Equation of Time as normally used by gnomonists, 
	# i.e. the correction to be applied to solar time to get LOCAL mean time
	EoT_Gnomon_Deg        = -EoT_Astro_Deg
	EoT_Gnomon_Min        = EoT_Gnomon_Deg * 4
	if Output_Type == "EOT"     : return  EoT_Gnomon_Min
	if Output_Type == "EOTD"    : return  EoT_Gnomon_Deg
	if Output_Type == "MICA1"   : return  str(UTC_Year) + "-" + str(UTC_Month) + "-" + str(UTC_Day) + " " + str(Hour) + ":" + str(Minute) + "\t" + str(RA_Hrs) + "\t" + str(Decl_Deg) + "\t" + str(-EoT_Gnomon_Min)

	#=========================================================================
	# FIND EQUATION OF TIME LONGITUDE CORRECTED = SOLAR > LOCAL STANDARD TIME
	Long_Corr_Deg	      = Zone * 15. - Long
	# EoT_Sol_Std is the Equation of Time in its most useful form,
	# i.e. the correction to be applied to solar time to local STANDARD time (what you read on your watch)
	EoT_Sol_Std_Deg	      = EoT_Gnomon_Deg + Long_Corr_Deg
	EoT_Sol_Std_Min	      = EoT_Sol_Std_Deg * 4
	if Output_Type == "EOTL"    : return  EoT_Sol_Std_Min

	#=========================================================================
	# FIND APPROX SOLAR NOON
	Solar_Noon_Hrs = 12 + EoT_Sol_Std_Min/60.
	if Output_Type == "SOL_N"   : return  Solar_Noon_Hrs

	#=========================================================================
	# CONVERT RA,DECL & LHA TO ALTITUDE & AZIMUTH
	Lat_Rad				  = radians(Lat)
	Alt_Rad				  = asin(sin(Lat_Rad) * sin(Decl_Rad) + cos(Lat_Rad) * cos(Decl_Rad) * cos(LHA_Rad))
	Alt_Deg				  = degrees(Alt_Rad)
	Zenith_Deg            = 90. - Alt_Deg
	if Output_Type == "ALT"     : return  Alt_Deg
	if Output_Type == "ZEN"     : return  Zenith_Deg

	cos_Azim		      = ( sin(Decl_Rad) - sin(Alt_Rad) * sin(Lat_Rad) ) / ( cos(Alt_Rad) * cos(Lat_Rad))
	sin_Azim              = - cos(Decl_Rad) * sin(LHA_Rad)                  /   cos(Alt_Rad)
	Azim_Rad              = atan2(sin_Azim, cos_Azim)
	Azim_Deg			  = degrees(Azim_Rad) % 360
	if Output_Type == "AZ"      : return  Azim_Deg
	if Output_Type == "ALT_AZ"  : return  [Alt_Deg, Azim_Deg]
	if Output_Type == "MICA4"   : return  str(UTC_Year) + "-" + str(UTC_Month) + "-" + str(UTC_Day) + " " + str(Hour) + ":" + str(Minute) + "\t" + str(Alt_Deg) + "\t" + str(Azim_Deg)

	#=========================================================================
	# DO REFRACTION ESTIMATIONS
	# These are standard empiracle routines
	Refrac_Corr_Deg      = 0
	if (Alt_Deg > 15)  : Refrac_Corr_Deg = 0.00452 * tan(radians(Zenith_Deg)) * Pressure / (273 + Temperature)
	if (Alt_Deg <= 15) : Refrac_Corr_Deg = Pressure * (0.1594 + 0.0196 * Alt_Deg + 0.00002 * Alt_Deg * Alt_Deg) / ((273 + Temperature)  * (1 + 0.505 * Alt_Deg + 0.0845 * Alt_Deg * Alt_Deg) )
	Alt_Deg_Corr         = Alt_Deg - Refrac_Corr_Deg
	if Output_Type == "ALT_C"   : return  Alt_Deg_Corr
	if Output_Type == "REF_C"   : return  Refrac_Corr_Deg
	if Output_Type == "ALTR_AZ" : return  [Alt_Deg_Corr, Azim_Deg]

	#=========================================================================
	# FIND ESTIMATES OF SUNRISE & SUNSET - these (quite good) estimates assume that the Declination & EOT 
	# just calculayred are applicable throughout the day AND sunset/sunrise occur when the altitude of the
	# Sun's centre is zero without allowance for refraction
	Argument = -tan(Lat_Rad) * tan(Decl_Rad)
	if abs(Argument) > 1 : # Catch the uncalculable in Arctic Winter & Summer
		LHA_Sunset_rise      = 999.
		Sunrise_Hrs          = 999.
		Sunset_Hrs           = 999.
		Sunrise_Azimuth_Rad  = 999.
		Sunrise_Azimuth_Deg  = 999.
		Sunset_Azimuth_Deg   = 999.
	else :	 
		LHA_Sunset_rise      = degrees(acos(Argument))
		Sunrise_Hrs          = 12 - LHA_Sunset_rise/15 + EoT_Sol_Std_Min/60
		Sunset_Hrs           = 12 + LHA_Sunset_rise/15 + EoT_Sol_Std_Min/60
		Sunrise_Azimuth_Rad  = acos(sin(Decl_Rad)/cos(Lat_Rad))
		Sunrise_Azimuth_Deg  = degrees(Sunrise_Azimuth_Rad)
		Sunset_Azimuth_Deg   = 360 - Sunrise_Azimuth_Deg
	if Output_Type == "SUN_R"   : return  Sunrise_Hrs
	if Output_Type == "SR_AZ"   : return  Sunrise_Azimuth_Deg
	if Output_Type == "SUN_S"   : return  Sunset_Hrs
	if Output_Type == "SS_AZ"   : return  Sunset_Azimuth_Deg

	#=========================================================================
	# FIND ESTIMATES OF DAYLIGHT HOURS
	if Sunrise_Hrs <> 999. and Sunset_Hrs <> 999.:
		Daylight_Hrs = Sunset_Hrs - Sunrise_Hrs
	else :
		Daylight_Hrs = 999.
	if Output_Type == "ADAYL"   : return  Daylight_Hrs

	#=========================================================================
	# FIND DIRECTION COSINES FOR SUN'S POSITION
	Sun_cos_a				  = cos(Alt_Rad)  * cos(Azim_Rad)
	Sun_cos_b				  = sin(Alt_Rad)
	Sun_cos_c				  = cos(Alt_Rad)  * sin(Azim_Rad)
	if Output_Type == "COSIN"   : return  [Sun_cos_a, Sun_cos_b, Sun_cos_c]

	#=========================================================================
	# FIND X,Y,Z COORDINATES OF SUN
	Sun_Dist				 = 1000000.
	Sun_x					 = -Sun_Dist  * cos(Alt_Rad)  * sin(Azim_Rad)
	Sun_y					 =  Sun_Dist  * cos(Alt_Rad)  * cos(Azim_Rad)
	Sun_z					 =  Sun_Dist  * sin(Alt_Rad)
	if Output_Type == "XYZ"     : return [Sun_x, Sun_y, Sun_z]

	return "HELP - Got Output Type = " + Output_Type + ", which I do not understand"

# ===============================================================================================================================
# This Subroutune performs a Newton Raphson iteration to find the times, when . . .
# the solar azimuth reaches eother 0 or 180 for midday & midnight
# (it should be noted that 12 + EoT Longtitude corrected is NOT exactly solar noon. This routine gets over this problem)
# It calls the sub-routine 'Solar_Position' for thr values of Altitude & Azimuth needed in the iteration
# ===============================================================================================================================
def Solar_Noon_Midnight(Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day,Which):
	# Which = -1 Previous Midnight
	# Which =  0 Noon
	# Which =  1 Next Midnight
	Inc                 = .00000001
	Iteration_Precision = .000000001
	Max_Iterations      = 10
	The_Message         = "OK"
	if Which == 0 :      # Solar Noon requested
		The_Last_Hour   = 12.
		Calc_Day        = Set_Day
	else :                # Solar Midnight requested
		The_Last_Hour   = 0.
		if Which == -1 :  # Previous Solar Midnight requested
			Calc_Day    = Set_Day
		else :            # Next Solar Midnight (get previous midnight but on next day)
			Calc_Day    = Set_Day + 1
	The_Hour            = The_Last_Hour + Sun("EOTL",0,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Calc_Day,The_Last_Hour,0)/60.
	Azim                =                 Sun("Az"  ,0,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Calc_Day,The_Hour     ,0)

	if Azim > 270 :Azim = Azim - 360.
	if Azim  > 90 : # Check if Sun is North or South at Noon
		Target_Az       = 180.
	else :
		Target_Az       = 0.
	Closure             = 999.
	Azim_Last           = 999.
	counter             = 0
	
	while Closure > Iteration_Precision and counter < Max_Iterations: # Here's the Newton Raphson iteration
		# Check if the calculation time had crossed a day boundary
		if The_Hour + Inc > 0 : # Check if the calculation time had crossed a day boundary
			Day1        = Calc_Day
			Hour1       = The_Hour + Inc 
		else :
			Day1        = Calc_Day - 1
			Hour1       = 24. + The_Hour + Inc 
		if The_Hour - Inc > 0 :
			Day2        = Calc_Day
			Hour2       = The_Hour - Inc 
		else :
			Day2        = Calc_Day - 1
			Hour2       = 24. + The_Hour - Inc 
		zz1             = Sun("ALT_AZ",0,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month, Day1, Hour1,0,0,0)
		zz2             = Sun("ALT_AZ",0,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month, Day2, Hour2,0,0,0)
		Alt1 = zz1[0]
		Azi1 = zz1[1]
		Alt2 = zz2[0]
		Azi2 = zz2[1]
		if Azi1 > 270. : Azi1 = Azi1 - 360.
		if Azi2 > 270. : Azi2 = Azi2 - 360.
		Azim_Slope      = (Azi1 - Azi2) / (2 * Inc)
		Azim            = (Azi1 + Azi2) / 2.
		Alt             = (Alt1 + Alt2) / 2.
		The_Last_Hour   = The_Hour
		The_Hour        = The_Last_Hour - (Azim - Target_Az) / Azim_Slope
		Closure         = abs(Azim - Azim_Last)
		Azim_Last  	    = Azim
		counter         = counter + 1
	if Alt < 0 and Which == 0    : The_Message  = "Sun below Horizon at Solar Noon (Artic Winter)"
	if Alt > 0 and Which <> 0    : The_Message  = "Sun above Horizon at Solar Midnight (Artic Summer)"
	if counter == Max_Iterations : return[999.,999.,999.,"No Iteration Closure"]
	return [The_Hour,Alt, Azim ,The_Message,counter ]

# ===============================================================================================================================
# This Subroutune performs a Newton Raphson iteration to find the times at which the solar altitude reaches a specific value
# the specific value in question depends on the type of calculation requested
# It calls the sub-routine 'Solar_Noon_Midnight' to check if the requested condition is reached (e.g. is there sunset in the arctic summer
# It calls the sub-routine 'Solar_Position' for the values of Altitude & Azimuth needed in the iteration
# ===============================================================================================================================
def Solar_Dawns_Dusks(Output_Type,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day) :
	Output_Type        = Output_Type.upper()
	Output_Index       = {'DAWN0':0,'DAWNS':1,'DAWNC':2,'DAWNN':3,'DAWNA':4,'DUSK0':5,'DUSKS':6,'DUSKC':7,'DUSKN':8,'DUSKA':9,'SOLAN':10,'SOLMP':11,'SOLMN':12}
	The_Alt            = [0.,-0.833333333,-6.,-12.,-18.,0.,-0.833333333,-6.,-12.,-18.]
	The_Name           = ["Sunrise 0 Alt", "Civil Sunrise" , "Civil Dawn", "Nautical Dawn", "Astronomical Dawn","Sunset 0 Alt", "Civil Sunset" , "Civil Dusk", "Nautical Dusk", "Astronomical Dusk"]
	# The_Alt gives to degrees below zero altitude for Sunrise and Sunset 
	#              i.e. 0         for DAWN0 or DUSK0  = zero altitude of centre of Sun's disk without refraction
	#                   0.833 deg for DAWNS or DUSKS  = Civil Sunrise (n.b. .833333 deg = 34' for standard refraction + 16' for Sun's semi diameter
	#                   6 deg     for DAWNC or DUSKC  = start or end of civil twilight
	#                  12 deg     for DAWNN or DUSKN  = start or end of Nautical twilight   
	#                  18 deg     For DAWNA or DUSKA  = start or end of Astronomical twilight 	#                  Th  twilight   values of  6,12 & 18 were chosen to match these results with that of MICA
	#                  The 'standard' values are 5,10 & 15 (NOT SURE WHY THERE IS A MISMATCH HERE)
	Inc                 = .00000001
	Iteration_Precision = .000000001
	Max_Iterations      = 10
	
	Noon_Time           = 12. + Solar_Position("EOTL",Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day,12       ,0,0,0)/60.
	Noon_Alt            =       Solar_Position("ALT" ,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day,Noon_Time,0,0,0)
	if Output_Type[0:4] == "DAWN" :
		Midnight_Time   =       Solar_Position("EOTL",Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day,0       ,0,0,0)/60.
		Midnight_Alt    =       Solar_Position("ALT" ,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day, Midnight_Time,0,0,0)
	else :
		Midnight_Time   = 24 +  Solar_Position("EOTL",Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day+1,0    ,0,0,0)/60.
		Midnight_Alt    =       Solar_Position("ALT" ,Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day+1, Midnight_Time,0,0,0)
	
	Target_Alt          = The_Alt[Output_Index[Output_Type]] # This is the altitude we are iterating to find
	
	if Noon_Alt < Target_Alt : # abandon, no result possible"
		return [999.,999.,999.,"Throughout this Day of the Year, Sun is too far below the horizon to give " + str(Target_Alt)]
	if Midnight_Alt > Target_Alt :
		return [999.,999.,999.,"Throughout this Day of the Year, Sun's Altitude is too high to give "       + str(Target_Alt)]
	
	if Output_Type[0:4] == "DAWN" :
		The_Hour            = 6.
	else :
		The_Hour            = 18.
	Closure             = 999.
	counter             = 0
	# Here's the Newton Raphson iteration looking for a specific solar altitude
	while Closure > Iteration_Precision and counter < Max_Iterations:
		zz1             = Solar_Position("ALT_AZ",Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day,The_Hour + Inc,0,0,0)
		zz2             = Solar_Position("ALT_AZ",Set_Place,Set_Lat,Set_Long,Set_Zone,Set_DST,Set_Year,Set_Month,Set_Day,The_Hour - Inc,0,0,0)
		#print zz1, zz1
		Alt1 = zz1[0]
		Azi1 = zz1[1]
		Alt2 = zz2[0]
		Azi2 = zz2[1]
		Alt_Slope       = (Alt1 - Alt2) / (2 * Inc)
		The_Alt         = (Alt1 + Alt2) / 2.
		The_Azim        = (Azi1 + Azi2) / 2.
		The_Hour_Last   = The_Hour
		The_Hour        = The_Hour_Last - (The_Alt - Target_Alt) / Alt_Slope
		Closure         = abs(The_Hour - The_Hour_Last)
		counter         = counter + 1
	if counter          == Max_Iterations : return[999.,999.,999.,"No Iteration Closure"]
	return [The_Hour, The_Alt, The_Azim,"OK"]

# ===============================================================================================================================
# This Subroutune just formats the output requested - either with an explanatory string
# or as a single number or as an array of values
# ===============================================================================================================================
def Output_Formatting(Output_Type, Output_Unit,The_Values) :
	The_Months            = {1:"Jan",2:"Feb",3:"Mar",4:"Apr",5:"May",6:"Jun",7:"Jul",8:"Aug",9:"Sep",10:"Oct",11:"Nov",12:"Dec"}
	Output_String = {'?'      : "?       Help                                                ",
	                 'GMSTH'  : "GMSTH   Greenwich Mean Sideral Time           -Hrs-         ",
                    'GMSTD'  : "GMSTD   Greenwich Mean Sideral Time           -Deg-         ",
                    'DECL'   : "DECL    Sun's Declination                     -Deg-         ",
                    'RAH'    : "RA      Sun's Right Ascension                 -Hrs-         ",
                    'RAD'    : "RA      Sun's Right Ascension                 -Deg-         ",
                    'EOT'    : "EOT     Equation of Time - Gnomonical         -Min-         ",
                    'EOTD'   : "EOTD    Equation of Time - Gnomonical         -Deg-         ",
                    'EOTL'   : "EOTL    Equation of Time - Long  Corr         -Min-         ",
                    'ALT'    : "ALT     Sun's Altitude                        -Deg-         ",
                    'AZ'     : "AZ      Sun's   Azimuth                       -Deg-         ",
                    'ALT_C'  : "ALT_C   Sun's Altitude - Refract Corr         -Deg-         ",
                    'ZEN'    : "ZEN     Zenith Distance = 90-Alt              -Deg-         ",
                    'REF_C'  : "REF_C   Refraction Correction                 -Deg-         ",
                    'LHA'    : "LHA     Sun's Local Hour Angle                -Deg-         ",
                    'ECCENT' : "ECCENT  Earth's Eccentricity                  -             ",
                    'OBL'    : "OBL     Earth's Obliquity                     -Deg-         ",
                    'M_LONP' : "M_LONP  Sun's Mean Longitude at Perigee       -Deg-         ",
                    'M_LONS' : "M_LONS  Sun's Mean Longitude                  -Deg-         ",
                    'M_ANOM' : "M_ANOM  Sun's Mean Anomaly                    -Deg-         ",
                    'T_ANOM' : "T_ANOM  Sun's True Anomaly                    -Deg-         ",
                    'E_ANOM' : "E_ANOM  Sun's Eccentric Anomaly               -Deg-         ",
                    'T_LONS' : "T_LONS  Sun's True Longitude                  -Deg-         ",
                    'DIST'   : "DIST    Sun's Geocentric Distance             -AU-          ",
                    'ANG_S'  : "ANG_S   Sun's Angular Size                    -Deg-         ",
                    'SUN_R'  : "SUN_R   Approx Time of Sunrise                -Hrs-         ",
                    'SR_AZ'  : "SR_AZ   Approx Sunrise Azimuth                -Deg-         ",
                    'SUN_S'  : "SUN_S   Approx Time of Sunset                 -Hrs-         ",
                    'SS_AZ'  : "SS_AZ   Approx Sunset Azimuth                 -Deg-         ",
                    'ADAYL'  : "ADAYL   Approx Length of Day                  -Hrs-         ",
                    'SOL_N'  : "SOL_N   Approx Solar Noon                     -Hrs-         ",
                    'COSIN'  : "COSIN   Direction Cosines of Sun              -             ",
                    'XYZ'    : "XYZ     x-y-z coodinates of Sun               -             ",
                    'RA_DEC' : "RA_DEC  RA -Hrs- Hours & Decl                 -Deg-         ",
                    'ALT_AZ' : "ALT_AZ  Altitude & Azimuth                    -Deg-         ",
                    'ALTR_AZ': "ALTR_AZ Refrac Corr. Alt & Azi                -Deg-         ",
                    'MICA1'  : "MICA1   Date, UTC Time,RA,Decl,EoT-y,m,d,h,m,Hrs,Deg,Mins-  ",
                    'MICA2'  : "MICA2   Date, UTC Time,GMST       -y,m,d,h,m,Hrs-           ",
                    'MICA3'  : "MICA3   Date, UTC Time,LHA        -y,m,d,h,m,Deg-           ",
                    'MICA4'  : "MICA4   Date, UTC Time,Alt,Az     -y,m,d,h,m,Deg,Deg-       ",
                    'LOC'    : "LOC     Location                              -String-      ",
                    'CT'     : "CT      Civil Time                            -String-      ",
                    'UTC'    : "UTC     UTC                                   -String-      ",
                    'SOLAN'  : "SOLAN   Solar Noon                Time,Alt,Az -Hrs,Deg,Deg- ",
                    'SOLAP'  : "SOLAP   Prev. Solar midnight      Time,Alt,Az -Hrs,Deg,Deg- ",
                    'SOLAF'  : "SOLAF   Foll. Solar midnight      Time,Alt,Az -Hrs,Deg,Deg- ",
                    'SOLAL'  : "SOLAL   Length Solar Day (solar noon to next) -Hrs-         ",
                    'DAWN0'  : "DAWN0   SunRise - 0 Alt, no Refr. Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DAWNS'  : "DAWNS   Civil SunRise             Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DAWNC'  : "DAWNC   Start of Civil Dawn       Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DAWNN'  : "DAWNN   Start of Nautical Dawn    Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DAWNA'  : "DAWNA   Start of Astron Dawn      Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DUSK0'  : "DUSK0   Sunset - 0 Alt, no Refr,  Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DUSKS'  : "DUSKS   Civil SunSet              Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DUSKC'  : "DUSKC   End of Civil Dusk         Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DUSKN'  : "DUSKN   End of Nautical Dusk      Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DUSKA'  : "DUSKA   End of Astronomical Dusk  Time,Alt,Az -Hrs,Deg,Deg- ",
                    'DAYL0'  : "DAYL0   Length of Daylight - 0 Alt            -Hrs-         ",
                    'DAYLS'  : "DAYLS   Length of Civil Daylight              -Hrs-         ",
                    'DAYLC'  : "DAYLC   Daylight + Civil Twilight             -Hrs-         ",
                    'DAYLN'  : "DAYLN   Daylight + Nautical Twilight          -Hrs-         ",
                    'DAYLA'  : "DAYLA   Daylight + Astronom. Twilight         -Hrs-         "
                    }
	if Output_Type == "LOC" :
		string = The_Values[0] + " : Lat. " + Dec_Deg_DMS(The_Values[1]) + " : Long. " + Dec_Deg_DMS(The_Values[2]) + " : Zone " + str(The_Values[3]) + " : DST " + str(The_Values[4])
		if Output_Unit == 0  : return string
		else                 : return Output_String[Output_Type] + " = " + string
	elif Output_Type == "CT" :
		string = '%00.0f' % The_Values[3] + ":" + '%0.0f' % The_Values[4] + " hrs on " + '%.0d' % The_Values[2] + "-" + The_Months[The_Values[1]] + "-" +'%.0f' % The_Values[0]
		if Output_Unit == 0  : return string
		else                 : return Output_String[Output_Type] + " = " + string
	elif Output_Type == "UTC" :
		string = '%00.0f' % The_Values[3] + ":" + '%0.0f' % The_Values[4] + " hrs on " + '%.0d' % The_Values[2] + "-" + The_Months[The_Values[1]] + "-" +'%.0f' % The_Values[0]
		if Output_Unit == 0  : return string
		else                 : return Output_String[Output_Type] + " = " + string
	else:
		if Output_Unit == 0  : return The_Values
		else                 : return Output_String[Output_Type] + " = " + str(The_Values) 
		
# ===============================================================================================================================
# This Subroutune just prints the Help - which you can get just by calling Sun("?")
# ===============================================================================================================================
def Help() :
	print "======================================================================================="
	print "YOU MUST CONFIGURE Using_Pythonista to True or False depending whether you are using Pythonista"
	print "This is on line 3 of these routines"
	print
	print "======================================================================================="
	print "You need to put the file with these routines in the same directory as the calling file"
	print "and include the following three lines at the top of the calling file"
	print "     import Astro_Subs"
	print "     reload(Astro_Subs)"
	print "     from Astro_Subs import*"
	print
	print "======================================================================================="
	print "These Routines are called thus. . . "
	print "     Sun(Output_Type, Output_Unit, . . . further parameters giving locational & time values)"
	print "e.g. Sun('DECL'     ,1           ,      'Llandogo',-2.,54.,0,0.,              2014,2,2,11,30)"
	print "     >>> Sun's Declination at 11:30 2nd Feb 2012 in Llandogo"
	print "or   Sun('DAYLS'    ,0           ,      'Llandogo',-2.,54.,0,0.,              2014,2,2)"
	print "     >>> Length of Time between between Civil Sunrise & Sunset on 2nd Feb 2012 in Llandogo"
	print "or   Sun('?')"
	print "     >>> Dumps this Help"
	print "See end of Help for expected errors"
	print 
	print "======================================================================================="
	print "PARAMETER 1 : OUTPUT TYPE : THE FOLLOWING KEYWORDS ARE VALID, upper or lowercase. . . ."
	print "======================================================================================="
	print "'?'        : Prints all of this . . . "
	print "'LOC'      : String with Location name, Latitude, Longitude, Time Zone & DST"
	print "'CT'       : String with Civil Time"
	print "'UTC'      : String with UTC Time & Date"
	print "'GMSTH'    : Greenwich Mean Sideral Time - Hrs"
	print "'GMSTD'    : Greenwich Mean Sideral Time - Deg"
	print "'DECL'     : Sun's Declination -Deg"
	print "'RAH'      : Sun's Right Ascension - Hrs"
	print "'RAD'      : Sun's Right Ascension - Deg"
	print "'EOT'      : Equation of Time - Gnomonical          i.e. local solar time > local mean time - Mins"
	print "'EOTL'     : Equation of Time - Longitude Corrected i.e. local solar time > standard   time - Mins"
	print "'ALT'      : Sun's Altitude - Deg"
	print "'AZ'       : Sun's Azimuth - Deg"
	print "'ZEN'      : Sun's Zenith Distance (90-Alt) - Deg"
	print "'ALT_C'    : Sun's Altitude - Refraction Corrected - Deg"
	print "'REF_C'    : Refraction Correction - Deg"
	print "'LHA'      : Sun's Local Hour Angle - Deg"
	print
	print "'ECCENT'   : Earth's Eccentricity -"
	print "'OBL'      : Earth's Obliquity - Deg"
	print "'M_LONP'   : Sun's Mean Longitude at Perigee - Deg"
	print "'M_LONS'   : Sun's Mean Longitude - Deg"
	print "'M_ANOM'   : Sun's Mean Anomaly - Deg"
	print "'T_ANOM'   : Sun's True Anomaly - Deg"
	print "'E_ANOM'   : Sun's Eccentric Anomaly - Deg"
	print "'T_LONS'   : Sun's True Longitude - Deg"
	print
	print "'SUN_R'    : Approximate Time of Sunrise - Hrs"
	print "'SR_AZ'    : Approximate Sun's Azimuth at Sunrise - Deg"
	print "'SUN_S'    : Approximate Time of SunSet - Hrs"
	print "'SS_AZ'    : Approximate Sun's Azimuth at Sunset - Deg"
	print "'ADAYL'    : Approximate Length of Day - Hrs"
	print "'SOL_N'    : Approximate Solar Noon - Hrs"
	print 
	print "'COSIN'    : Trio of Values giving Direction Cosines of Sun"
	print "'XYZ'      : Trio of Values giving x-y-z coodinates of Sun from  observer N=x-axis, E=y-axis, Vert = z-axis, Sun distance of 1 million units"
	print "             ==============================================================================="
	print "             The following return a pair of values . . ."
	print "'RA_DEC'   : Pair of Values - RA - Hrs & Decl - Deg" 
	print "'ALT_AZ'   : Pair of Values - Altitude & Azimuth - Deg"
	print "'ALTR_AZ'  : Pair of Values - Refraction Corrected Altitude & Azimuth - Deg"
	print
	print "             ==============================================================================="
	print "             The following return 5 values . . "
	print "             ----- the time     - Hrs (999. >>> not available e.g. sunset in high artic summer)"
	print "             ----- the altitude - Deg       (ditto)"
	print "             ----- the azimuth  - Deg       (ditto)"
	print "             ----- a message - usually 'OK',otherwise arctic winter and summer remarks"
	print "             ----- number of iterations required to find the time"
	print "'DAWN0'    : Time when Sun's Altitude is zero (not refraction corrected) in the morning"
	print "'DAWNS'    : Time of Civil Sunrise"
	print "'DAWNC'    : Time of start of Civil Twilight"
	print "'DAWNN'    : Time of start of Nautical Twilight"
	print "'DAWNA'    : Time of start of Astronomical Twilight"
	print 
	print "'DUSK0'    : Time when Sun's Altitude is zero in the evening"
	print "'DUSKS'    : Time of Civil Sunset"
	print "'DUSKC'    : Time of end of Civil Twilight"
	print "'DUSKN'    : Time of end of Nautical Twilight"
	print "'DUSKA'    : Time of end of Astronomical Twilight"
	print
	print "'DAYL0'    : Length of Time from zero altitude Sunrise & Sunset"
	print "'DAYLS'    : Length of Time between Civil Sunrise & Sunset"
	print "'DAYLC'    : Length of Time of Daylight & Civil Twilight"
	print "'DAYLN'    : Length of Time of Daylight & Nautical Twilight"
	print "'DAYLA'    : Length of Time of Daylight & Astronomical Twilight"
	print 
	print "'SOLAN'    : Pair of Values for Solar Noon - Time & Altitude"
	print "'SOLAP'    : Pair of Values for (P) Previous  Solar Midnight - Time & Altitude"
	print "'SOLAF'    : Pair of Values for (F) Following Solar Midnight - Time & Altitude"
	print "'SOLAD'    : Length of Solar Day - Solar Midnight to Midnight"
	print "             ==============================================================================="
	print "             The following return a formatted text string including tabs:"
	print "             --- these can be pasted in Excel for comparison with the MICA program"
	print "'MICA1'    : Date, UTC Time & RA, Decl & EoT"
	print "'MICA2'    : Date, UTC Time & Greenwich Mean Sidereal Time"
	print "'MICA3'    : Date, UTC Time & Local Hour Angle"
	print "'MICA4'    : Date, UTC Time & Altitude & Azimuth"
	print
	print "======================================================================================="
	print "PARAMETER 2 : OUTPUT UNITS : THE FOLLOWING NUMBERS ARE VALID. . . . "
	print "======================================================================================="
	print " 0         : Output single number or array of numbers (as appropriate)"
	print " 1         : String with explanatory text & unit + number (or array)"
	print "             If a result of 999 is returned, the parameter in uncalculable" 
	print "             e.g sunset in high arctic summer" 
	print
	print "======================================================================================="
	print "MAIN INPUT ARGUMENTS. . . . "
	print "           : You must enter 10, 12 or 14 main input arguments. . . . "
	print "           : Parameters 3 to 10 for all calculations"
	print "           : Parameters 11 to 12 for are not needed for dawn, dusk, solar noon. . ."
	print "           :       . . . solar midnight, length of solar day, length of daylight"
	print "           : Parameters 13 to 13 for are only needed for refraction calculations"
	print 
	print "           : If too few parameters are added, times default to local 12 midday"
	print "           : and 20 deg 1020 millibars"
	print "======================================================================================="
	print "Param 3    : String with Place Name"
	print "Param 4    : Latitude  - degs  : +ve North"
	print "Param 5    : Longitude - degs  : +ve East of Greenwich"
	print "Param 6    : Time Zone - hours : +ve East of Greenwich"
	print "Param 7    : DayLight Savings  : number of saved hours"
	print "Param 8    : Year"
	print "Param 9    : Month"
	print "Param 10   : Day"
	print "Param 11   : Hour    - NOT NEEDED FOR DAWN,DUSK,LENGTH OR DAY/DAYLIGHT CALCULATIONS"
	print "Param 12   : Minute  - NOT NEEDED FOR DAWN,DUSK,LENGTH OR DAY/DAYLIGHT CALCULATIONS"
	print "Param 13   : Temperature - deg Centigrade - ONLY NEEDED FOR REFRACTION CALCULATIONS"
	print "Param 14   : Pressure    - millibars      - ONLY NEEDED FOR REFRACTION CALCULATIONS"
	print "NOTE for Params 8 - 12"
	print "  for Apple IOS devices, using Pythonista, if you enter '*' - with the quotes, the current location is automatically entered"
	print "  for all devices,                         if you enter '*' for the date, time, DST & time zone, these are entered automatically"
	print "======================================================================================="
	print
	print "======================================================================================="
	print "YOU CAN EXPECT TO GET THE FOLLOWING ERRORS IF YOU USE THESE ROUTINES. . . . "
	print "======================================================================================="
	print "Errors expected . . .               Max   Min    Av   StdD"
	print "Greenw Mean Sid Time MicroSec      0.02 - 0.03  0.00  0.01" # Checked
	print "Right Ascension      Secs Time     2.4  - 3.1  -0.4   1.0"  # Checked
	print "Declination          Secs Arc     16.8  -17.7  -0.1   6.4"  # Checked
	print "Equation of Time     Secs Time     2.2  - 1.5   0.4   0.7"  # Checked
	print "Local Hour Angle     Secs Arc     38.1  -28.2   5.9  10.4"  # Checked
	print "Altitude             Secs Arc     19.9  -40.1  -7.7   7.3"  # Checked
	print "Azimuth              Secs Arc     58.4  -43.5   6.2  11.4"  # Checked
	print "Errors by comparison with MICA Version 2.2.2 at 6 hourly intervals from 2000 - 2050"
	print "======================================================================================="
	print

# ===============================================================================================================================
# This Subroutune just prints everything the set of routines will do
# Just call Lets_See(0) for numerical output or Lets_See(1) numbers + explanation
# ===============================================================================================================================
def Lets_See(Which) :
	The_Types = ['LOC','CT','UTC','GMSTH','GMSTD','DECL','RAH','RAD','EOT','EOTL','ALT','ALT_C','AZ','ZEN','LHA',
	             'ECCENT','OBL','M_LONP','M_LONS','M_ANOM','T_ANOM','E_ANOM','T_LONS','REF_C','DIST','ANG_S',
	             'SUN_R','SR_AZ','SUN_S','SS_AZ','ADAYL','SOL_N','COSIN','XYZ',
	             'RA_DEC','ALT_AZ','ALTR_AZ','MICA1','MICA2','MICA3','MICA4',
	             'SOLAN','SOLAP','SOLAF','SOLAL',
	             'DAWN0','DAWNS','DAWNC','DAWNN','DAWNA',
	             'DUSK0','DUSKS','DUSKC','DUSKN','DUSKA',
	             'DAYL0','DAYLS','DAYLC','DAYLN','DAYLA']
	Output_Unit = Which
	Place       = "Llandogo"
	Lat         =   51. + 44./60. +  6.4/3600.
	Long        = -( 2. + 41./60. + 42.4/3600.)
	Zone        = 0
	DST         = 0
	Year        = 2014
	Month       = 2
	Day         = 2
	Hour        = 11
	Minute      = 15
	Temperature = 20.
	Pressure    = 1010.
	#Help()
	counter = 0
	while counter <= len(The_Types) - 1:
		Output_Type      = The_Types[counter]
		if Using_Pythonista :
			print Sun(Output_Type,Output_Unit,Place,"*","*","*","*","*","*","*","*","*",Temperature,Pressure)
		else:
			print Sun(Output_Type,Output_Unit,Place,Lat,Long,Zone,DST,Year,Month,Day,Hour,Minute,Temperature,Pressure)
		counter = counter + 1

#===========================================
# ROUTINE TO GET JULIAN DAY FROM DATE & TIME
def Get_Julian_Day(Year, Month, Day, Hour, Minute) :
#Reference: Astronomical Algorithms 2nd Edition 1998 by Jean Meeus - Page 60-61	print "herer"
	if Month <= 2 :
		YYear					= Year - 1
		MMonth					= Month + 12
	else :
		YYear					= Year
		MMonth					= Month
	a							= int(YYear / 100)
	if YYear > 1582 :
		Switcher				= 1
	else :
		if YYear < 1582 :
			Switcher			    = 0
		else :
			if MMonth > 10 :
				Switcher			= 1
			else :
				if MMonth < 10 :
					Switcher		= 0
				else :
					if Day >= 15 :
						Switcher    = 1
					else :
						Switcher    = 0
	if Switcher == 0 :
		b						    = 0
	else :
		b						    = 2 - a + int(a / 4)

	c							    = int(365.25 * YYear) ;
	d							     = int(30.6001 * (MMonth + 1)) ;
	return b + c + d + Day + 1720994.5 + (Hour + Minute/60.)/24.

#======================================================
# ROUTINE TO GET CALENDAR DATE AND TIME FROM JULIAN DAY
def Get_Calendar_Day(This_JD) :
	#Reference: Practical Astronomy with your Calculator 3rd Edn : Duffet Smith - Page 8
	Month_List = ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"]
	JDD                          = This_JD + .5
	III                          = int(JDD)
	FFF                          = JDD - III
	if III > 2299160 :
		AA                        = int((III - 1867216.25) / 36524.25)
		BB                        = III + 1 + AA - int(AA / 4)
	else :
		BB                       = III
	CC                           = BB + 1524
	DD                           = int((CC - 122.1) / 365.25)
	EE                           = int(365.25 * DD)
	GG                           = int((CC - EE) / 30.6001)
	Day_inc_Frac                 = CC - EE + FFF - int(30.6001 * GG)
	Dayo                       = int(Day_inc_Frac)
	HH                           = 24 * (Day_inc_Frac - Dayo)
	Houro                      = int(HH)
	Minuteo                    = 60 * (HH - Houro)
	if GG < 13.5 :
		Montho                  = GG - 1
	else :
		Montho                  = GG - 13
	if Montho > 2.5 :
		Yearo                   = DD - 4716
	else :
		Yearo                   = DD - 4715

	return [Yearo, Montho, Dayo, Houro, Minuteo, Month_List[Montho-1]]