Привет, Хабр! После красочных метеорных потоков мы плавно движемся к астрономической осени. В этом году она предвещает нам лунное затмение, соединение Венеры и Юпитера, а также полеты ярких рукотворных спутников. Мой сегодняшний рассказ – о том, как моделировать отражение света от таких спутников, и что необычного ожидает нас в этом октябре.

С чего все началось

Весной на Гиктаймсе появился пост про спутник «Маяк». Это кубсат размером 30х10х10 см, основная задача которого после выхода на орбиту – раскрыть зеркальную пленку, которая будет отражать солнечный свет в сторону Земли. По замыслу авторов, это должно сделать Маяк одним из ярчайших спутников в ночном небе. Пикантности добавляет тот факт, что Маяк – первый в России спутник, средства на который собирали краудфандингом, причем весьма успешно.

(картинка отсюда)

Хабр, несомненно, все еще торт: в комментариях разгорелась оживленная дискуссия по поводу реализуемости проекта, к которой подключились и разработчики. И если с конструкцией спутника все было в порядке, то заявления про «самую яркую звезду на небе» очень удивили. Например, из-за того, что пленочное зеркало редко получается хорошим. Или потому, что спутник должен крутиться вокруг своей оси со скоростью один оборот в секунду – тогда скорость движения солнечного зайчика по Земле оказывается огромной, и вспышки должны быть слишком короткими. Разработчики признались, что у них и правда есть проблемы с расчетом отражений, и предложили желающим помочь им с этим.

Первым аналитический расчет предложил Quiensabe. (Его пост с не очень оптимистичными результатами был опубликован на ГТ, но его удалили из-за неаккуратной ссылки на краудфандинг — впрочем, это уже другая история.) Такой расчет давал адекватную оценку, но не мог учесть массы факторов – кривизны поверхности Земли, изменение яркости вспышки по мере вращения зеркала и прочего. Нужна была аккуратная численная модель, которая считает яркость отражения с небольшим шагом по времени в зависимости от взаимного расположения спутника, наблюдателя, Солнца и всего прочего.

Задача казалась не слишком сложной – просто нужно было аккуратно разобраться во всей геометрии и ничего не напутать. Мотивировала и возможность посчитать что-нибудь эдакое космическое. Ну и, само собой, радость от написания собственного велосипеда – бесценна.

(Дальше будет много математики. Результаты начинаются здесь, ну а совсем tl;dr — тут.)

Алгоритм

Основная задача модели – рассчитать максимальную яркость вспышки для заданных параметрах спутника и заданного положения наблюдателя на Земле. При этом предсказание вспышек на несколько дней уходит на второй план – для этого есть, например, heavens-above.com. Логично будет взять один пролет спутника над наблюдателем, и уже для него вычислить яркость отражения с небольшим шагом по времени. После этого, варьируя начальные параметры, можно будет вычислить максимально возможную яркость вспышки.

Итак, шаг первый: чтобы увидеть вспышку, необходимо, чтобы спутник: а) пролетел над горизонтом и б) находился не в тени Земли. Шаг второй: если эти условия выполнены, можно начинать считать яркость вспышки. Да, очень вероятно, что спутник будет вращаться – в этом случае вспышки могут быть периодическими. (Мигающий спутник!) Поэтому в общем случае нужен шаг третий – искать отдельные вспышки, разделенные интервалами темноты. Выстраивается логика работы программы, которую реализуют три функции:

• Первая функция считает координаты спутника и наблюдателя в течение одного витка с большим шагом по времени (1-5 секунд) и проверяет, в какие моменты спутник виден над горизонтом и не находится в тени. Если спутник появляется над горизонтом, функция рассчитывает интервалы времени, когда он виден, и передает их второй функции.
• Вторая функция работает с маленьким шагом по времени (0.001 – 1 с) и на каждом шаге рассчитывает яркость вспышки, которую видит наблюдатель.
• Третья функция – самая некосмическая. Она находит отдельные вспышки в результатах работы второй функции и считает их параметры: яркость, длительность, момент прохождения максимума и так далее.

