Здравствуйте, товарищи! Хочу с вами поделиться Python-скриптом для расчета времени восходазахода солнца.
Скрипт писал давно исключительно из любопытства, не помню точно, что меня сподвигло, но после недолгих поисков в интернете я наткнулся на довольно понятный алгоритм, вот и решил все по нему слепить (:
Интерпретатор: Python3
Скрипт несколько раз переписывался, в итоге он состоит из трех файлов:
- cities.py — описание городовточек, содержит экземпляры класса Cities с координатами и смещением по времени;
# coding=utf-8 from tools import * class City: def __init__(self, coordinates, timeOffset): self.latitude, self.longitude = [time2deg(x) for x in coordinates] self.timeOffset = timeOffset Moscow = City(coordinates=((55, 45, 20.83), (37, 37, 3.48)), timeOffset=+4)
- tools.py — всяческие инструменты, в основном для работы с координатами;
# coding=utf-8 from math import radians, degrees from datetime import timedelta def deg2rad(deg): return radians(deg) def deg2time(deg): result = list() rem = deg for d in [1, 60, 60]: resFrac = d * rem resInt = int(resFrac) rem = round(resFrac - resInt, 4) result.append(resInt) return result def rad2deg(rad): return degrees(rad) def rad2time(rad): return deg2time(rad2deg(rad)) def time2deg(timeDeg): res = 0 for t, d in zip(timeDeg, [1, 60, 3600]): res += t / d return round(res, 4) def time2rad(timeDeg): return deg2rad(time2deg(timeDeg)) def time2timeD(angTime): d, h = divmod(angTime[0], 24) m, s = angTime[1:] return timedelta(d, seconds=h * 3600 + m * 60 + s) def deg2timeD(deg): return time2timeD(deg2time(deg)) def rad2timeD(rad): return time2timeD(rad2time(rad)) def chSign(value): return -value if value > 0 else abs(value) def iterPair(iterable): for i in range(1, len(iterable)): yield iterable[i - 1], iterable[i] if __name__ == '__main__': print('Self test is running') print('Converting degrees to radians and angular time') deg = 180 rad = deg2rad(deg) angTime = deg2time(deg) print('tDegreest=>t{}ntRadianst=>t{}ntAng. timet=>t{}'.format(deg, rad, angTime)) print('Converting radians to degrees and angular time') rad = 3.14 deg = rad2deg(rad) angTime = rad2time(rad) print('tRadianst=>t{}ntDegreest=>t{}ntAng. timet=>t{}'.format(rad, deg, angTime)) print('Converting angular time to degrees and radians') angTime = [90, 50, 0] deg = time2deg(angTime) rad = time2rad(angTime) print('tAng. timet=>t{}ntDegreest=>t{}ntRadianst=>t{}'.format(angTime, deg, rad)) print('Converting all values to timedelta (human time)') angTimeD = time2timeD(angTime) degD = deg2timeD(deg) radD = rad2timeD(rad) print('tAng. timet{}t=>t{}ntDegreest{}t=>t{}ntRadianst{}t=>t{}'.format(angTime, angTimeD, deg, degD, rad, radD)) - sun.py — содержит все необходимые расчеты.
# coding=utf-8 from datetime import timedelta, date, datetime from math import floor, sin, cos, tan, asin, acos, atan from tools import * import cities class Sun: def __init__(self, city=cities.Moscow): #Словарь интервала дат + 2 даты (крайние для расчета солнечного дня), со значениями восходазахода self.dateRangeSun = dict() #Объект - город, содержащий параметры координат и часового смещения self.city = city #Углы зенита для разных положения солнца self.zenith = {'Offical': time2deg([90, 50, 0]), 'Civil': time2deg([96, 0, 0]), 'Nautical': time2deg([102, 0, 0]), 'Astronomical': time2deg([108, 0, 0])} #Главная функция def getResult(self, dates=date.today().strftime('%d.%m.%Y')): #Получаем диапазон дат, dstetime dates = self._datesRange(dates) #Для всез дат из диапазона for focusDate in dates: #Получаем параметры солнца self._getSun(focusDate) #Создания списка дат из диапазона def _datesRange(self, dates): dates = {datetime.strptime(focusDate, '%d.%m.%Y').date() for focusDate in dates.split('-')} if len(dates) == 2: dateBegin, dateEnd = sorted(dates) daysRange = dateEnd-dateBegin for i in range(daysRange.days+1): dates.add(dateBegin+timedelta(i)) return dates #Получение параметров солнца def _getSun(self, focusDate): #Список смещений для соседей directDisp = range(-1, 2) #Срез с данными мз 3х дат для расчета светового и темного дня section = list() #Для всех смещений for disp in directDisp: #Получаем дату со смещением относительно фокусной даты dispDate = focusDate + timedelta(disp) #Если такой даты в глобальном списке нет, то if not dispDate in self.dateRangeSun: #Определяем параметры солнца в глобальном словаре дат self.dateRangeSun[dispDate] = self._sunInfo(dispDate) #Попутно, составляем список-срез из 3х дат для расчета светлого и темного времени суток section.append([self.dateRangeSun[dispDate], chSign(disp)]) #ПРоверяем, если для фокусной даты, уже расчитаны эти параметры, то опускаем расчеты if not (self.dateRangeSun[focusDate].get('Light') and self.dateRangeSun[focusDate].get('Dark')): self.dateRangeSun[focusDate].update(self._lightDark(section)) #Общая функция вызова расчетов восхода и закаат солнца def _sunInfo(self, focusDate): return {'Sunrise': self._sunrise(focusDate), 'Sunset': self._sunset(focusDate)} #индивидуальные параметры для расчета восхода def _sunrise(self, focusDate): ft = lambda dayOfYear, lngHour: dayOfYear + ((6 - lngHour) / 24) fH = lambda cosH: 360 - rad2deg(acos(cosH)) fUT = lambda T, lngHour: T - lngHour #Вызываем калькулятор с индивидуальными параметрами return self._calculate(focusDate, ft, fH, fUT) #Индивидуальные параметры для расчета заката def _sunset(self, focusDate): ft = lambda dayOfYear, lngHour: dayOfYear + ((18 - lngHour) / 24) fH = lambda cosH: rad2deg(acos(cosH)) fUT = lambda T, lngHour: 24 - abs(T - lngHour) #Вызываем калькулятор с индивидуальными параметрами return self._calculate(focusDate, ft, fH, fUT) #Функция, для расчета светлого и емного времени суток def _lightDark(self, section): #Список - "отрезок", со всеми событиями за сегодняшний день (т.к. их возможно больше 2х) line = list() #Условимся, что изначально, у нас 24 часа - темно. (Одни сутки) lTime, dTime = timedelta(0), timedelta(1) #Возможные типы солнца sunTypes = ['Sunrise', 'Sunset'] #Разбераем срез. В нем, dateSunInfo - инф. о Солнце, и величина обратного смещения даты. for dateSunInfo, backDisp in section: #Перебераем типы, для навигации по словарю for sunType in sunTypes: #Получаем дату события, [тип][тип зенита], и пребавляем к ней обратное смещения focus = dateSunInfo[sunType]['Offical'] + timedelta(backDisp) #Прибавляем, для того, чтобы понять, возможно это событие произошло в нашу дату? #Например, есть дата со смещение -1 от нужной, прибавив к ней обратное смещение получим 0 #И если это так, то, вчерашнее событие, например - закат, произошел сегодня if focus.days == 0: #Собераем список сегодняшних событий line.append((sunType, dateSunInfo[sunType]['Offical'])) #Определяем начало и конец списка-отрезка #Если первая величина - закат, значит начало - восход, ибо до заката, всегда восход и наоборот first = [('Sunrise', timedelta(0)) if line[0][0] == 'Sunset' else ('Sunset', timedelta(0))] #То же самое для конца отрезка last = [('Sunrise', timedelta(1)) if line[-1][0] == 'Sunset' else ('Sunset', timedelta(1))] #Собераем отрезок line = first + line + last #iterPair - разберает отрезок на пары значений, со смещением=1 for pair in iterPair(line): #Получаем тип события pairType = [x[0] for x in pair] #И его значение pairVal = [x[1] for x in pair] #Если событие - Восход-Закат if pairType == ['Sunrise', 'Sunset']: #Значит - в это время было светл #ПРибавляем к свету, разницу во времени между Закатом и Восходом lTime += pairVal[1] - pairVal[0] #А темно было, ровно столько, сколько небыло светло (= dTime -= lTime return {'Light': lTime, 'Dark': dTime} #Функция-калькулятор, для расчета времени события def _calculate(self, focusDate, ft, fH, fUT): calculation = dict() dayOfYear = int(focusDate.strftime('%j')) latitude = self.city.latitude longitude = self.city.longitude lngHour = longitude / 15 t = ft(dayOfYear, lngHour) sunPos = (0.9856 * t) - 3.289 lngSun = sunPos + (1.916 * sin(deg2rad(sunPos))) + (0.020 * sin(deg2rad(2 * sunPos))) + 282.634 - 360 RA = rad2deg(atan(0.91764 * tan(deg2rad(lngSun)))) Lquadrant = floor(lngSun / 90) * 90 RAquadrant = floor(RA / 90) * 90 RA = (RA + (Lquadrant - RAquadrant)) / 15 sinDec = 0.39782 * sin(deg2rad(lngSun)) cosDec = cos(asin(sinDec)) for name, deg in self.zenith.items(): cosH = (cos(deg2rad(deg)) - (sinDec * sin(deg2rad(latitude)))) / (cosDec * cos(deg2rad(latitude))) if cosH > 1: localT = None elif cosH < -1: localT = None else: H = fH(cosH) / 15 T = H + RA - (0.06571 * t) - 6.622 UT = fUT(T, lngHour) localT = deg2timeD(UT + self.city.timeOffset) calculation[name] = localT return calculation if __name__ == '__main__': MoscowSun = Sun(cities.Moscow) MoscowSun.getResult('06.03.2013-07.03.2013') for d in MoscowSun.dateRangeSun: print(d) for s in MoscowSun.dateRangeSun[d]: if s in ['Dark', 'Light']: print(s, MoscowSun.dateRangeSun[d][s]) else: for t in MoscowSun.dateRangeSun[d][s]: print(s, t, MoscowSun.dateRangeSun[d][s][t])
При расчете значений на один день, скрипт считает параметры для трех дней (-1; 0; +1), это необходимо для расчетов протяженности светового дня.
На выходе мы имеем словарь:
{datetime.date(2014, 1, 1):
{'Dark': datetime.timedelta(0, 60743),
'Sunrise': {'Offical': datetime.timedelta(0, 35965),
'Nautical': datetime.timedelta(0, 30270),
'Civil': datetime.timedelta(0, 33178),
'Astronomical': datetime.timedelta(0, 27559)},
'Light': datetime.timedelta(0, 25657),
'Sunset': {'Offical': datetime.timedelta(0, 61622),
'Nautical': datetime.timedelta(0, 67311),
'Civil': datetime.timedelta(0, 64406),
'Astronomical': datetime.timedelta(0, 70021)}},
datetime.date(2013, 12, 31):
{'Sunrise': {'Offical': datetime.timedelta(0, 35973),
'Nautical': datetime.timedelta(0, 30266),
'Civil': datetime.timedelta(0, 33179),
'Astronomical': datetime.timedelta(0, 27552)},
'Sunset': {'Offical': datetime.timedelta(0, 61531),
'Nautical': datetime.timedelta(0, 67233),
'Civil': datetime.timedelta(0, 64322),
'Astronomical': datetime.timedelta(0, 69945)}},
datetime.date(2014, 1, 2):
{'Sunrise': {'Offical': datetime.timedelta(0, 35950),
'Nautical': datetime.timedelta(0, 30266),
'Civil': datetime.timedelta(0, 33170),
'Astronomical': datetime.timedelta(0, 27558)},
'Sunset': {'Offical': datetime.timedelta(0, 61695),
'Nautical': datetime.timedelta(0, 67373),
'Civil': datetime.timedelta(0, 64472),
'Astronomical': datetime.timedelta(0, 70080)}}} Добавлять особо нечего, если кому пригодится буду рад, если нужно подогнать вывод скрипта — помогу по мере возможностей.
