Вот на 12м году работы программистом микроконтроллеров мне наконец-то пригодились численные методы из ВУЗ(овского) курса по математике.
Причем не сколько для работы сколько для pet-проекта на котором я отлаживаю накопленную кодовую базу.
Я делаю солнечный навигатор по фазам освещения. Про это у меня есть отдельный подробный текст. Там надо, в частности, по измеренной секундомером кварцевого RTC продолжительности светового дня вычислить географическую широту наблюдения. Вычисления, конечно же, следует производить с учетом рефракции Солнца в атмосфере.
Для этого давно существует формула (1). Это трансцендентное уравнение. Его вывод можно найти в любом советском учебнике астрономии для моряков.
тут
Обозначение | Назначение переменной | Ед. измерения |
Day | Продолжительность светового дня | часы |
phi | Географическая широта | радианы |
delta | Солнечное склонение | радианы |
pi | число пи | радианы |
Вот график этой функции для солнечного склонения -18.74 градусов. Это 26 января каждого года.
В математике бывают случаи, когда есть функция, которая просто выражается элементарными функциями. А вот обратную функцию выразить аналитически либо очень трудно, либо вообще невозможно. Либо аналитическое решение есть, но формула уж очень громоздкая. Это как раз похоже на этот случай.
Вот простой Си-код для вычисления формулы (1) на микроконтроллере.
double calc_day_length_duration_refraction_h(double phi_deg, double delta_deg) {
double day_length_duration_h = 0.0;
double phi_rad = DEG_2_RAD(phi_deg);
double delta_rad = DEG_2_RAD(delta_deg);
double refrac_rad = DEG_2_RAD(51.0 / 60.0);
double numerator = -sin(refrac_rad) - sin(phi_rad) * sin(delta_rad);
double denominator = cos(phi_rad) * cos(delta_rad);
double cos_arg = numerator / denominator;
if(cos_arg<=1.0){
if(-1.0<=cos_arg){
day_length_duration_h = (24.0 * acos(cos_arg)) / M_PI;
LOG_DEBUG(PLANETARIUM, "Phi:%7.2f Deg,Delta:%7.2f Deg,Dur:%5.2f h", phi_deg, delta_deg,day_length_duration_h);
}else{
day_length_duration_h = 25.0;
LOG_ERROR(PLANETARIUM, "AcosSmallErr %f", cos_arg);
}
}else{
day_length_duration_h = 0;
LOG_ERROR(PLANETARIUM, "AcosBigErr %f", cos_arg);
}
return day_length_duration_h;
}Однако для задачи навигации надо по известной продолжительности светового дня вычислить широту
Фоторезистор измерил мне долготу дня 8 часов 14 минут. Какая же широта соответствует этой продолжительности? Эту задачу можно решить графически. Взять рейсшину, угольник и карандашом провести две линии.
Однако как решить эту задачу на микроконтроллере? Внутри микроконтроллера нет бумаги, рейсшины и угольника. Надо написать код на языке Си.
Можно конечно рассчитать формулу (1) для всех широт с шагом нужной нам точности (0.001 градус), но этот код будет очень долго исполняться на микроконтроллере. Трудоёмкость перебора огромна. Это не наши методы. Особенно когда частота ядра специально занижена для увеличения срока службы батареи.
Вот тут-то нам как раз и помогут численные методы. Заметьте, что функция (1) монотонная. При отрицательном солнечном склонении функция убывает. Это хорошая новость. Можно воспользоваться бинарным поиском, чтобы найти нужную нам широту(X) по долготе для (Y).
Это функция вычисление широты по известной продолжительности дня, вычисленная методом бинарного поиска.
double planetarium_day_dur_to_phi(double day_length_duration_h, double delta_deg){
double phi_deg = -180.0;
LOG_INFO(PLANETARIUM, "CalcPhiRefr: DayDur: %f H, SunDec: %f Deg", day_length_duration_h, delta_deg);
if(0.0 < day_length_duration_h) {
bool res = is_double_equal_absolute(0.0, delta_deg, 0.001);
if(false == res) {
double cur_phi_deg = 0.0;
double max_phi_deg = 88.0;//1.2*ARCTIC_CIRCLE_DEG;
double min_phi_deg = -88.0; //1.2*ANTARCTIC_CIRCLE_DEG;
double cur_day_length_duration_diff_abs_h=0.0;
double prev_day_length_duration_diff_abs_h=0.0;
double cur_day_length_duration_diff_h=0.0;
double cur_day_length_duration_h=0.0;
uint32_t step = 0;
for(;;) {
step++;
cur_phi_deg = AVERAGE_2(min_phi_deg,max_phi_deg);
cur_day_length_duration_h = calc_day_length_duration_refraction_h(cur_phi_deg, delta_deg);
cur_day_length_duration_diff_h = day_length_duration_h -cur_day_length_duration_h;
LOG_DEBUG(PLANETARIUM, "Step:%u,Phi:[%7.3f..%7.3f..%7.3f] Deg->CurDur:%f h,Diff:%f h,",
step,
min_phi_deg, cur_phi_deg,max_phi_deg,
cur_day_length_duration_h,
cur_day_length_duration_diff_h
);
cur_day_length_duration_diff_abs_h = fabs(cur_day_length_duration_diff_h);
if(cur_day_length_duration_diff_abs_h<prev_day_length_duration_diff_abs_h){
}
if(cur_day_length_duration_diff_abs_h<MIN_2_HOUR(1.0/60.0)) {
phi_deg = cur_phi_deg;
LOG_INFO(PLANETARIUM, "SpotSol,DayLen:%f h,Phi:%f Deg,Steps:%u",day_length_duration_h, phi_deg,step);
break;
}else{
if(delta_deg<0){
if(0.0<cur_day_length_duration_diff_h){
max_phi_deg = cur_phi_deg;
}else {
//cur_day_length_duration_diff_h<0.0
min_phi_deg = cur_phi_deg;
}
}else{
//0<delta_deg
if(0.0<cur_day_length_duration_diff_h){
min_phi_deg = cur_phi_deg;
}else {
//cur_day_length_duration_diff_h<0.0
max_phi_deg = cur_phi_deg;
}
}
}
prev_day_length_duration_diff_abs_h=cur_day_length_duration_diff_abs_h;
}
}
}
return phi_deg;
}
Функция planetarium_day_dur_to_phi перестает работать как только приближение к решению меньше приемлемой для нас погрешности. В данном случае нас устраивает, если мы найдем решение с погрешностью +/- 1на угловая минута.
Отладка
Для исходных данных возьмем продолжительность светового дня 8.23 часов и солнечное склонение -18.742 deg, алгоритм нашел значение широты 55.907...+/-0.016 градусов. Буквально 15 итераций потребовалось на то, чтобы решить уравнение (1).
Плюсы численных методов:
++Можно достигнуть любую наперед заданную точность.
++Меньше вычислений, чем в методе перебора.
Минусы численных методов:
--Больше вычислений чем в вычислении формулы.
Итоги
Численные методы позволяют получать решения уравнений, для которых нет аналитического решения с любой точностью.
Вот так просто с помощью фоторезистора за один световой день можно выдать навигационное решение.
Надеюсь этот текст сподвигнет кого-нибудь тоже определить свою широту по датчику естественного освещения.