Физическая модель

Моделировать задачу будем в системе отсчета ECI (Earth-centered inertial) – ее начало совпадает с центром Земли, ось Z направлена на северный полюс, а Солнце находится в плоскости X-Z.

Почему именно ECI? Потому, что Солнце в ней неподвижно, наблюдатель вращается вместе с Землей вокруг оси Z, а орбита спутника не меняет своего положения в пространстве. Это значит, что координаты и спутника, и наблюдателя легко параметризуются через углы. Подробности – под спойлером.

Больше координат богу координат!

Координаты наблюдателя задаются широтой fi_obs и текущей долготой theta_obs.

По сути theta_obs является местным временем: полдень соответствует theta_obs = 0, а полночь – theta_obs = 180. Еще нужно знать местное время в начальный момент времени theta_obs0 и скорость вращения Земли вокруг своей оси. Тогда координаты наблюдателя в ECI выражаются так:

Орбиту спутника для облегчения расчетов можно считать круговой. В этом случае ее можно задать тремя параметрами: радиусом R_orb, наклонением fi и долготой восходящего узла theta_AN. Восходящий узел (ascending node) – это точка, в которой орбита пересекает экватор по направлению к северу.

Гравитация, бессердечная ты с. Так как орбита круговая, то спутник движется по ней с неизменной угловой скоростью

Теперь, когда мы знаем все это, положение спутника на орбите легко параметризовать через угол alpha:

Здесь dir_AN и dir_B – направления на восходящий узел и самую высокую точку орбиты, которые считаются простой тригонометрией:

Солнце. Фишка системы отсчета ECI в том, что Солнце всегда находится в плоскости X-Z. В зависимости от времени года оно меняет положение: когда в северном полушарии лето, Солнце находится выше оси Z, а зимой – наоборот. Нас интересует не само положение Солнца, а направление на него. Вернее, от него, то есть вдоль направления движения солнечных лучей. На картинках он называется k_sun.

Посчитали координаты? Теперь надо проверить, виден ли спутник над горизонтом – иначе считать что-либо бессмысленно. Это определяется сонаправленностью векторов от наблюдателя на спутник -dir_{sat_obs} и в зенит dir_obs: если их скалярное произведение положительно, то спутник где-то над нами.

Если спутник над горизонтом, то нужно узнать его координаты на небе. Для этого направление на спутник надо спроецировать на плоскость наблюдателя. Плоскость наблюдателя задается векторами на север dir_nord и запад dir_west. Первый находится через направление на зенит и северный полюс Земли d_NP, а уже по нему и направлению в зенит вычисляется вектор на запад.

Теперь нужно проверить, не находится ли спутник в тени Земли (иначе отражения, очевидно, не будет). Из картинки ниже легко понять, что спутник находится в тени, если одновременно выполняются два условия:

• угол между k_sun и dir_sat (вектор на спутник) острый;
• d < R

Аналогичным способом надо проверить, что наблюдатель находится в тени, то есть у него сейчас ночь. На солнечном свету спутник увидеть будет непросто, хотя в дневной астрономии, несомненно, тоже есть свой шарм.

Наконец, если спутник освещен Солнцем и находится где-то над нами, остается посчитать яркость вспышки. Точнее, ее видимую звездную величину, которая показывает яркость по сравнению с другими звездами. Шкала звездных величин логарифмическая: она определяется как

то есть изменение яркости в 100 раз соответствует изменению звездной величины на 5 единиц. А еще отсюда видно, что шкала звездных величин инвертирована: более ярким объектам соответствуют меньшие значения. Например, звездная величина Арктура и Веги – ярчайших звезд Северного полушария – составляет +0.0, Сириуса – -1.7, а вот Венера достигает -4.6 звездной величины. Глаз различает звезды вплоть до +6 звездной величины, хотя в крупных городах засветка не дает увидеть что-то ярче +4.

Еще немного фотометрии

