Привет! После красочных метеорных потоков мы плавно движемся к астрономической осени. В этом году она предвещает нам лунное затмение, соединение Венеры и Юпитера, а также полеты ярких рукотворных спутников. Мой сегодняшний рассказ – о том, как моделировать отражение света от таких спутников, и насколько яркие вспышки ожидают нас в этом октябре.
Вспышка Иридиума, первое фото своими руками – навелся не туда, затвор открыл поздно, горизонт завалил :)
С чего все началось
Весной на Гиктаймсе появился пост про спутник «Маяк». Это кубсат размером 30х10х10 см, основная задача которого после выхода на орбиту – раскрыть зеркальную пленку, которая будет отражать солнечный свет в сторону Земли. По замыслу авторов, это должно сделать Маяк одним из ярчайших спутников в ночном небе. Пикантности добавляет тот факт, что Маяк – первый в мире спутник, средства на который собирали краудфандингом, причем весьма успешно.
(картинка отсюда)
Хабр, несомненно, все еще торт: в комментариях разгорелась оживленная дискуссия по поводу реализуемости проекта, к которой подключились и разработчики. И если с конструкцией спутника все было в порядке, то заявления про «самую яркую звезду на небе» очень удивили. Например, из-за того, что пленочное зеркало редко получается хорошим. Или потому, что спутник должен крутиться вокруг своей оси со скоростью один оборот в секунду – тогда скорость движения солнечного зайчика по Земле оказывается огромной, и вспышки должны быть слишком короткими. Разработчики признались, что у них и правда есть проблемы с расчетом отражений, и предложили желающим помочь им с этим.
Первым аналитический расчет предложил Quiensabe. (Его пост с не очень оптимистичными результатами был опубликован на ГТ, но его удалили из-за неаккуратной ссылки на краудфандинг — впрочем, это уже другая история.) Такой расчет давал адекватную оценку, но не мог учесть массы факторов – кривизны поверхности Земли, изменение яркости вспышки по мере вращения зеркала и прочего. Нужна была аккуратная численная модель, которая считает яркость отражения с небольшим шагом по времени в зависимости от взаимного расположения спутника, наблюдателя, Солнца и всего прочего.
Задача казалась не слишком сложной – просто нужно было аккуратно разобраться во всей геометрии и ничего не напутать. Мотивировала и возможность посчитать что-нибудь эдакое космическое. Ну и, само собой, радость от написания собственного велосипеда – бесценна.
Алгоритм
Основная задача модели – рассчитать максимальную яркость вспышки для заданных параметрах спутника и заданного положения наблюдателя на Земле. При этом предсказание вспышек на несколько дней уходит на второй план – для этого есть, например, heavens-above.com. Логично будет взять один пролет спутника над наблюдателем, и уже для него вычислить яркость отражения с небольшим шагом по времени. После этого, варьируя начальные параметры, можно будет вычислить максимально возможную яркость вспышки.
Итак, шаг первый: чтобы увидеть вспышку, необходимо, чтобы спутник: а) пролетел над горизонтом и б) находился не в тени Земли. Шаг второй: если эти условия выполнены, можно начинать считать яркость вспышки. Да, очень вероятно, что спутник будет вращаться – в этом случае вспышки могут быть периодическими. (Мигающий спутник!) Поэтому в общем случае нужен шаг третий – искать отдельные вспышки, разделенные интервалами темноты. Выстраивается логика работы программы, которую реализуют три функции:
- Первая функция считает координаты спутника и наблюдателя в течение одного витка с большим шагом по времени (1-5 секунд) и проверяет, в какие моменты спутник виден над горизонтом и не находится в тени. Если спутник появляется над горизонтом, функция рассчитывает интервалы времени, когда он виден, и передает их второй функции.
- Вторая функция работает с маленьким шагом по времени (0.001 – 1 с) и на каждом шаге рассчитывает яркость вспышки, которую видит наблюдатель.
- Третья функция – самая некосмическая. Она находит отдельные вспышки в результатах работы второй функции и считает их параметры: яркость, длительность, момент прохождения максимума и так далее.
Физическая модель
Моделировать задачу будем в системе отсчета ECI (Earth-centered inertial) – ее начало совпадает с центром Земли, ось Z направлена на северный полюс, а Солнце находится в плоскости X-Z.
Почему именно ECI? Потому, что Солнце в ней неподвижно, наблюдатель вращается вместе с Землей вокруг оси Z, а орбита спутника не меняет своего положения в пространстве. Это значит, что координаты и спутника, и наблюдателя легко параметризуются через углы. Подробности – под спойлером.
Координаты наблюдателя задаются широтой fiobs и текущей долготой thetaobs.
По сути thetaobs является местным временем: полдень соответствует thetaobs = 0, а полночь – thetaobs = 180. Еще нужно знать местное время в начальный момент времени thetaobs0 и скорость вращения Земли вокруг своей оси. Тогда координаты наблюдателя в ECI выражаются так:
Орбиту спутника для облегчения расчетов можно считать круговой. В этом случае ее можно задать тремя параметрами: радиусом Rorb, наклонением fi и долготой восходящего узла thetaAN. Восходящий узел (ascending node) – это точка, в которой орбита пересекает экватор по направлению к северу.
Гравитация, бессердечная ты с. Так как орбита круговая, то спутник движется по ней с неизменной угловой скоростью
Теперь, когда мы знаем все это, положение спутника на орбите легко параметризовать через угол alpha:
Здесь dirAN и dirB – направления на восходящий узел и самую высокую точку орбиты, которые считаются простой тригонометрией:
Солнце. Фишка системы отсчета ECI в том, что Солнце всегда находится в плоскости X-Z. В зависимости от времени года оно меняет положение: когда в северном полушарии лето, Солнце находится выше оси Z, а зимой – наоборот. Нас интересует не само положение Солнца, а направление на него. Вернее, от него, то есть вдоль направления движения солнечных лучей. На картинках он называется ksun.
Посчитали координаты? Теперь надо проверить, виден ли спутник над горизонтом – иначе считать что-либо бессмысленно. Это определяется сонаправленностью векторов от наблюдателя на спутник -dirsat_obs и в зенит dirobs: если их скалярное произведение положительно, то спутник где-то над нами.
Если спутник над горизонтом, то нужно узнать его координаты на небе. Для этого направление на спутник надо спроецировать на плоскость наблюдателя. Плоскость наблюдателя задается векторами на север dirnord и запад dirwest. Первый находится через направление на зенит и северный полюс Земли dNP, а уже по нему и направлению в зенит вычисляется вектор на запад.
Теперь нужно проверить, не находится ли спутник в тени Земли (иначе отражения, очевидно, не будет). Из картинки ниже легко понять, что спутник находится в тени, если одновременно выполняются два условия:
- угол между ksun и dirsat (вектор на спутник) острый;
- d < R
Аналогичным способом надо проверить, что наблюдатель находится в тени, то есть у него сейчас ночь. На солнечном свету спутник увидеть будет непросто, хотя в дневной астрономии, несомненно, тоже есть свой шарм.
Наконец, если спутник освещен Солнцем и находится где-то над нами, остается посчитать яркость вспышки. Точнее, ее видимую звездную величину, которая показывает яркость по сравнению с другими звездами. Шкала звездных величин логарифмическая: она определяется как
то есть изменение яркости в 100 раз соответствует изменению звездной величины на 5 единиц. А еще отсюда видно, что шкала звездных величин инвертирована: более ярким объектам соответствуют меньшие значения. Например, звездная величина Арктура и Веги – ярчайших звезд Северного полушария – составляет +0.0, Сириуса – -1.7, а вот Венера достигает -4.6 звездной величины. Глаз различает звезды вплоть до +6 звездной величины, хотя в крупных городах засветка не дает увидеть что-то ярче +4.
Для расчета звездной величины нужно знать яркость, а вернее освещенность – то есть мощность светового потока, проходящая через единицу площади. Измеряется она в Вт/м2 или в люксах, то есть физический смысл освещенности в том, насколько ярко будет освещена поверхность. Из определения следует, что освещенность не зависит от площади. Поэтому наблюдение невооруженным глазом, фотоаппаратом или в телескоп дадут один и тот же результат измерения звездной величины.
Как посчитать освещенность, которую создает отражение от спутника? Прежде всего, мощность отраженного света (то есть сколько фотонов за секунду отражаются от спутника) находится перемножением площади зеркала на освещенность, создаваемую Солнцем. Нужно не забыть, что зеркало может быть под углом к солнечным лучам. А еще про то, что зеркало неоднородно и отражает свет в большой телесный угол – «зайчик» на поверхности Земли будет протяженным, и мы можем быть вдали от его центра. Наконец, освещенность квадратично зависит от расстояния до спутника – чем дальше спутник, тем меньше света достигнет наблюдателя. Конечное выражение для освещенности выглядит так:
где Eorbit – освещенность, создаваемая Солнцем на орбите Земли, второй сомножитель – площадь зеркала под углом к солнечным лучам, третий – доля телесного угла, который видит наблюдатель, четвертый – квадратичная зависимость от расстояния до спутника.
Сравнивая это число освещенностью, соответствующей нулевой звездной величине (2,5∙10-8 Вт/м2), находим звездную величину спутника.
Увидим мы вспышку или нет, определяется положением зеркала. Его удобно задавать вектором нормали nmir: это позволяет легко посчитать направление отраженного луча по старому-доброму правилу про угол падения, равный углу отражения:
Легко и просто? Отнюдь! На практике быстрее считать, как должно располагаться зеркало, чтобы при заданном положении спутника и наблюдателя отражение падало бы прямо на наблюдателя. И затем сравнивать с текущим положением зеркала. Преимущество такого подхода в том, что он позволяет изменять шаг по времени. Пока зеркало далеко от оптимального положения, шаг можно увеличить – вспышки в ближайшее время все равно не будет. Если же вспышка уже началась, то шаг по времени нужно уменьшить, чтобы посчитать ее с максимальным разрешением.
Наконец, спутник вращается, а вместе с ним вращается и зеркало. Поэтому нормаль к зеркалу nmir удобно разложить на две компоненты – параллельную и перпендикулярную оси вращения спутника dmir. Первая остается неподвижной, вторая прецессирует – математика этого точь-в-точь такая же, как и движение спутника по орбите.
Код
Код написан на матлабе. Как я рассказал выше, он разбит на три функции, которые по очереди вызываются основным скриптом calc_one. В основном скрипте задаются константы, начальные условия и параметры симуляции, которые передаются функциям через структуру param следующим образом:
% satellite orbit
param.incl_sat = 97.75; % inclination, deg (ISS: 51.6, SSO for 600 km: 97.75)
param.theta_asc = -7.5; % longitude of the ascending node, deg
param.alpha_sat_0 = 0; % angle between AN and the sat in the plane of orbit at t=0, deg
param.h_orbit = 600 * 1E3; % orbit height, m%%%%%%%%%%%%%%%%%% SIMULATION PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%
% world
param.rad_earth = 6400 * 1E3; % Earth radius, m
param.GM = 6.67*1E-11 * 5.972*1E24; % grav. constant times Earth mass, m^3/s^2
param.omega_obs = 1/(60*60*24); % Earth angular velocity, turns/s
param.k_sun = [-0.99, 0, 0.16]; % sunlight wavevector (22 Jun: [-0.92,0,-0.4]; 21 Dec: [-0.92,0,0.4], 15 Oct: [-0.99,0,0.16])
param.atmosphere_transmission = 0.7; % atmosphere transmission
param.flux_sun = 1367; % Sun irradiant flux on the Earth orbit, W/m^2
param.flux_m0 = 2.5*1E-8; % flux corresponding to the apparent magnitude = 0, W/m^2
param.min_m = 6; % minimal visible apparent magnitude
% observer
param.phi_obs = 55.75; % latitude, deg (Moscow: 55.75)
param.theta_obs_0 = 170; % longitude, deg (noon: 0, midnight: 180)
% satellite orbit
param.incl_sat = 97.75; % inclination, deg
param.theta_asc = -7.5; % longitude of the ascending node, deg
param.alpha_sat_0 = 0; % angle between AN and the sat in the plane of orbit at t=0, deg
param.h_orbit = 600 * 1E3; % orbit height, m
% satellite mirror
param.omega_rot = 1; % rotation speed of the satellite, turns/s
param.d_rot = [0 0.7 1]; % rotation axis of the satellite
param.n_mirror_0 = [1 0 0]; % normal to the mirror at t=0
param.gauss_w = 17.2; % divergence angle of a normal distrubution, deg
param.sq_mirror = 3.8; % mirror square, m^2
% simulation parameters
param.dt1 = 1; % time step if the satellite is below horizon or in the shadow, s
param.dt2 = 1; % same if satellite is above horizon, but reflection is far away
param.dt3 = 1 * 1E-2; % same if satellite is above horizon, reflection is close
param.dt4 = 5 * 1E-3; % same if reflection is seen
param.frame_duration = 0.03; % exposition time of a camera/eye, sПервая функция fn_pass считает координаты спутника и ищет один пролет над горизонтом. Вначале она вычисляет всякие вспомогательные величины (направление на северный полюс, скорость вращения спутника и так далее).
% calculate some aux values
rad_orbit = rad_earth + h_orbit; % radius of the orbit from the center of Earth
r_np = [0 0 rad_earth]; % coordinates of the North pole
dir_an = [cosd(theta_asc), sind(theta_asc), 0]; % direction on the acsending node
dir_b = [-sind(theta_asc)*cosd(incl_sat),...
cosd(theta_asc)*cosd(incl_sat), sind(incl_sat)]; % direction on the highest point of the sat
omega_sat = sqrt(2*GM/rad_orbit)/(2*pi*rad_orbit); % satellite angular velocity on the orbit, s^-1
dir_obs_z = sind(phi_obs); % z-coordinate of the observer, does not changeПосле этого запускается while-цикл по времени с постоянным шагом, равным 1 секунде. Условие выхода из цикла – уход спутника за горизонт (или же два оборота спутника без единого появления над горизонтом, что нереально на низких орбитах). В цикле последовательно рассчитываются координаты наблюдателя и спутника:
% coordinates of the observer
theta_obs = mod(theta_obs_0 + 360*omega_obs*current_time, 360);
dir_obs = [cosd(phi_obs)*cosd(theta_obs), cosd(phi_obs)*sind(theta_obs), dir_obs_z];
r_obs = rad_earth * dir_obs;
% coordinates of the satellite
theta_sat = mod(alpha_sat_0 + 360*omega_sat*current_time, 360);
dir_sat = cosd(theta_sat)*dir_an + sind(theta_sat)*dir_b;
r_sat = rad_orbit*dir_sat;
% vector from the satellite to the observer
r_sat_obs = r_obs-r_sat;
dir_sat_obs = fn_norm(r_sat_obs); проверяется, что спутник находится над горизонтом и рассчитываются его координаты на небе:
% CHECK1: check if the satellite is above the horizon
if fn_dot(dir_obs, dir_sat_obs) < 0 % basis on the skychart
r_obs_np = r_np - r_obs; % vector from the observer to the North pole
dist_obs_np = sqrt(sum(r_obs_np.^2));
cos_angle_obs_np = fn_dot(r_obs_np, -dir_obs)/dist_obs_np;
dir_nord = r_obs_np + cos_angle_obs_np*dist_obs_np*dir_obs;
dir_nord = fn_norm(dir_nord); % vector to the nord
dir_west = fn_cross(dir_obs, dir_nord);
dir_west = fn_norm(dir_west); % vector to the west
% coordinates of the satellite on the sky chart
pass_path.alt(pass_step,1) = 90 - acosd(fn_dot(-dir_sat_obs, dir_obs));
map_proj_nord = fn_dot(-dir_sat_obs, dir_nord);
map_proj_west = fn_dot(-dir_sat_obs, dir_west);
pass_path.azimuth(pass_step,1) = mod(-atan2d(map_proj_west, map_proj_nord), 360);и, наконец, проверяется, находится ли спутник в тени:
% CHECK 2: check that the satellite is not in the shadow
cos_angle_sat_ksun = fn_dot(dir_sat,k_sun);
dist_sat_ksun = rad_orbit*sqrt(1-cos_angle_sat_ksun^2); % distance from the sat to the central axis of shadow
if (dist_sat_ksun > rad_earth)
% satellite is not in the shadowВсе данные сохраняются в структуре pass_path, которую функция и возвращает. Но перед этим нужно найти все интервалы, в течение которых спутник находится над горизонтом и освещен Солнцем – они понадобятся для работы второй функции.
% STEP 2: if the pass was above horizon,
% select the regions when sat is above horizon and not in shadow
pass_path.reflection_intervals = [];
if isempty(pass_path.time) == 0
flag_in_shadow = 1;
num_regions_not_in_shadow = 0;
for counter_step = 1:size(pass_path.time,1)
if pass_path.in_shadow(counter_step) == 0
% satellite is not in shadow
if flag_in_shadow == 0
% sat is still not in shadow
else
% sat was in shadow and just went out of it
flag_in_shadow = 0;
num_regions_not_in_shadow = num_regions_not_in_shadow + 1;
pass_path.reflection_intervals(num_regions_not_in_shadow,1) = pass_path.time(counter_step);
end
else
% satellite is in shadow
if flag_in_shadow == 0
% sat was not in shadow and just went into it
flag_in_shadow = 1;
pass_path.reflection_intervals(num_regions_not_in_shadow,2) = pass_path.time(counter_step);
else
% sat is still in shadow
end
end
end
% if at the end sat was still not in shadow,
% time of the end of the last interval will be time of the last
% point above the horizon
if flag_in_shadow == 0
pass_path.reflection_intervals(num_regions_not_in_shadow,2) = pass_path.time(counter_step);
end
endИнтервалы сохраняются в pass_path.reflection_intervals в массиве Nx2, где N – количество интервалов видимости, 2 числа – время начала и конца интервала видимости. Вот теперь структура pass_path возвращается в главный скрипт, откуда вызывается вторая функция.
Вторая функция fn_reflection работает только в интервалах, где спутник находится над горизонтом и освещен Солнцем. На каждом шаге (который может быть короче одной секунды) он также рассчитывает координаты наблюдателя и спутника, после чего проверяет, находится ли зеркало близко к оптимальной позиции:
% CHECK 3: if angle (d_rot, n_opt) is close to angle (d_rot, n_mir0), proceed
% this means that geometry can in principle give a reflection
if (angle_drot_nopt > (angle_drot_nmirror-angle_div))...
&& (angle_drot_nopt < (angle_drot_nmirror+angle_div))
% current orientation of the mirror
n_mir = n_mir_parallel...
+ n_mir_perp0*cosd(360*omega_rot*current_time)...
+ n_mir_perp1*sind(360*omega_rot*current_time);
% angle between current orientation of the mirror and n_opt
angle_nmir_nopt = acosd(fn_dot(n_mir, n_opt));
% CHECK 4: if n_mir is close to n_opt, proceed and calculate the reflection
if angle_nmir_nopt < angle_divКак видно из кода, расположение зеркала проверяется по двум параметрам – азимутальному (CHECK 3) и радиальному (CHECK 4) углам по отношению к оси вращения спутника. Фишка в том, что радиальный угол меняется при вращении спутника. А вот азимутальный – нет!
Поэтому сначала имеет смысл проверить, находится ли азимутальный угол близко к оптимальному положению. Если нет, то в ближайшее время ничего не поможет. Если да, то можно посчитать и радиальный угол, а заодно уменьшить шаг по времени.
Наконец, если зеркало близко к оптимальной позиции, то можно посчитать и освещенность, которую создает солнечный зайчик. Ну и звездную величину заодно
% calculate the angle between the direction on the observer
% and direction of the main reflection
cos_angle_nmir_ksun = fn_dot(n_mir, -k_sun);
k_refl = 2*cos_angle_nmir_ksun*n_mir + k_sun;
cos_angle_refl_obs = abs(fn_dot(k_refl, dir_sat_obs));
angle_refl_obs = acosd(cos_angle_refl_obs);
% flux at the observer's place
refl_fraction = 2/(pi*degtorad(gauss_w)^2) * exp (-2*angle_refl_obs^2/gauss_w^2); % now it is Gaussian
int_reflection = flux_sun * sq_mirror * cos_angle_nmir_ksun;
flux = max (int_reflection * refl_fraction ...
* atmosphere_transmission / distance^2, flux_min); % W/m^2Видно, что если яркость слишком мала (тусклее +6), то в результат запишется +6. Это значит, что мы ничего не увидим. Все результаты записываются в структуру pass, которую функция и возвращает.
Третья функция fn_flares самая простая: она считает вспышки, их максимальную яркость, длительность и расположение на небе. Она реально примитивная, но длинная, поэтому я убрал ее под спойлер.
function pass = fn_flares(pass, param)
% this functions counts the flares into the 'pass' structure
% summarizes data about their magnitudes, durations, etc.
% and saves them in 'pass.flares' field
% initialization before the main loop
counter_flare = 0;
flare_is_now = 0;
clear flares
flares = struct('time',[],'dist',[],'alt',[],'azimuth',[],'inst_magn',[],'vis_magn',[],'dur',[]);
% main loop
for step = 1:length(pass.time)
if pass.magnitude(step) < param.min_m
% flare is present now
if flare_is_now == 0
% flare just began
flare_is_now = 1;
counter_flare = counter_flare + 1;
clear current_flare
current_flare.time_beg = pass.time(step);
% integrate the energy of the flare, [J/m^2]
if step < length(pass.time)
current_flare.energy = pass.flux(step)*...
(pass.time(step+1) - pass.time(step));
else
% this was the last data point, nothing to add
end
% create field 'current_flare.brightest'
% with the data about the brightest point of the flare
current_flare.brightest.magn = pass.magnitude(step);
current_flare.brightest.time = pass.time(step);
current_flare.brightest.azimuth = pass.azimuth(step);
current_flare.brightest.alt = pass.alt(step);
current_flare.brightest.distance = pass.distance(step);
else
% flare still goes on
% increment the energy of the flare
if step < length(pass.time)
current_flare.energy = current_flare.energy + pass.flux(step)*...
(pass.time(step+1) - pass.time(step));
else
% this was the last data point, nothing to add
end
% if here flare is brighter, renew the data about the brightest point
if pass.magnitude(step) < current_flare.brightest.magn
current_flare.brightest.magn = pass.magnitude(step);
current_flare.brightest.time = pass.time(step);
current_flare.brightest.azimuth = pass.azimuth(step);
current_flare.brightest.alt = pass.alt(step);
current_flare.brightest.distance = pass.distance(step);
end
end
else
% no flare now
if flare_is_now == 0
% still no flare
else
% flare just ended, sum up the results
flare_is_now = 0;
current_flare.time_end = pass.time(step-1);
current_flare.duration = current_flare.time_end - current_flare.time_beg;
if current_flare.duration > param.frame_duration
% if the flare is long, we can see it directly
current_flare.vis_magn = current_flare.brightest.magn;
else
% if the flare is short, one should divide its energy
% over the exposure time of an eye/camera to get the
% apparent magnitude
current_flare.vis_magn = ...
min(-2.5*log10(current_flare.energy/param.frame_duration/param.flux_m0), param.min_m);
end
% save the data in the 'flares' structure
flares.time(counter_flare,1) = current_flare.brightest.time;
flares.dist(counter_flare,1) = current_flare.brightest.distance;
flares.azimuth(counter_flare,1) = current_flare.brightest.azimuth;
flares.alt(counter_flare,1) = current_flare.brightest.alt;
flares.inst_magn(counter_flare,1) = current_flare.brightest.magn;
flares.vis_magn(counter_flare,1) = current_flare.vis_magn;
flares.dur(counter_flare,1) = current_flare.duration;
end
end
end
% return the result
flares.num_of_flares = counter_flare;
pass.flares = flares;
endЗдесь появляется одна проблема. Если вспышка очень короткая, то суммарное число попавших в глаз фотонов может быть низким даже если моментальная яркость была большой. Тогда видимая яркость вспышки будет гораздо меньше максимальной.
Что значит «короткая вспышка»? Очевидно, такая, которая короче времени экспозиции камеры. А если мы наблюдаем ее глазом? Что считать экспозицией глаза, я не имел ни малейшего понятия, поэтому спросил у Meklon. Он рассказал много всего интересного – но сводилось это к тому, что глаз работает ооочень сложно, чувствительность у него сильно нелинейная, и поэтому четкого ответа на вопрос дать нельзя. Даже примерно. В итоге было решено выбрать какое-нибудь более-менее осмысленное число — я выбрал 30 мс, что соответствует примерно 30 кадрам в секунду.
Итак, если вспышка оказывается короче времени экспозиции, то ее длительность округляется вверх до времени экспозиции, суммарное количество света остается неизменным, ну а яркость, соответственно, падает. Например, если в реальности вспышка длится 10 мс и за это время прилетело 3000 фотонов, то ее моментальная «яркость» составляет 300 фотонов/мс. Мы же увидим, что она продолжалась 30 мс с «яркостью» 100 фотонов/мс – в 3 раза меньше моментальной.
Наконец, главный скрипт calc_one сначала вызывает первую функцию fn_pass, и если спутник виден над горизонтом и не находится в тени, то считает яркость отражений функцией fn_reflection и параметры вспышек функцией fn_flares:
% calculate trajectory of the pass
pass_path = fn_pass(param);
% proceed if the satellite is above the horizon and not always in shadow
if isempty(pass_path.reflection_intervals) == 0
% calculate the reflections and magnitudes of flares
pass = fn_reflection(param, pass_path.reflection_intervals);
pass = fn_flares(pass, param);
endОптимизация
Первый прием оптимизации я уже упомянул: это переменный шаг по времени. Функция fn_pass работает с большим шагом (1 с) – для проверки того, виден ли спутник над горизонтом, этого более, чем достаточно. А вот у функции fn_reflection целых три возможных шага по времени: самый маленький (до 1 мс) нужен для расчета профиля вспышки, два других – когда вспышка не видна, но вскоре может появиться. Они определяются в структуре с параметрами:
param.dt1 = 1; % time step if the satellite is below horizon or in the shadow, s
param.dt2 = 1; % same if satellite is above horizon, but reflection is far away
param.dt3 = 1 * 1E-2; % same if satellite is above horizon, reflection is close
param.dt4 = 5 * 1E-3; % same if reflection is seenВторой момент связан со встроенными матлабовскими функциями. Они задумывались универсальными – скажем, один и тот же косинус можно считать и от вещественных, и от комплексных чисел, как одиночных, так и в составе массивов. Естественно, что для этой универсальности приходится проверять входные параметры. В случае простых функций типа скалярного произведения эта проверка может занимать больше времени, чем собственно математика:
function c = dot(a,b,dim)
if isinteger(a) || isinteger(b)
error(message('MATLAB:dot:integerClass'));
end
% Special case: A and B are vectors and dim not supplied
if ismatrix(a) && ismatrix(b) && nargin<3
if min(size(a))==1, a = a(:); end
if min(size(b))==1, b = b(:); end
end;
% Check dimensions
if any(size(a)~=size(b)),
error(message('MATLAB:dot:InputSizeMismatch'));
end
if nargin==2,
c = sum(conj(a).*b);
else
c = sum(conj(a).*b,dim);
end
В нашем коде чаще всего вызываются как раз векторное и скалярное произведения, ну и еще нормировка вектора. Поэтому имеет смысл написать для них свои функции без проверок — ведь все векторы заведомо трехмерны, а их компоненты вещественны:
function c = fn_dot(a, b)
c = a(1)*b(1) + a(2)*b(2) + a(3)*b(3);
endПосле этого скорость вычислений выросла примерно вдвое.
Вообще в матлабе есть такая фишка как анонимные функции. Идея в том, чтобы объявить функцию прямо в тексте другой функции, не занимая дополнительный файл. Утверждается, что они делают код более удобочитаемым и работают быстрее, нежели отдельные функции.
Оказалось, что на деле все наоборот. Когда я пытался объявлять скалярные и векторные произведения в виде анонимных функций, вычисления занимали заметно больше времени. В тонкостях работы компилятора матлаба я не силен, поэтому было бы здорово, если бы кто-нибудь мог бы объяснить происходящее.
Код можно найти на гитхабе. Расчет вспышек проверялся на Иридиуме и давал осмысленные результаты – яркость вспышек составляла около -9 звездной величины.
Маяк
Спутник «Маяк» планируется запустить в этом октябре на круговую орбиту высотой 600 км. В средней полосе России его будет лучше всего видно примерно с 22:00 до 02:00 по местному времени. Пленочные зеркала Маяка отражают свет гораздо хуже антенн Иридиума, зато заметно превосходят их по площади.
| Параметр | Значение |
|---|---|
| Высота орбиты | 400 – 600 км |
| Наклонение орбиты | 97.7° |
| Направление на восходящий узел | -7.5° |
| Площадь зеркала | 3.8 м2 |
| Расходимость отражения (1σ) | 17.2° |
| Скорость вращения | неизвестна (ожидается 0.5 – 1 об/с) |
| Ориентация оси вращения | неизвестна |
После выхода на орбиту спутник должен раскрыть зеркала и закрутиться со скоростью около 1 оборота в секунду. Проблема в том, что направление оси закрутки неизвестно в принципе – на спутнике нет ни систем ориентации, ни телеметрии, краудфандинг их бы просто не потянул. Ко всему этому добавляется закрутка при отделении кубсата от носителя. Вот, скажем, два кубсата, запущенные с МКС; видно, что сразу же после запуска они уже вращаются по-разному:
Остается задать случайное направление оси вращения и ориентации зеркала в скрипт calc_one и посмотреть, увидит ли наблюдатель пролет спутника:
Ух ты, что-то видно! Слева – траектория полета спутника в московском небе, примерно в 11 часов вечера. Красные точки – видимые вспышки спутника, черная линия – спутник в тени Земли и не виден. Графики справа – это яркость вспышек (сверху) и их длительность (снизу) в зависимости от времени. Видно, что вспышки бывают яркими (ярче Веги), но короткими. А еще судя по обилию точек их реально очень много (скрипт насчитал 330):
По вспышке в секунду! Капитан Очевидность подсказывает, что так и должно быть для спутника, закрученного со скоростью 1 об/с. А профиль одиночной вспышки выглядит так:
Вспышка действительно короткая и обрывается на +6 звездной величине, когда она становится тусклее фона.
Главный вопрос, который нас волнует – а когда лучше всего наблюдать вспышки? Для этого можно запустить скрипт для разных положений (= разного местного времени) наблюдателя. Как быть с тем, что мы не знаем ориентацию спутника? Можно каждый раз проводить вычисления для нескольких случайных ориентаций зеркала и оси вращения, после чего выбирать самую яркую вспышку. Код изменится примерно так:
line_theta_obs_0 = 134:2:236; % values of theta_obs_0 to be simulated
number_of_iterations = 50; % number of different mirror orientations for one data point
for counter_theta_obs_0 = 1:size(line_theta_obs_0,2)
param.theta_obs_0 = line_theta_obs_0(counter_theta_obs_0);
% calculate trajectory of the pass
pass_path = fn_pass(param);и если спутник пролетает на горизонтом,
% loop over different random orientations of the mirror
for counter_mirror = 1:number_of_iterations
% create randomly oriented mirror and rotation axis
param.d_rot = [rand, rand, rand];
param.n_mirror_0 = [rand, rand, rand];
% calculate the reflections and magnitudes of flares
pass = fn_reflection(param, pass_path.reflection_intervals);
pass = fn_flares(pass, param);
endЗапускаем его для 50 положений наблюдателя и 50 случайных ориентаций спутника, ждем полночи (да, времени на симуляцию уходит много), получаем результат:
Самые яркие вспышки приходятся на полночь – довольно удобное время для наблюдений. Вечерние вспышки тоже могут быть яркими. К сожалению, светить ярче Сириуса и планет не получится, но при удачном раскладе видно спутник будет довольно неплохо. Видно, что яркость вспышек падает по мере удаления спутника от зенита (нижний синий график), что связано с увеличением расстояния до спутника, когда он летит вблизи горизонта. При этом длительность вспышек (средний фиолетовый график) практически не меняется.
Какова вероятность увидеть яркую вспышку? Единственный случайный элемент в задаче – это ориентация спутника на орбите. В цикле по ориентациям спутника стоило сохранять яркость вспышки для каждой итерации, а потом строить гистограмму. (Черт, надо было сразу об этом подумать!) Добавляем нужный код и запускаем скрипт для 100 ориентаций спутника. Наутро получаем результат:
Здесь разным цветом показаны вспышки разной яркости. Столбик на диаграмме – конкретное местное время; вероятности вспышек всевозможных яркостей в нем составляют 100%. Видно, что результаты получаются неутешительные:
- вспышки ярче Веги (звездная величина +0) можно увидеть где-то в течение часа с 11:15 до 00:15; вероятность этого около 10%;
- вспышки ярче звезд Большой медведицы и Кассиопеи (ярче +2) видны чуть дольше двух часов (с 10 вечера до половины первого), вероятность их увидеть около 20%;
- с вероятностью более 50% вспышек не будет вообще.
Теперь хорошие новости: зеркал на спутнике три. Это не значит, что вероятность вырастет втрое (хотя бы потому, что тогда мы получали бы бессмысленные вероятности в 150%). Но при удачном расположении зеркал увеличение вероятности в два раза вполне реально, что весьма неплохо. Если очень повезет, то можно будет наблюдать вспышки от двух зеркал по очереди. Да, на яркость это никак не повлияет – увидеть что-либо ярче Веги можно будет только в «звездный час» около полуночи.
Давайте еще что-нибудь посчитаем. Например, как будет изменяться яркость вспышки в зависимости от скорости закрутки спутника.
Если спутник вращается слишком быстро, то вспышки становятся слишком короткими (нижний график). Тогда камера/глаз округлит время экспозиции до минимально возможного (1 кадр или 0.03 с). Во столько же раз упадет яркость вспышки. В нашем случае это происходит при закрутке 3 об/с.
А если скорость вращения наоборот, слишком маленькая? Тоже ничего страшного вплоть до ~0.01 об/с. А вот ниже вспышки будут идти слишком редко – может сложиться так, что первая вспышка закончится задолго до зенита, а вторая начнется сильно после. Насколько мне известно, «Маяк» будет вращаться со скоростью 0.1 – 1 об/с, что далеко и от верхнего, и от нижнего пределов – а значит, при любой закрутке мы должны видеть частые яркие вспышки.
Ну и теперь о грустном. Век Маяка очень недолог – его огромные зеркала будут тормозиться об остаточную атмосферу даже на высоте 600 км. На протяжении месяца спутник будет снижаться до 400 км, где торможение еще сильнее, и в течение считанных дней после этого он войдет в плотные слои атмосферы. Но при этом он станет ближе к нам, а значит, и ярче аж на 0.6 звездной величины:
Вместо послесловия
Несомненно, код можно было улучшать – добавить GUI, эллиптичность орбиты, количество зеркал, параллельные вычисления на CUDA и много всего прочего. Впрочем, совершенству нет предела. К тому же нынешняя версия успешно решает все поставленные перед нею задачи.
Наконец, хочется сказать спасибо Quiensabe за плодотворные дискуссии в процессе разработки модели. А команде Маяка – спасибо за занимательную задачу и, конечно же, успешного полета!
Автор: qbertych