Для расчета звездной величины нужно знать яркость, а вернее освещенность – то есть мощность светового потока, проходящая через единицу площади. Измеряется она в Вт/м² или в люксах, то есть физический смысл освещенности в том, насколько ярко будет освещена поверхность. Из определения следует, что освещенность не зависит от площади. Поэтому наблюдение невооруженным глазом, фотоаппаратом или в телескоп дадут один и тот же результат измерения звездной величины.

Как посчитать освещенность, которую создает отражение от спутника? Прежде всего, мощность отраженного света (то есть сколько фотонов за секунду отражаются от спутника) находится перемножением площади зеркала на освещенность, создаваемую Солнцем. Нужно не забыть, что зеркало может быть под углом к солнечным лучам. А еще про то, что зеркало неоднородно и отражает свет в большой телесный угол – «зайчик» на поверхности Земли будет протяженным, и мы можем быть вдали от его центра. Наконец, освещенность квадратично зависит от расстояния до спутника – чем дальше спутник, тем меньше света достигнет наблюдателя. Конечное выражение для освещенности выглядит так:

где E_orbit – освещенность, создаваемая Солнцем на орбите Земли, второй сомножитель – площадь зеркала под углом к солнечным лучам, третий – доля телесного угла, который видит наблюдатель, четвертый – квадратичная зависимость от расстояния до спутника.

Сравнивая это число освещенностью, соответствующей нулевой звездной величине (2,5•10^-8 Вт/м²), находим звездную величину спутника.

Увидим мы вспышку или нет, определяется положением зеркала. Его удобно задавать вектором нормали n_mir: это позволяет легко посчитать направление отраженного луча по старому-доброму правилу про угол падения, равный углу отражения:

Легко и просто? Отнюдь! На практике быстрее считать, как должно располагаться зеркало, чтобы при заданном положении спутника и наблюдателя отражение падало бы прямо на наблюдателя. И затем сравнивать с текущим положением зеркала. Преимущество такого подхода в том, что он позволяет изменять шаг по времени. Пока зеркало далеко от оптимального положения, шаг можно увеличить – вспышки в ближайшее время все равно не будет. Если же вспышка уже началась, то шаг по времени нужно уменьшить, чтобы посчитать ее с максимальным разрешением.

Наконец, спутник вращается, а вместе с ним вращается и зеркало. Поэтому нормаль к зеркалу n_mir удобно разложить на две компоненты – параллельную и перпендикулярную оси вращения спутника d_mir. Первая остается неподвижной, вторая прецессирует – математика этого точь-в-точь такая же, как и движение спутника по орбите.

Код

Код написан на матлабе. Как я рассказал выше, он разбит на три функции, которые по очереди вызываются основным скриптом 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

Наконец, если зеркало близко к оптимальной позиции, то можно посчитать и освещенность, которую создает солнечный зайчик. Ну и звездную величину заодно

                    % 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 по местному времени. Пленочные зеркала Маяка отражают свет гораздо хуже антенн Иридиума, зато заметно превосходят их по площади.

После выхода на орбиту спутник должен раскрыть зеркала и закрутиться со скоростью около 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 звездной величины:

Резюмируя:

• Яркость Маяка в районе полуночи будет сопоставима с Вегой. "Ярчайшей звездой на небе" быть не получится — тягаться с Сириусом, планетами, МКС и Иридиумом шансов, к сожалению, нет.
• Вероятность наблюдения ярких вспышек (при условии, что на дворе полночь, а спутник пролетает где-то над головой) — около 20 %.
• При этом сами вспышки будут периодическими — спутник будет лететь и "мигать" примерно раз в секунду.

Вместо послесловия

Несомненно, код можно было улучшать – добавить GUI, эллиптичность орбиты, количество зеркал, параллельные вычисления на CUDA и много всего прочего. Впрочем, совершенству нет предела. К тому же нынешняя версия успешно решает все поставленные перед нею задачи.

Наконец, хочется сказать спасибо Quiensabe за плодотворные дискуссии в процессе разработки модели. А команде Маяка – спасибо за занимательную задачу и, конечно же, успешного полета!