Реферат Исследование движения центра масс межпланетных космических аппаратов

Работа добавлена на сайт bukvasha.net: 2015-10-28

Поможем написать учебную работу

Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.

Предоплата всего

от 25%

Подписываем

договор

Имя

Выберите тип работы:

Принимаю Политику конфиденциальности

Скидка 25% при заказе до 27.12.2024

1. Оглавление.................................................................................... 1

2. Исследовательская часть..................................................... 3

2.1. Введение........................................................................................... 3

2.2. Краткие сведения об орбите............................................................ 4

2.2.1. Характеристика орбиты................................................................ 4

2.2.2. Связь МКА с наземными пунктами управления.......................... 5

2.2.3. Выведение на рабочую орбиту..................................................... 6

2.3. Исходные данные и цели работы................................................... 10

2.3.1. Исходные данные........................................................................ 10

2.3.2. Цели работы................................................................................ 12

2.4. Моделирование движения центра масс МКА................................ 13

2.4.1. Уравнения движения МКА.......................................................... 13

2.4.2. Возмущающие ускорения, действующие на МКА...................... 15

2.4.3. Расчет параметров текущей орбиты МКА.................................. 22

2.5. Коррекция траектории МКА.......................................................... 24

2.5.1. Коррекция приведения................................................................ 24

2.5.2. Расчет потребного топлива......................................................... 26

2.5.3. Коррекция поддержания............................................................. 27

2.6. Движение МКА относительно центра масс................................... 28

2.6.1. Уравнения движения относительно центра масс МКА.............. 28

2.6.2. Стабилизация углового положения при коррекции................... 28

3. ОРГАНИЗАЦИОННО-ЭКОНОМИЧЕСКАЯ ЧАСТЬ....................... 31

3.1. Организация и планирование выполнения темы........................... 31

3.2. Определение затрат труда.............................................................. 31

3.3. Расчет сметы затрат на разработку программного продукта....... 35

4. ПроМЫШЛЕННАЯ ЭКОЛОГИЯ И БЕЗОПАСНОСТЬ.................. 39

4.1. Введение......................................................................................... 39

4.2. Анализ вредных факторов............................................................. 39

4.3. Требования к видеотерминальным устройствам.......................... 44

4.4. Расчет вредных излучений............................................................. 46

4.5. Рациональная организация рабочего места................................... 46

4.6. Рекомендации по снижению утомляемости.................................. 47

4.7. Защита от напряжения прикосновения. Зануление....................... 48

4.8. Пожарная безопасность................................................................. 49

5. Список литературы................................................................. 53

6. Приложение. ТекстЫ Программ для Borland C++ и Matlab 4.0 for windows............................................................................................ 54

6.1. Основной программный модуль main.cpp.................................... 54

6.2. Подпрограмма расчета возмущающих ускорений, параметров орбиты и коррекции sfun.cpp................................................................................................. 57

6.3. Файл начальной инициализации init.h............................................ 77

6.4. Файл описания переменных def.h.................................................. 79

6.5. Файл sfun.h...................................................................................... 81

6.6. Файл rk5.h....................................................................................... 81

6.7. Программа построения временных диаграмм............................... 82

2. ИССЛЕДОВАТЕЛЬСКАЯ ЧАСТЬ
2.1. ВВЕДЕНИЕ
В данной работе проводится исследование движения центра масс МКА под действием различных возмущающих ускорений (от нецентральности гравитационного поля Земли, сопротивления атмосферы, притяжения Солнца и Луны, из-за давления солнечных лучей) и создание математической модели движения ЦМ МКА, позволяющей учесть при интегрировании уравнений движения ЦМ МКА эволюцию орбиты МКА.

В работе разрабатывается алгоритм коррекции, ликвидирующий ошибки выведения МКА и рассчитывается масса топлива, необходимая для проведения коррекции, необходимой из-за эволюции параметров орбиты и из-за ошибок выведения МКА на рабочую орбиту.

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

В работе приводятся программы, реализующие интегрирование уравнений движения ЦМ МКА, процесс осуществления коррекции и расчет топлива для коррекции.

2.2. КРАТКИЕ СВЕДЕНИЯ ОБ ОРБИТЕ
Основными показателями эффективности космической группировки, являются:

- предельная производительность МКА в сутки на освещенной стороне Земли не менее 400-500 объектов.

- периодичность наблюдения районов съемки не реже одного раза в сутки.

Расположение плоскости орбиты по отношению к Солнцу выбрано таким образом, чтобы угол между линией узлов и следом терминатора на плоскости экватора Земли составлял D_т = 30°. При этом северный полувиток орбиты должен проходить над освещенной частью земной поверхности. Для определенности углу D_тприписывается знак «+» в том случае, если восходящий узел орбиты находится над освещенной частью Земли, и знак «-», если ВУ находится над неосвещенной частью. При выборе баллистического построения оперируют углом D, однозначно определяющимся прямым восхождением Солнца a₀ и долготой восходящего узла орбиты в абсолютном пространстве W: D = a₀ - W. Соотношение между углом D_т и углом D: Dº D_т- 90°.
2.2.1. ХАРАКТЕРИСТИКА ОРБИТЫ
Для решения задач наблюдения Земли из космоса с хорошим разрешением при жестких ограничениях на массу КА и минимизации затрат на выведение целесообразно использовать низкие круговые орбиты. В этом классе орбит выделяют солнечно-синхронные орбиты со следующими свойствами:

- скорость прецессии плоскости орбиты в пространстве составляет примерно 1° в сутки, что практически обеспечивает постоянство ориентации ее относительно терминатора Земли в течении всего срока активного существования КА.

- близость наклонения плоскости орбиты к полярному, что обеспечивает глобальность накрытия полюсами обзора поверхности Земли.

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

Всем этим условиям удовлетворяют солнечно-синхронные орбиты с высотами от нескольких сот до полутора тысяч километров. На больших высотах наклонение солнечно-синхронной орбиты отличается от полярного, и глобальность накрытия поверхности Земли не обеспечивается. Для повышения эффективности наблюдения целесообразно выбрать орбиты с изомаршрутной трассой, у которых следы орбит ежесуточно проходят на одними и теми же районами Земли, что позволяет обеспечивать периодичность наблюдения одного и того же объекта, как минимум, раз в сутки с одного КА.

Предварительные расчеты показали, что целесообразно использовать орбиту с высотой Н = 574 км и наклонением плоскости орбиты к плоскости экватора Земли i = 97,6°.

Масса МКА может составить от 500 до 800 кг (что зависит от вида целевой аппаратуры, устанавливаемой на борту МКА). Для выведения МКА на орбиту используется РН СС-19 («Рокот») с разгонным блоком «Бриз».
2.2.3. СВЯЗЬ МКА С НАЗЕМНЫМИ ПУНКТАМИ УПРАВЛЕНИЯ
Управление МКА осуществляется с наземных пунктов управления на территории России. Их количество и место расположения выбирается таким образом, чтобы на любом витке можно было организовать сеанс связи с МКА хотя бы с одного пункта управления. Угол возывшения МКА над горизонтом наземного пункта управления должен быть не менее 7°, а дальность до МКА не должна превышать 2200 км.

В расчете зон связи были использованы следующие исходные данные:

- высота орбиты - 574 км.

- наклонение орбиты - 97,6°.

- географическая долгота восходящего узла первого витка - 4° в.д.

- минимальный угол возвышения МКА над местным горизонтом - 7°.

Из рассматривавшихся возможных наземных пунктов управления (Москва, Новосибирск, Хабаровск, Мурманск, Калининград, Диксон, Комсомольск-на-Амуре, Петропавловск-Камчатский), было выбрано три (Москва, Диксон, Петропавловск-Камчатский), обеспечивающие возможности связи с МКА на любом витке орбиты. При этом зоны связи с МКА составляют от 3 до 9 минут на витке.

Интергральные характеристики возможности связи с МКА:

- высота орбиты - 574 км.

- число витков, видимых из Москвы, вит/сутки - 6.

- суммарное время видимости из Москвы, мин - 41.

- суммарное время видимости с трех пунктов, мин - 153.

- максимальное время видимости одного витка, мин - 9,1.
2.2.4. ВЫВЕДЕНИЕ МКА НА РАБОЧУЮ ОРБИТУ
Выведение МКА на орбиту с наклонением i = 97,6° и высотой Н = 574 км осуществляется ракетой-носителем «Рокот» с разгонным блоком «Бриз». При выведении для каждой отделяющейся части РН (отработанная первая ступень, обтекатель, отработанная вторая ступень) существует свой район падения.

Возможные варианты старта:

1. Полигон Байконур.

Из-за отсутствия зон падения отделяющихся частей возможно сформировать опорную орбиту с наклонением i порядка 65°. Для формирования опорной орбиты с наклонением близким полярному при использовании трассы с азимутом стрельбы более 180° (направление стрельбы на юг) - первая ступень падает в районе Ашхабада, обтекатель сбрасывается на высоте Н порядка 100 км, вторая ступень падает за Аравийским полуостровом. С точки зрения энергетики, выведение осуществляется не по оптимальной схеме, в результате чего на круговую орбиту высотой Н порядка 700 км выводится МКА массой менее 600 кг.

2. Полигон Ледяная (Свободный).

Из-за отсутствия зон падения отделяющихся частей возможно сформировать опорную орбиту с наклонением i порядка 54° и 65°. При северном запуске РН первая ступень падает в районе заповедника в устье реки Олейма (приток Лены).

3. Космодром Плесецк.

Азимуты пуска с космодрома Плесецк обеспечивают наклонения орбит i от 72° до 93°. Формирование требуемового наклонения i = 97,6° осуществляется с помощью разгонного блока «Бриз».

В результате работы двух ступеней РН формируется баллистическая траектория с наклонением i = 93°. Высота в момент окончания работы двигателя второй ступени составляет Н = 190 км, наклонная дальность L = 300 км. Приблизительно через 1,2 секунды после прохождения команды на выключение двигателя второй ступени проходит команда на запуск ДУ РБ. После выключения двигателя второй ступени РН происходит отделение от ракеты связки РБ с КА. Время расцепки t = 318 секунд. Абсолютная скорость в момент отделения V = 5550 м/с. Отделяемая масса 6700 кг.

Двигательная установка РБ «Бриз» выполняет задачу доразгона КА при формировании опорной орбиты.

Характеристики двигателя РБ «Бриз»:

- тяга R, кг - 2000.

- удельный импульс R_уд, сек - 324.

- количество включений, р - 7.

- интервал между включениями, сек - 20.

- время функционирования, час - 7.

В результате работы двигателя РБ «Бриз» при первом включении происходит увеличение высоты баллистической траектории с Н = 190 км до Н = 270 км и к моменту окончания работы двигателя (t = 905,5 сек) в точке с аргументом широты u = 104,1° формируется опорная эллиптическая орбита с параметрами:

- высота в перигее Н_п, км - 190.

- высота в апогее Н_а, км - 574.

- большая полуось орбиты а, км - 6747.

- эксцентриситет e - 0,02548

- наклонение i, ° - 93,4.

- период обращения Т, час - 1,53.

- аргумент перигея w, ° - 128,38.

- долгота восходящего узла в гринвичской СК, фиксированной на момент старта W_г, ° - 48,37.

Величина импульса характеристической скорости, отрабатываемого при первом включении ДУ РБ dV₁ = 2,36 км/с, время работы порядка 600 сек.

Работа двигателя при первом включении происходит вне зоны видимости НПУ на территории России. Географические координаты, соответствующие этому моменту:

- широта j » 76°.

- долгота l » 238°.

В момент прохождения МКА перигея опорной эллиптической орбиты (t = 1231 сек) географические координаты составляют:

- широта j » 53°.

- долгота l » 227°.

На опорной эллиптической орбите МКА совершает пассивный полет до апогея. В районе апогея (t = 1,12 час) осуществляется второе включение ДУ РБ.

В результате приложения второго компланарного импульса характеристической скорости dV₂ = 0,12 км/с, при втором включении (время работы 20 сек) формируется круговая орбита с параметрами:

- высота Н, км - 574.

- наклонение i, ° - 93,4.

- период обращения Т, час - 1,6.

Работа двигателя при втором включении происходит вне зоны видимости НПУ на территории России. Географические координаты, соответствующие этому моменту:

- широта j » 1,5°.

- долгота l » 35,8°.

Для создания круговой, солнечно-синхронной орбиты необходимо изменить наклонение до i = 97,6°. С этой целью осуществляется третье включение ДУ РБ при первом прохождении восходящего узла орбиты (t = 1,3 час).

В результате приложения ортогонального импульса характеристической скорости dV₃ = 0,62 км/с, при третьем включении (время работы 90 сек) формируется солнечно-синхронная круговая орбита с параметрами:

- высота Н, км - 574.

- наклонение i, ° - 97,6.

- период обращения Т, час - 1,6.

- число оборотов в сутки N - 15.

Работа двигателя при третьем включении происходит вне зоны видимости НПУ на территории России. Географические координаты, соответствующие этому моменту:

- широта j » 0°.

- долгота l » 28,1°.

После выключения двигателя при третьем запуске происходит отделение МКА от РБ «Бриз».

Кинематические параметры в гринвичской СК, фиксированной на момент старта РН и оскулирующие элементы орбиты на момент отделения от РБ:

Параметр	Значение
t, сек	4946,5
X, м	4638800
Y, м	5120280
Z, м	689680
V_x, м/с	241,23
V_y, м/с	-1233
V_z, м/с	7473,5
l, °	28,1
T, c	5761,67
e	0,0009
i, °	97,595
R_a, м	6940000
R_п, м	6952000

2.3. ИСХОДНЫЕ ДАННЫЕ И ЦЕЛИ РАБОТЫ
2.3.1. ИСХОДНЫЕ ДАННЫЕ
Номинальная орбита, необходимая для выполнения задач МКА, имеет следующие параметры:

- круговая, e = 0.

- солнечно-синхронная, скорость прецессии линии узлов орбиты W равна скорости обращения Солнца относительно Земли

W = 2p / 365,2422 = 0,0172 рад/сут = 0,98 °/сут.

- изомаршрутная, за сутки МКА совершает целое количество оборотов (n = 15).

Это обеспечивает прохождение МКА над одними и теми же районами в одно и тоже местное время.

- период Т = 5765 с.

- высота орбиты Н = 574 км.

- наклонение орбиты i = 97,6°.

- географическая долгота восходящего узла орбиты l_э= 28,1°.

Долгота восходящего узла в геоцентрической экваториальной (абсолютной) системе координат OXYZ определяется как разность

l_э- s₀,

где s₀- часовой угол, отсчитывающийся от гринвичского меридиана до оси X, направленной в точку весеннего равноденствия.

Часовой угол зависит от даты старта и выбирается из астрономического ежегодника. В данной задаче для моделирования выбран часовой угол = 0.

Следовательно долгота восходящего узла орбиты W = l_э= 28,1°.

Исходя из ТЗ, начальная точка выведения имеет следующие координаты в гринвичской системе координат, фиксированной на момент старта РН:

Параметр	Значение
t, сек	4946.5
X, м	4638800
Y, м	5120280
Z, м	689506,95
V_x, м/с	241,23
V_y, м/с	-1233
V_z, м/с	7472,65

Элементы орбиты:

l, °	28,1
T, c	5761,67
e	0,0009
i, °	97,595
R_a, м	6940000
R_п, м	6952000

Кинематические параметры в геоцентрической экваториальной системе координат:

t, сек	4946.5
X, м	6137262,9
Y, м	3171846,1
Z, м	689506,95
V_x, м/с	-201,3
V_y, м/с	-1247,03
V_z, м/с	7472,65
l, °	28,1

Точность выведения:

- предельная ошибка по координате (3s) - 7 км.

- предельная ошибка по скорости (3s) - 5 м/с.

Пересчитав ошибку по координате на ошибку по периоду выведения орбиты получим предельную ошибку по периоду DT - 10 сек.

Корреляционная матрица ошибок выведения на момент выведения составляет:

Члены, стоящие на главной диагонали представляют собой квадраты предельных ошибок - (3s)².

K₁₁ = K₂₂ = K₃₃ = (3s)²= 7² = 49 км.

K₄₄ = K₅₅ = K₆₆ = (3s)²= 5² = 25 м/с.

Остальные члены представляют собой вторые смешанные моменты K_ij= K_ji = r_ijs_is_j или K_ij = K_ji = r_jj(3s_i)(3s_j), где r_jj - коэффициенты связи величин i и j. В данном случае вторые смешанные моменты K_ij= K_ji = 0.

Кинематические параметры в геоцентрической экваториальной системе координат на момент выведения с учетом ошибок выведения:

t, сек	4946.5
X, м	6144262,9
Y, м	3178846,1
Z, м	696506,95
V_x, м/с	-206,3
V_y, м/с	-1252,03
V_z, м/с	7477,65
l, °	28,1

Параметры орбиты с учетом ошибок выведения:

l, °	28,13
T, c	5795,7
W, °	28,13
p, км	6973,5
а, км	6973,6
e	0,00314
i, °	97,637

2.3.2. ЦЕЛИ РАБОТЫ
1) Исследование и моделирование движения ЦМ МКА при воздействии на КА возмущающих ускорений.

2) Разработка алгоритмов проведения коррекции траектории МКА, моделирования процесса, и расчет потребного топлива для проведения коррекции траектории.

3) Исследование динамики системы коррекции траектории при стабилизации углового положения в процессе проведения коррекции траектории МКА.

2.4. МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ ЦЕНТРА МАСС МКА
2.4.1.УРАВНЕНИЕ ДВИЖЕНИЯ КА
Рассмотрим невозмущенное движение материальных точек М и m в некоторой инерциальной системе координат. Движение совершается под действием силы притяжения F_z. Сила F_z для материальной точки m определяется формулой:

,

где ¦ - постоянная притяжения,

r_o
- единичный вектор, направленный от М к m,

,

где - радиус-вектор, проведенный из т.М до т.m.

r - относительное расстояние от М до m.

На точку М действует сила F_z, равная по величине и направленная в противоположную сторону.

На основе второго закона Ньютона уравнения движения материальных точек М и m имеют вид:

(1), (2)

или

(3), (4)

где p₁- радиус-вектор, проведенный из начала инерциальной системы координат в точку m.

p₂- радиус-вектор, проведенный из начала инерциальной системы координат в точку М.

.

Вычитая из уравнения (3) уравнение (4), получим уравнение движения материальной точки m относительно притягивающего центра М:

Так как m<<М, следовательно, можно пренебречь ускорением, которое КА с массой m сообщает притягивающему центру М. Тогда можно совместить начало инерциальной системы координат с притягивающим центром М. Следовательно, .

Таким образом, уравнение невозмущенного движения КА относительно притягивающего центра М в инерциальной системе координат, центр которой находится в М, имеет вид

,

где m
= fM - гравитационная постоянная Земли.

Рассмотрим возмущенное движение КА в геоцентрической экваториальной (абсолютной) системе координат OXYZ:

- начало О - в центре масс Земли.

- ось X направлена в точку весеннего равноденствия g.

- ось Z совпадает с осью вращения Земли и направлена на Северный полюс Земли.

- ось Y дополняет систему до правой.

Движение КА в абсолютной системе координат OXYZ происходит под действием центральной силы притяжения Земли F_z, а также под действием возмущающих сил F_в. Уравнение движения имеет вид

или

где m = 597 кг - масса КА.

В проекциях на оси абсолютной системы координат OXYZ получим

или

или

или

где ax
_в
, ay
_в
, az
_в- возмущающиеся ускорения.

Основные возмущающиеся ускорения вызываются следующими причинами:

- нецентральностью поля притяжения Земли.

- сопротивлением атмосферы Земли.

- влиянием Солнца.

- влиянием Луны.

- давлением солнечного света.
2.4.2. ВОЗМУЩАЮЩИЕ УСКОРЕНИЯ, ДЕЙСТВУЮЩИЕ НА МКА
1) Возмущающееся ускорение, вызванное нецентральностью гравитационного поля Земли.

Рассмотрим потенциал поля притяжения Земли. При точном расчете параметров орбиты спутников, в качестве хорошего приближения к действительной поверхности Земли принимают геоид. Геоид - это гипотетическая уровенная поверхность, совпадающая с поверхностью спокойного океана и продолженная под материком.

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

Потенциал притяжения Земли можно представить в виде разложения по сферическим функциям.

где m
_z
= fM
_z - гравитационная постоянная Земли.

r₀
- средний экваториальный радиус Земли.

с
_nm
, d_nm - коэффициенты, определяемые из гравиметрических данных, а также по наблюдениям за движением ИСЗ.

L - долгота притягивающей точки.

j - широта притягивающей точки.

P_nm(sin
j
) - присоединенные функции Лежандра степени m и порядка n (при m ¹ 0).

P_nm(sin
j
) - многочлен Лежандра порядка n (при m = 0).

Составляющие типа (
m
_z
/r)(r₀/r)ⁿc_n0P_n0(sin
j
) - называют зональными гармониками n-порядка. Т.к. полином Лежандра n-го порядка имеет n действительных корней, функция P_n0(sin
j
) будет менять знак на n широтах, сфера делится на n+1 широтную зону, где эти составляющие имеют попеременно «+» или «-» значения. Поэтому их называют зональными гармониками.

Составляющие типа

(
m
_z
/r)(r₀/r)ⁿc_nmcos(mL)P_nm(sin
j
)
и
(
m
_z
/r)(r₀/r)ⁿd_nmsin(mL)P_nm(sin
j
)

- называют тессеральными гармониками n-порядка и степени m. Они обращаются в 0 на 2m меридианах, где cos(mL) = 0 и sin(mL)
= 0 и на n-m параллелях, где P_nm(sin
j
)

= 0 или d^mP_nm(sin
j
)/d(sin
j
)^m= 0
, сфера делится на n+m+1 трапецию, где эти составляющие сохраняют знак.

Составляющие типа и

(
m
_z
/r)(r₀/r)ⁿc_nncos(nL)P_nn(sin
j
)
и
(
m
_z
/r)(r₀/r)ⁿd_nnsin(nL)P_nn(sin
j
)

- называют секториальными гармониками n-порядка и степени m. Эти составляющие меняю знак только на меридианах, cos(nL) = 0 и sin(nL)
= 0, на сфере выделяют 2n меридиональных секторов, где эти составляющие сохраняют знак.

Многочлен Лежандра степени n находится по следующей формуле:

P_n0(z) = 1/(2ⁿn!)
´
(dⁿ(z² - 1)ⁿ/dzⁿ)

Присоединенная функция Лежандра порядка n и степени m находится по следующей формуле:

P_nm(z) = (1-z²)^m/2
´
d^mP_n0(z)/dz^m

Возмущающая часть гравитационного потенциала Земли равна

U
_в

= U’ +
D
U’ = (U -
m
_z
/r) +
D
U’

где D
U’ - потенциал аномалий силы тяготения Земли.

U’ - часть потенциала Земли, которая учитывает несферичность Земли.

Следовательно,

Первая зональная гармоника в разложении потенциала учитывает полярное сжатие Земли.

Зональные гармоники нечетного порядка и тессеральные гармоники, где n-m нечетное число - учитывают ассиметрию Земли относительно плоскости экватора.

Секториальные и тессеральные гармоники - учитывают ассиметрию Земли относительно оси вращения.

Первая зональная гармоника имеет порядок 10^-3, а все остальные - порядок 10^-6и выше. Поэтому будем учитывать в разложении потенциала притяжения только зональную гармонику (n=2, m=0) и секторальную гармонику (n=2, m=2). Также не будем учитывать потенциал аномалий силы тяготения Земли DU’.

Таким образом,

U
_в
= (
m
_z
/r)(r₀/r)²[c₂₀P₂₀(sin
j
) + (c₂₂cos(2L) + d₂₂sin(2L))P₂₂(sin
j
)],

где c₂₀= - 0,00109808,

c₂₂ = 0,00000574,

d₂₂ = - 0,00000158.

P₂₀(x) = 1/2²2!
´
d²(x² - 1)²/dx².

Следовательно P₂₀(x) = (3x² - 1)/2.

Так как sinj
= z/r, следовательно P₂₀(sinj
) = (3(z/r)² - 1)/2.

P₂₂(x) = (1 - x²)^2/2
´
d²P₂₀(x)/dx²= 1/2
´
(1 - x²)
´
d²(3x² - 1)/dx²

Следовательно P₂₂(x) = 3(1 - x²).

Так как sinj
= z/r, следовательно P₂₂(sinj
) = 3(1 - (z/r)²).

Значит

Чтобы найти возмущающее ускорение от нецентральности поля тяготения Земли в проекциях на оси абсолютной системы координат OXYZ, надо взять производные от возмущающего потенциала U
_в по координатам X, Y, Z, причем r =
Ö
(x² + y² + z²).

Следовательно,

2) Возмущающее ускорение, вызванное сопротивлением атмосферы.

При движении в атмосфере на КА действует сила аэродинамического ускорения R_x, направленная против вектора скорости КА относительно атмосферы:

где C_x
= 2 - коэффициент аэродинамического сопротивления.

S
_м = 2,5 м² - площадь миделевого сечения - проекция КА на плоскость, перпендикулярную направлению скорости полета.

V - скорость КА.

r - плотность атмосферы в рассматриваемой точке орбиты.

Так как исследуемая орбита - круговая с высотой Н = 574 км, будем считать, что плотность атмосферы одинакова во всех точках орбиты и равна плотности атмосферы на высоте 574 км. Из таблицы стандартной атмосферы находим плотность наиболее близкую к высоте Н = 574 км. Для высоты Н = 580 км r = 5,098´10^-13 кг/м³.

Сила аэродинамического ускорения создает возмущающее касательное ускорение a_a:

Найдем проекции аэродинамического ускорения на оси абсолютной системы координат a_xa, a_ya, a_za:

a_a направлено против скорости КА, следовательно единичный вектор направления имеет вид

e_a= [V_x/|V|, V_y|V|, V_z/|V|], |V| =
Ö
(V_x²+V_y²+V_z²)

Таким образом,

Значит

, ,
3) Возмущающее ускорение, вызванное давлением солнечного света.

Давление солнечного света учитывается как добавок к постоянной тяготения Солнца - D
m
_c
. Эта величина вычисляется следующим образом:

D
m
_c
= pS
_м
A²/m

где p = 4,64´10^-6 Н/м² - давление солнечного света на расстоянии в одну астрономическую единицу А.

A = 1,496´10¹¹ м - 1 астрономическая единица.

m - масса КА.

S
_м
= 8 м² - площадь миделевого сечения - проекция КА на плоскость, перпендикулярную направления солнечных лучей.

Таким образом,

D
m
_c= 1,39154´10¹⁵ м³/c².
4) Возмущающее ускорение, возникающее из-за влияния Солнца.

Уравнение движения КА в абсолютной системе координат OXYZ относительно Земли при воздействии Солнца:

где m
_z - постоянная тяготения Земли.

m
_c
- постоянная тяготения Солнца.

r - радиус-вектор от Земли до КА.

r_c
- радиус-вектор от Земли до Солнца.

Таким образом, возмущающее ускорение, возникающее из-за влияния Солнца:

.

Здесь первое слагаемое есть ускорение, которое получил бы КА, если он был непритягивающим, а Земля отсутствовала.

Второе слагаемое есть ускорение, которое сообщает Солнце Земле, как непритягивающему телу.

Следовательно, возмущающее ускорение, которое получает КА при движении относительно Земли - это разность двух слагаемых.

Так как r_c>>r, то в первом слагаемом можно пренебречь r. Следовательно

| r_c- r| =
Ö
((x_c-x)²+(y_c-y)²+(z_c-z)²)

где x_c, y_c, z_c- проекции радиуса-вектора Солнца на оси абсолютной системы координат.

Моделирование движения Солнца проводилось следующим образом: за некоторый промежуток времени t Солнце относительно Земли сместится на угол J
=
J
_н

+
w
_c
t
,

где J
_н

=
W
+ (90 -
D
) - начальное положение Солнца в эклиптической системе координат.

W = 28,1° - долгота восходящего узла первого витка КА.

D = 30° - угол между восходящим узлом орбиты КА и терминатором.

w
_c - угловая скорость Солнца относительно Земли.

w
_c
= 2
p
/T = 2p/365,2422´24´3600 = 1,991´10^-7 рад/c = 1,14´10^-5°/c

Таким образом, в эклиптической системе координат проекции составляют:

x_c^e = r_ccos
J

y_c^e = r_csin
J

z_c^e = 0

r_c= 1,496´10¹¹ м (1 астрономическая единица) - расстояние от Земли до Солнца

Плоскость эклиптики наклонена к плоскости экватора на угол e
= 23,45°, проекции r_c на оси абсолютной системы координат можно найти как

x_c = x_c^e = r_ccos
J

y_c^e = y_c^ecos
e
= r_ccos
J
cos
e

z_c^e = r_csin
J
sin
e

Таким образом, проекции возмущающего ускорения на оси абсолютной системы координат:

a_xc= -
m
_c
x/(
Ö
((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_yc= -
m
_c
y/(
Ö
((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_zc= -
m
_c
z/(
Ö
((x_c-x)²+(y_c-y)²+(z_c-z)²))³

С учетом солнечного давления

a_xc= - (
m
_c
-
D
m
_c
)x/(
Ö
((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_yc= - (
m
_c
-
D
m
_c
)y/(
Ö
((x_c-x)²+(y_c-y)²+(z_c-z)²))³

a_zc= - (
m
_c
-
D
m
_c
)z/(
Ö
((x_c-x)²+(y_c-y)²+(z_c-z)²))³
5) Возмущающее ускорение, возникающее из-за влияния Луны.

Уравнение движения КА в абсолютной системе координат OXYZ относительно Земли при воздействии Луны:

где m
_л = 4,902´10⁶ м³/c²- постоянная тяготения Луны.

r
_л - радиус-вектор от Земли до Луны.

Таким образом, возмущающее ускорение, возникающее из-за влияния Луны:

Так как r_л>>r, то в первом слагаемом можно пренебречь r. Следовательно

|r
_л

- r| =
Ö
((x
_л
-x)²+(y
_л
-y)²+(z
_л
-z)²)

где x
_л
, y
_л
, z
_л- проекции радиуса-вектора Луны на оси абсолютной системы координат.

Движение Луны учитывается следующим образом: положение Луны в каждый момент времени рассчитывается в соответствии с данными астрономического ежегодника. Все данные заносятся в массив, и далее этот массив считается программой моделирования движения КА. В первом приближении принимается:

- орбита Луны - круговая.

- угол наклона плоскости орбиты Луны к плоскости эклиптики i = 5,15°.

- период обращения линии пересечения плоскостей лунной орбиты и эклиптики (по ходу часовой стрелки, если смотреть с северного полюса) = 18,6 года.

Угол между плоскостями экватора Земли и орбиты Луны можно найти по формуле

cos(
h
_л
) = cos(
e
)cos(i) - sin(
e
)sin(i)cos(
W
_л
)

где W
_л - долгота восходящего узла лунной орбиты, отсчитывается от направления на точку весеннего равноденствия.

e
- угол между плоскостями эклиптики и экватора Земли.

Величина h_л колеблется с периодом 18,6 лет между минимумом при h
_л
=
e
- i = 18°18’ и максимумом при
h
_л
=
e
+ i = 28°36’ при W
= 0.

Долгота восходящего узла лунной орбиты W_л изменяется с течением времени t на величину W
_л= t´360/18,6´365,2422´24´3600.

Положение Луны на орбите во время t определяется углом

J
_л= t´360/27,32´24´3600.

По формулам перехода найдем проекции вектора положения Луны на оси абсолютной системы координат:

x
_л
= r
_л
(cos
J
_л
cos
W
_л

- cos
h
_л
sin
J
_л
sin
W
_л
)

y
_л
= r
_л
(cos
J
_л
sin
W
_л

+ cos
h
_л
sin
J
_л
cos
W
_л
)

z
_л
= r
_л
sin
h
_л
sin
J
_л

r
_л
= 3,844´10⁸ м - среднее расстояние от Земли до Луны

Таким образом, проекции возмущающего ускорения на оси абсолютной системы координат:

a_x
_л

= -
m
_л
x/(
Ö
((x
_л!
-x)²+(y
_л
-y)²+(z
_л
-z)²))³

a_y
_л

= -
m
_л
y/(
Ö
((x
_л!
-x)²+(y
_л
-y)²+(z
_л
-z)²))³

a_z
_л

= -
m
_л
z/(
Ö
((x
_л!
-x)²+(y
_л
-y)²+(z
_л
-z)²))³

Уравнения возмущенного движения при действии корректирующего ускорения имеют вид:

или

d²x/dt² = - (
m
_z
/r²)x + a_xu+ a_xa+ a_xc+ a_x
_л

+ a_x
_к

d²y/dt² = - (
m
_z
/r²)y + a_yu+ a_ya+ a_yc+ a_y
_л

+ a_y
_к

d²z/dt² = - (
m
_z
/r²)z + a_zu+ a_za+ a_zc+ a_z
_л

+ a_z
_к

2.4.3. РАСЧЕТ ПАРАМЕТРОВ ТЕКУЩЕЙ ОРБИТЫ КА
Полученная система уравнений движения ЦМ КА интегрируется методом Рунге-Кутта 5-го порядка с переменным шагом. Начальные условия x₀, y₀, z₀, V_x0, V_y0, V_z0 - в абсолютной системе координат, соответствуют начальной точке вывода при учете ошибок выведения. После интегрирования мы получаем вектор состояния КА (x, y, z, V_x, V_y, V_z) в любой момент времени.

По вектору состояния можно рассчитать параметры орбиты. соответствующие этому вектору состояния.

а) Фокальный параметр - р.

р
= C²/
m
_z, где С - интеграл площадей.

C = r
´
V, |C| = C =
Ö
(C_x²+C_y²+C_z²)

C_x = yV_z - zV_y

C_y = zV_x - xV_z- проекции на оси абсолютной СК

C_z = xV_y - yV_x

б) Эксцентриситет - е.

e = f/
m
_z
, где f
- вектор Лапласа

f = V
´
C -
m
_z
r/r, |f| = f =
Ö
(f_x²+f_y²+f_z²)

f_x = V_yC_z - V_zC_y -
m
_z
x/r

f_y = V_zC_x - V_xC_z -
m
_z
y/r - проекции на оси абсолютной СК

f_z = V_xC_y - V_yC_x -
m
_z
z/r

в) Большая полуось орбиты.

a = p/(1 - e²)

г) Наклонение орбиты - i.

C_x = Csin(i)sin
W

C_y = - Csin(i)cos
W

C_z = Ccos(i)

можно найти наклонение i = arccos(C_z/C)

д) Долгота восходящего узла -
W
.

Из предыдущей системы можно найти

sin
W
= C_x/Csin(i)

cos
W
= - C_y/Csin(i)

Так как наклонение орбиты изменяется несильно в районе i = 97,6°, мы имеем право делить на sin(i).

Если sinW
=> 0,
W
= arccos (-C_y/Csin(i))

Если sinW
< 0,
W
= 360 - arccos (-C_y/Csin(i))

е) Аргумент перицентра - w.

f_x = f(cos
w
cos
W
- sin
w
sin
W
cos(i))

f_y = f(cos
w
sin
W
+ sin
w
cos
W
cos(i))

f_z = fsin
w
sin(i)

Отсюда найдем

cos
w
= f_xcos
W
/f + f_ysin
W
/f

sin
w
= f_z/fsin(i)

Если sinw
> 0,
w
= arccos (f_xcos
W
/f + f_ysin
W
/f)

Если sinw
< 0,
w
= 360 - arccos (f_xcos
W
/f + f_ysin
W
/f)

ж) Период обращения - Т.

T = 2
p
Ö
(a³/
m
_z
)

Графики изменения элементов орбиты при действии всех, рассмотренных выше, возмущающих ускорений в течение 2-х периодов (Т = 5765 с) приведены на рис. 1-12.

Графики изменения во времени возмущающих ускорений приведены на рис. 13-18.

2.5. ПРОВЕДЕНИЕ КОРРЕКЦИИ ТРАЕКТОРИИ МКА
Существующие ограничения на точки старта РН и зоны падения отработавших ступеней РН, а также ошибки выведения не позволяют сразу же после пуска реализовать рабочую орбиту. Кроме того, эволюция параметров орбит под действием возмущающих ускорений в процессе полета МКА приводит к отклонению параметров орбиты КА от требуемых значений. Для компенсации воздействия указанных факторов осуществляется коррекция орбиты с помощью корректирующей двигательной установки (КДУ), которая располагается на борту МКА.

В данной работе проведена разработка алгоритма коррекции, моделирование процесса коррекции и расчет топлива, необходимого для проведения коррекции.

Из-за различных причин возникновения отклонений элементов орбиты проводится:

- коррекция приведения - ликвидация ошибок выведения и приведение фактической орбиты к номинальной с заданной точностью.

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

Для того, чтобы орбита отвечала заданным требованиям, отклонения параметров задаются следующим образом:

- максимальное отклонение наклонения орбиты Di = 0,1°

- предельное суточное смещение КА по долготе Dl = 0,1°

Следовательно, максимальное отклонение периода орбиты DT = 1,6 сек.

Алгоритм коррекции следующий:

1) Коррекция приведения.

2) Коррекция поддержания.
2.5.1. КОРРЕКЦИЯ ПРИВЕДЕНИЯ
После окончания процесса выведения МКА, проводятся внешне-траекторные измерения (ВТИ). Эти измерения обеспечивают, по баллистическим расчетам, знание вектора состояния с требуемой точностью через 2 суток. После этого начинается коррекция приведения.

Предложена следующая схема проведения коррекции:

а) Коррекция периода.

б) Коррекция наклонения.

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

Осуществляется в два этапа:

- коррекция перицентра

- коррекция апоцентра

Сначала осуществляется коррекция перицентра - приведение текущего расстояния до перицентра r_p к номинальному радиусу r_н = 6952137 м. После измерения вектора состояния рассчитываются параметры орбиты. Далее определяется нужный корректирующий импульс DV_к. Направление импульса (тормозящий или разгоняющий) зависит от взаимного расположения перицентра орбиты и радиуса номинальной орбиты. Для этого вычисляется Dr_p= r_p- r_н.

Возможны ситуации:

1)

Dr_p< 0 - прикладывается разгоняющий импульс

2) Dr_p> 0 - прикладывается тормозящий импульс

КА долетает до апоцентра и в апоцентре прикладывается корректирующий импульс. Время работы КДУ - 20 сек.

Так как время работы КДУ ограничено, а DV_к может быть большим, следовательно, далее рассчитывается максимальный импульс скорости DV_max за 20 сек работы двигателя:

DV_max = Pt/m = 25´20/597 = 0,8375 м/с

Если DV_к > DV_max в апоцентре прикладывается импульс DV_к = DV_max. В результате этого r_pнемного корректируется. На следующем витке опять рассчитывается DV_к, и если на этот раз DV_к < DV_max, в апоцентре прикладывается импульс DV_к. КДУ включается не на полную мощность P = (DV_к/DV_max)P_max.

Время включения = 20 сек.

Это происходит до тех пор, пока не приблизится к r_pс заданной точностью.

После того, как скорректирован перицентр, начинается коррекция апоцентра. Рассчитываются параметры орбиты и нужный корректирующий импульс, такой, чтобы r_a=r_н = 6952137 м. Направление корректирующего импульса также зависит от величин r_aи r_н.

Вычисляется Dr_a= r_a- r_н.

Возможна ситуация:

Dr_a> 0 - в перицентре прикладывается тормозящий импульс.

КА долетает до перицентра и в перицентре прикладывается корректирующий импульс. Время работы КДУ - 20 сек.

Так как время работы КДУ ограничено, а DV_к может быть большим, следовательно, далее рассчитывается максимальный импульс скорости DV_max за 20 сек работы двигателя:

DV_max = Pt/m = 25´20/597 = 0,8375 м/с

Если DV_к > DV_max, в перицентре прикладывается импульс DV_к = DV_max. В результате этого немного корректируется r_a. На следующем витке опять рассчитывается DV_к, и если на этот раз DV_к < DV_max, в перицентре прикладывается импульс DV_к. КДУ включается не на полную мощность P = (DV_к/DV_max)P_max.

Время включения = 20 сек.

Это происходит до тех пор, пока r_a не приблизится к r_нс заданной точностью.

Таким образом осуществляется коррекция перехода.
б) Коррекция наклонения.

После коррекции периода проводятся внешне-траекторные измерения и получают вектор состояния КА. Если снова необходима коррекция периода ее проводят еще раз и снова измеряют вектор состояния КА.

Далее проводится коррекция наклонения по такой же схеме. Коррекция производится в точке пересечения орбиты КА с линией узлов.
После того, как рассчитаны корректирующие импульсы скорости, по формулам перехода проекции вектора на оси абсолютной системы координат. Далее рассчитывается корректирующее ускорение и подставляется в уравнения движения центра масс КА. После этого уравнения интегрируются методом Рунге-Кутта 5-го порядка с переменным шагом.

Графики изменения элементов орбиты в процессе коррекции приведения приведены на рис.19-30.
2.5.2. РАСЧЕТ ПОТРЕБНОГО ТОПЛИВА
Масса топлива, необходимого для проведения коррекции траектории рассчитывается по формуле Циолковского:

m = m₀(1 - e^-^D^V^к^/W)

m₀= 597 кг - начальная масса МКА (кг)

W = 2200 м/с - скорость истечения газов из сопла двигателя.

Результаты проведения коррекции приведения:

	tн, с	tк, с	Dt, с	DV_к, м/c	Имп.	m, кг
Коррекция периода	176242	262592	300	12,1	15	3,26
Коррекция наклонения	273984	432298	580	24,11	29	6,48

2.5.3.КОРРЕКЦИЯ ПОДДЕРЖАНИЯ
Основная задача МКА - проведение съемки определенных районов Земли по крайней мере один раз в сутки, т.е. трасса КА должна проходить над заданным районом каждые сутки.

Требования для проведения коррекции:

- предельное суточное смещение орбиты по долготе Di = 0,1°

- предельное отклонение наклонения Dl = 0,1°.

В пересчете отклонения Dl на отклонение по периоду получим:

DT = 1,597 сек. - максимальное отклонение по периоду.

При помощи программы моделирования было просчитано 3 месяца и получено, что средний период изменился на 3,2 сек, а наклонение - на 0,001°.

Таким образом, коррекцию периода надо делать примерно 1 раз в 1,5 мес.

Нужный импульс скорости - 1 м/с за время активного существования - 5 лет - коррекцию периода надо провести 40 раз, DV = 40 м/с, масса топлива = 10,8 кг.

За 5 лет Di = 0,02° - коррекцию наклонения проводить не надо.

Графики изменения элементов орбиты за 3 месяца приведены на рис.31-42.

2.6. ДВИЖЕНИЕ МКА ОТНОСИТЕЛЬНО ЦЕНТРА МАСС
2.6.1. УРАВНЕНИЯ ДВИЖЕНИЯ ОТНОСИТЕЛЬНО ЦМ КА
При рассмотрении движения относительно ЦМ КА используют уравнения Эйлера:

J_x
w
_x
+ (J_z-J_y)
w
_y
w
_z
= M_xy + M_x
_в

J_y
w
_y
+ (J_x-J_z)
w
_x
w
_z
= M_yy + M_y
_в

J_z
w
_z
+ (J_y-J_x)
w
_y
w
_x
= M_zy + M_z
_в

где J_x
,
J_y
,
J_z - главные моменты инерции,

M_y - управляющий момент,

M
_в - возмущающий момент.

Так как угловые скорости КА малы, следовательно, можно пренебречь произведением угловых скоростей, значит, уравнения Эйлера имеют вид:

J_x
w
_x
= M_xy + M_x
_в

J_y
w
_y
= M_yy + M_y
_в

J_z
w
_z
= M_zy + M_z
_в

Главные моменты инерции:

J_x= 532 кг´м², J_y= 563 кг´м², J_z= 697 кг´м².

Центробежные моменты инерции принимаются равными 0.

Возмущающий момент M_в возникает из-за того, что двигатель коррекции расположен не в центре масс КА, и реактивная тяга, линия действия которой находится на удалении (плече) l от центра масс КА, создает паразитный крутящий момент M_в.

M_в= P´l,

где P = 25 H - тяга корректирующего двигателя,

l = 4 мм - плечо.

Таким образом, M_в= 25´0,0004 = 0,1 Нм.
2.6.2. СТАБИЛИЗАЦИЯ УГЛОВОГО ПОЛОЖЕНИЯ ПРИ КОРРЕКЦИИ
Основное требование, предъявляемое в этом режиме:

- точность поддержания направления импульса коррекции - не хуже 1 угл.мин.

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

Функциональная схема МКА состоит из следующих эелементов:

1) МКА - малый космический аппарат.

МКА описывается как абсолютно твердое тело.

2) ДУС - датчик угловой скорости.

В качестве ДУС используется командный гироскопический прибор. Он описывается колебательным звеном с параметрами T = 1/30 c^-1 и e = 0,7, а также нелинейным звеном с насыщением 2°/сек.

3) АЦП - аналогово-цифровой преобразователь.

Преобразует аналоговый сигнал с ДУС в цифровой сигнал.

4) ЦАП - цифро-аналоговый преобразователь.

Преобразует цифровой сигнал с ЦВМ в аналоговый.

5) ШИМ - широтно-импульсный модулятор.

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

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

6) Двигатель стабилизации.

Двигатель описывается нелинейностью с насыщением 0,127 Нм и звеном запаздывания с Т_д = 0,05 сек.

Тяга двигателя 0,1 Н

7) ЦВМ.

В ЦВМ формируется управление по углу и угловой скорости. Закон управления имеет вид:

e = K(K₁j +K₂j), К = 1, К₁ = 550, К₂ = 430.

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

Система была промоделирована по каналу х. Для других каналов схемы моделирования будут аналогичными.

Для разомкнутой системы были построены ЛАЧХ и ФЧХ. Эти графики представлены на рис.43.

Результаты моделирования замкнутой системы представлены на рис.44-46.

Таким образом, в результате моделирования получено, что процесс стабилизации углового положения происходит примерно за 15 сек., статическая точность поддержания углового положения - 0,62 угл.мин., что полностью удовлетворяет требованиям технического задания.

3. ОРГАНИЗАЦИОННО-ЭКОНОМИЧЕСКАЯ ЧАСТЬ
3.1. ОРГАНИЗАЦИЯ И ПЛАНИРОВАНИЕ ВЫПОЛНЕНИЯ ТЕМЫ
Сроки выполнения и затраты на исследования в большой мере зависят от организационных условий выполнения исследовательских работ Поэтому необходимо в первую очередь определить, хотя бы в общем виде, порядок и организацию проведения дипломной работы по заданной теме.

Организация дипломной работы по любой теме складывается из определённых этапов и подэтапов, каждый из которых хотя и может иметь разное содержание, однако структурно занимает равное положение для всех дипломных работ, выполняемых в данной отрасли.

Таким образом, структура дипломной работы может быть сформирована по типовой схеме, упорядоченной в соответствии с конкретным видом исследования. Состав дипломной работы по заданной теме, а также потребные категории исследований по этапам и подэтапам представлены в табл.1.

Этапы		Содержание	Исполнители
1.	Техническое задание	Постановка задачи. Определение состава программного продукта.	Руководитель Разработчик
2.	Эскизный проект	Разработка общего описания программного продукта.	Руководитель Разработчик
3.	Технический проект	Разработка структуры программного продукта.	Разработчик
4.	Рабочий проект	Программирования и отладка программы. Проверка результатов и внесение корректив в программу.	РуководительРазработчик
5.	Внедрение	Оформление необходимой документации.	Разработчик

3.2. ОПРЕДЕЛЕНИЕ ЗАТРАТ ТРУДА
Первым шагом при определении себестоимости программного комплекса является расчет трудоемкости создания и внедрения. Расчет производится по методике, приведенной в документе «Типовые нормы времени на программирование задач для ЭВМ». Типовые нормы времени предназначены для определения затрат времени на разработку программных средств вычислительной техники (ПСВТ).

Исходными данными для расчета трудоемкости, при разработке программы являются:

Количество разновидностей форм входной информации - 2,

в том числе:

информации, получаемой от решения смежных задач - 1,

справочной, условно постоянной информации (файл инициализации) - 1;

Количество разновидностей форм выходной информации - 2,

в том числе:

печатных документов (временные диаграммы) - 1,

информации, наносимой на магнитные носители (файл инициализации) - 1;

Степень новизны комплекса задач - Г (разработка программной продукции, основанной на привязке типовых проектных решений).

Сложность алгоритма - 3 (реализуются стандартные методы решения, не предусмотрено применение сложных численных и логических методов).

Вид используемой информации:

количество разновидностей форм переменной информации (ПИ) - 1, в том числе: информации, получаемой от решения смежных задач - 1;

количество разновидностей форм нормативно-справочной информации (НСИ) (файл инициализации) - 1;

Язык программирования - Borland С++.

Вид представления исходной информации - группа 11 (требуется учитывать взаимовлияние различных показателей).

Вид представления выходной информации - группа 22 (вывод информационных массивов на машинные носители).

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

t_пп = t_тз + t_эп + t_тп + t_рп + t_в,

где t_тз - трудоемкость разработки технического задания на создание программного продукта,

t_эп - трудоемкость разработки эскизного проекта программного продукта,

t_тп - трудоемкость разработки технического проекта программного продукта,

t_рп - трудоемкость разработки рабочего проекта программного продукта,

t_в - трудоемкость внедрения программного продукта.

Трудоемкость разработки технического задания рассчитывается по формуле

t_тз= Т^з_рз + Т^з_рп,

где Т^з_рз - затраты времени разработчика постановки задач на разработку ТЗ, чел.-дни,

Т^з_рп - затраты времени разработчика программного обеспечения на разработку ТЗ, чел.-дни.

Значения Т^з_рз и Т^з_рп рассчитываются по формуле

Т^з_рз = t_зК^з_рз, Т^з_рп = t_зК^з_рп,

где t_з- норма времени на разработку ТЗ для программного продукта в зависимости от функционального назначения и степени новизны разрабатываемового программного продукта, чел.-дни (t_з = 29),

К^з_рз- коэффициент, учитывающий удельный вес трудоемкости работ, выполняемых разработчиком постановки задач на стадии ТЗ (К^з_рз= 0,65),

К^з_рп- коэффициент, учитывающий удельный вес трудоемкости работ, выполняемых разработчиком программного обеспечения на стадии ТЗ (К^з_рп= 0,35).

Т^з_рз = 29*0,65 = 18,85 чел.-дней.

Т^з_рз = 29*0,35 = 10,15 чел.-дней.

t_тз= Т^з_рз + Т^з_рп =18,85 + 10,15 = 29 чел.-дней.

Трудоемкость разработки эскизного проекта рассчитывается по формуле

t_эп= Т^э_рз + Т^э_рп,

где Т^э_рз - затраты времени разработчика постановки задач на разработку ЭП, чел.-дни,

Т^э_рп - затраты времени разработчика программного обеспечения на разработку ЭП, чел.-дни.

Значения Т^з_рз и Т^з_рп рассчитываются по формуле

Т^э_рз = t_эК^э_рз, Т^э_рп = t_эК^э_рп,

где t_э- норма времени на разработку ЭП для программного продукта в зависимости от функционального назначения и степени новизны разрабатываемового программного продукта, чел.-дни (t_э = 41),

К^э_рз- коэффициент, учитывающий удельный вес трудоемкости работ, выполняемых разработчиком постановки задач на стадии ЭП (К^э_рз= 0,7),

К^э_рп- коэффициент, учитывающий удельный вес трудоемкости работ, выполняемых разработчиком программного обеспечения на стадии ЭП (К^э_рп= 0,3).

Т^з_рз = 41*0,7 = 28,7 чел.-дней.

Т^з_рз = 41*0,3 = 12,3 чел.-дней.

tэ_п= Т^з_рз + Т^з_рп = 28,7 + 12,3 = 41 чел.-дней.

Трудоемкость разработки технического проекта зависит от функционального назначения программного продукта, количества разновидностей входной и выходной информации и определяется как сумма времени, затраченного разработчиком постановки задач и разработчиком программного обеспечения:

t_тп= (t^т_рз + t^т_рп)К_вК_р,

где t^п_рз, t^п_рп- норма времени на разработку ТП разработчиком постановки задач и разработчиком программного обеспечения соответственно, чел.-дни (t^т_рз = 9, t^т_рп = 8),

К_в- коэффициент учета вида используемой информации,

К_р- коэффициент учета режима обработки информации (К_р = 1,1).

Значение коэффициента К_вопределяется по формуле

К_в= (К_пn_п + К_нсn_нс + К_бn_б)/(n_п + n_нс + n_б),

где К_п, К_нс, К_б - значения коэффициентов учета вида используемой информации для переменной, нормативно-справочной информации и баз данных соответственно (К_п = 0,5, К_нс = 0,43, К_б= 1,25),

n_п, n_нс, n_б- количество наборов данных переменной, нормативно-справочной информации и баз данных соответственно (n_п = 1, n_нс = 1, n_б= 0).

К_в= (0,5 + 0,43)/2 = 0,465

t_тп= (8 + 9)*0,465*1,1 = 8,6955 чел.-дней.

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

t_рп= (t^р_рз + t^р_рп)К_кК_рК_яК_зК_иа,

где t^р_рз, t^р_рп- норма времени на разработку РП разработчиком постановки задач и разработчиком программного обеспечения соответственно, чел.-дни (t^р_рз = 5, t^р_рп = 27),

К_к- коэффициент учета сложности контроля информации (К_к = 1,07),

К_я- коэффициент учета уровня используемового языка программирования (К_я = 1,0),

К_з- коэффициент учета степени использования готовых программных модулей (К_з = 0,8),

К_иа- коэффициент учета вида используемой информации, и сложности алгоритма программного продукта.

Значение коэффициента К_иаопределяется по формуле

К_иа= (К’_пn_п + К’_нсn_нс + К’_бn_б)/(n_п + n_нс + n_б),

где К’_п, К’_нс, К’_б - значения коэффициентов учета сложности алгоритма программного продукта и вида используемой информации для переменной, нормативно-справочной информации и баз данных соответственно (К’_п = 0,48, К’_нс = 0,29, К’_б= 0,24),

n_п, n_нс, n_б- количество наборов данных переменной, нормативно-справочной информации и баз данных соответственно (n_п = 1, n_нс = 1, n_б= 0).

К_иа= (0,48 + 0,29)/2 = 0,385

t_рп= (5 + 27)*1,07*1,1*0,8*0,385 = 11,6 чел.-дней.

Трудоемкость внедрения может быть рассчитана по формуле:

t_в= (t^в_рз + t^в_рп) К_кК_рК_з,

где t^в_рз, t^в_рп- норма времени на внедрение программного продукта разработчиком постановки задач и разработчиком программного обеспечения соответственно, чел.-дни (t^в_рз = 8, t^в_рп = 24).

t_тп= (8 + 24)*0,8*1,07 = 27,392 чел.-дней.

t_пп = 29 + 41 + 8,6955 + 11,6 + 27,392 = 117,6875 чел.-дней.

Продолжительность выполнения всех работ по этапам разработки программного продукта рассчитывается по формуле

T_i = (t_i + Q)/n_i,

где t_i - трудоемкость i-й работы, чел.-дни,

Q - трудоемкость дополнительных работ, выполняемых исполнителем, чел.-дни,

n_i - количество исполнителей, выполняющих i-ю работу, чел.

T_тз = t_тз/2 = 29/2 = 14,5 (15) дней

T_эп = t_эп/2 = 41/2 = 20,5 (21) дней

T_тп = t_тп = 8,6955 (9) дней

T_рп = t_рп/2 = 11,6/2 = 5,8 (6) дней

T_в = t_в = 27,392 (28) дней

T_пп = SТ_i = 15 + 21 + 9 + 6 + 28 = 79 дней.
3.3. РАСЧЕТ СМЕТЫ ЗАТРАТ НА РАЗРАБОТКУ ПРОГРАММНОГО ПРОДУКТА
Смета затрат на выполнение работ составляется по калькуляционным статьям. В общем случае статьи, учитывающие расходы, следующие:

1. Материалы (суммарные затраты на материалы, приобретаемые для разработки программного продукта).

2. Специальные оборудование (суммарные расходы на аренду приборов, требуемых для разработки программного продукта).

3. Основная заработная плата производственного персонала.

4. Дополнительная заработная плата.

5. Отчисление на социальное страхование.

6. Накладные расходы.

7. Производственные командировки.

8. Контрагентские расходы.

Однако затраты, связанные с разработкой программного обеспечения, носят специфический характер. Расходы по статьям 7, 8 обычно крайне незначительны. Статьи 1, 2 связаны с расходами на использование ЭВМ. Эти расходы определяются, исходя из затрат машинного времени и стоимости часа работ ЭВМ, а также стоимости необходимых материалов и покупных изделий, необходимых при работе на ЭВМ.

В результате можно определить следующие статьи расходов на разработку программного продукта:

1. Стоимость машинного времени, затраченного на разработку.

2. Стоимость материалов и покупных изделий.

3. Основная заработная плата исполнителей.

4. Дополнительная заработная плата.

5. Отчисление на социальное страхование.

6. Накладные расходы.
Расчет стоимости затраченного машинного времени.

С_эвм = t^эвмК^эвм_иЦ^эвмК^эвм_бдК^эвм_э,

где t^эвм- время использования ЭВМ для создания данного программного продукта, час (t^эвм= 10).

К^эвм_и- поправочный коэффициент учета времени использования ЭВМ (К^эвм_и= 0,6).

Ц^эвм- цена одного часа работы ЭВМ, руб (на январь 1997 года 8000 рублей).

К^эвм_бд- коэффициент учета степени использования СУБД (К^эвм_бд= 1,0).

К^эвм_э- коэффициент учета быстодействия ЭВМ (К^эвм_и= 1,2).

С_эвм = 10*0,6*1,2*8000 = 57600 рублей.
Расчет затрат на материалы.

В процессе разработки программного изделия используются следующие материалы:

бумага формата А4 в количестве 500 листов общей стоимости 60000 рублей.

дискеты 3,5 дюйма в количестве 2 штук по цене 5000 рублей за штуку.

С_м = 60000 + 5000*2 = 70000 рублей.
Таблица 4.2. Расходные материалы.

Материал	Цена, руб	Количество, шт.	Стоимость, руб.
дискета 1,44 Мб	5000	2	10000
бумага формата А4		500	60000
Итого:			70000

Расчет основной заработной платы исполнителя.

С_зо = Sз_it_i/d,

где з_i - средняя заработная плата i-го исполнителя (300000 рублей в месяц),

t_i - трудоемкость работ, выполняемых i-м исполнителем (чел.-дни),

d - среднее количество рабочих дней в месяце (d = 22).

С_зо = 79*300000/22 = 1077300 рублей.
Расчет дополнительной заработной платы.

В статье “дополнительная заработная плата” учитываются выплаты, предусмотренные законодательством о труде и коллективными договорами за непроработанное на производстве время.

Дополнительная заработная плата определяется по установленному нормативу от основной заработной платы по формуле

С_зд = С_зоa,

где a - коэффициент дополнительной заработной платы, a = 0,2.

С_зд = 1077300*0,2 = 215500 рублей.
Расчет отчислений на социальное страхование.

В статью “отчисления на социальное страхование” включено отчисление по единому установленному нормативу от суммы основной и дополнительной заработной платы.

Размер отчислений вычисляется по формуле

С_сс = (С_зд + С_зо)a_сс,

где a_сс - коэффициент, устанавливающий отчисление на социальное страхование и в фонд стабилизации, a_сс = 0,4.

С_сс = (1077300 + 215500)*0,4 = 517120 рублей.
Расчет накладных расходов.

В статье “накладные расходы” учитываются командировочные расходы, оплата подъемных при перемещениях, арендная плата, оплата услуг сторонних организаций.

С_н = С_зоa_н,

где a_н - коэффициент накладных расходов, a_н= 1,8.

С_н = 1,8*1077300 = 1939140 рублей
Расчет суммарных расходов.

С = С_эвм + С_м + С_зо + С_зд + С_сс + С_н =

= 57,6 +70 + 1077,3 + 215,5 + 517,12 + 1939,14 = 3876,66 тыс.рублей.
Смета затрат на разработку программного продукта приведена в таблице 4.3.

Таблица 4.3.

№ п/п	Наименование статей расходов	Затраты (тыс.руб.)	Удельный вес, %
1	Стоимость машинного времени	57,6	1,4
2	Материалы	70	1,9
3	Основная заработная плата	1077,3	27,7
4	Дополнительная заработная плата	215,5	5,5
5	Отчисления на социальное страхование	517,12	13,4
6	Накладные расходы	1939,14	50,1
	Итого:	3876,66

4. ПРОМЫШЛЕННАЯ ЭКОЛОГИЯ И БЕЗОПАСНОСТЬ
4.1. ВВЕДЕНИЕ
В результате развития производственных сил общества возникла проблема взаимодействия человека и машины. Охрана труда и эргономика позволяют с научной точки зрения подойти к этой проблеме, способствуют изучению влияния окружающей среды на человека, который непосредственно контактирует с ЭВМ, определению вредных и опасных производственных факторов, разрабатывают организационно-технические мероприятия, направленные на профилактику профессиональных заболеваний, создавая здоровые и безопасные условия труда для работающего.

Предметом исследования эргономики в этой области стало согласование психо-физических возможностей человека со свойствами современных технических систем. Только в этом случае можно рассчитывать на высокое качество и эффективность его труда. Особую актуальность эта проблема приобретает в связи с возросшим культурным уровнем современного персонала, предъявляюшего повышенные требования к содержанию и условиям труда на рабочем месте (РМ).

Под рабочим местом в эргатических системах (ЭС) согласно ГОСТ 26387-84 понимается «часть пространства в системе человек-машина (СЧМ), оснащенная средствами отображения информации, органами управления, вспомогательным оборудованием и предназначенная для осуществления деятельности оператора СЧМ». Соответственно среда на РМ определяется этим же ГОСТом как «совокупность физических, химических и психологических факторов, воздействующих на оператора СЧМ, на его РМ в ходе его деятельности».
4.2. АНАЛИЗ ВРЕДНЫХ ФАКТОРОВ
Нормальная и безопасная работа инженера-программиста за экраном дисплея во многом зависит от того, в какой мере условия его работы соответствуют оптимальным. При этом под условиями работы подразумевают комплекс физических, химических, биологических и психофизических факторов, установленных стандартами по безопасности труда (ССТБ).

К физическим факторам относятся:

- вибрация и шум из-за движущихся машин, механизмов и их элементов, запыленность и загазованность воздуха, температура поверхностей оборудования, материалов и воздуха;

- плотность воздуха, ее резкое изменение, подвижность и ионизация воздуха;

- ионизирующие и электромагнитные излучения, статические заряды и повышение напряжения в цепи, электрические и магнитные поля;

- отсутствие или недостаток естественного света, повышенная или пониженная освещенность, яркость и контрастность, блесткость поверхности, пульсация светового потока;

- ультрафиолетовое или инфракрасное излучение.

К химическим факторам относятся:

- общетоксические, раздражающие, сенсибилизирующие, канцерогенные, мутагенные;

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

К биологическим факторам относятся:

- микроорганизмы (бактерии, вирусы, грибы и т.д.);

- макроорганизмы (растения и животные).

К психофизическим факторам относятся перегрузки:

- физические (статические, динамические, гиподинамия);

- нервно-психические (умственное перенапряжение, монотонность труда, эмоциональные перегрузки).

При проектировании рабочего места инженера-программиста необходимо учитывать и нормировать все указанные группы факторов, поскольку при определенных условиях они могут вызвать нежелательные функциональные сдвиги в организме оператора, снизить качество и эффективность его работы, оказать отрицательное влияние на его здоровье.

Наиболее значительным фактором является микроклимат, особенно температура и влажность воздуха. Исследования показывают, что высокая температура в сочетании с высокой влажностью воздуха оказывают большое влияние на работоспособность человека. Резко увеличивается время сенсорных и моторных реакций, нарушается координация движений, увеличивается количество ошибок. Высокая температура отрицательно сказывается и на ряде психологических функций человека. Уменьшается объем оперативной памяти, резко суживается способность к ассоциациям. При +11⁰С начинается окоченение конечностей, такая температура минимально допустима. Наиболее благоприятный диапазон температур в летнее время от +18⁰С до +24⁰С, в зимнее время от +17⁰ до +22⁰С.

Движение воздуха позволяет увеличить рабочий диапазон температур. Так при скорости движения воздуха 0.1, 0.5, 0.9 м/с верхняя допустимая граница рабочего диапазона сдвигается соответственно до +22⁰, +24⁰, +26⁰С при интенсивном расходе энергии человеком порядка 1000 Дж/ч.

Атмосферное давление в пределах 80-106 кПа легко переносимо человеком. При давлениях, выходящих за эти пределы, человеку требуется предварительная акклиматизация.

Результаты работы инженера-программмиста в большой степени зависят и от освещенности рабочего места. Чтобы правильно спланировать рациональную систему освещения, необходимо учитывать яркость источников света, их расположение в помещении, яркостной контраст между устройствами ЭВМ и фоном, блесткость поверхностей, качество и цвет светильников и поверхностей. Для малой и средней контрастности поверхностей ЭВМ при темном фоне наименьший уровень освещенности должен быть 150 лк. Для большой контрастности при светлом или темном фоне наименьший уровень освещенности 100 лк.

В помещениях, где эксплуатируют ЭВМ, необходимо предусматривать систему искусственного освещения из люминисцентных ламп дневного света или ламп накаливания. Существуют прямая, отраженная и диффузная системы искусственного освещения. При прямом освещении свет попадает на объект непосредственно от источников света. При этом 90-100% мощности светильника направлено на рабочую поверхность, что вызывает яркостные контрасты, резкие тени и блесткость (свойство ярко освещенной поверхности вызывать ослепление или дезадаптацию наблюдателя). При освещении отраженным светом 90-100% света направляется на потолок и верхнюю часть стен, от которых свет более или менее равномерно отражается по всему помещению. При этом достигается равная освещенность без теней и блесткости. Диффузное освещение обеспечивает рассеянный свет, одинаково распределенный по всем направлениям. Такая система освещения требует меньшей мощности, чем две предыдущие, но вызывает частичное образоование теней и блесткости.

Кроме освещенности, большое влияние на деятельность человека оказывает цвет окраски помещения и спектральные характеристики используемого цвета. Рекомендуется, чтобы потолок отражал 80-90%, стены - 50-60%, панели - 15-20%, а пол - 15-30% падающего на них света. Кроме того, цвет обладает некоторым психологическим и физиологическим действием. Так, например, применение тонов теплой гаммы (красный, оранжевый, желтый) создает впечатление бодрости, возбуждения и замедленного течения времени. Эти же цвета вызывают у человека ощущение тепла.

Большое влияние на деятельность инженера-программиста оказывает и уровень акустического шума. Шум резко снижает производительность труда и увеличивает травматизм. Физиологически шум воздействует на органы зрения и слуха, повышает кровяное давление, при этом притупляется внимание.

Шум оказывает также и эмоциональное воздействие: он является причиной возникновения таких отрицательных эмоций, как досада, раздражение. Особенно неприятны высокочастотные и прерывистые шумы.

Основным из механических факторов производственной среды являются вибрации. Они не только вредно воздействуют на организм, но и мешают человеку выполнять как мыслительные так и двигательные операции. Под действием вибраций ухудшается зрительное восприятие, в осообенности на частотах между 25 и 40 Гц и между 60 и 90 Гц. Наиболее опасна вибрация с частотой 6-8 Гц, так как в этом диапазоне лежит собственная резонансная частота тела, головы и брюшной полости человека.

К числу неблагоприятных факторов относятся злектромагнитные поля (ЭМП) высоких частот. Их воздействие на человека может вызвать функциональные сдвиги в организме: быструю утомляемость, головные боли, нарушение сна, раздражительность, утомление зрения и т.п.

Предельно допустимые уровни ЭМП следующие:

- в СВЧ-диапазоне - мкВт/см;

- в диапазоне до 300 МГц по электрической составляющей - 5 В/м, по магнитной составляющей - 5 А/м. С учетом этого стандарта было исследовано свыше 150 мониторов различных типов.

На жизнедеятельность человека большое влияние оказывает газовый состав воздуха. Здесь обычно исследуется две группы факторов: изменение обычного состава воздуха (кислорода и углекислого газа) и посторонние добавки к нему в результате работы техники.

Благоприятными условиями газового состава воздуха считается содержание кислорода 19-20%, углекислого газа около 1%; допустимые значения, при которых не происходит выраженного снижения работоспособности составляют: кислорода - 18-29%, углекислого газа - 1-2%. Снижение содержания кислорода ниже 16% и повышение содержания углекислого газа выше 3% являются недопустимыми и могут привести к нежелательным последствиям. Важнейшим способом борьбы с неблагоприятным воздействием на человека химических факторов является соблюдение их предельно допустимых концентраций в производственных помещениях. Предельно допустимыми считаются такие максимальные концентрации вредных веществ, которые при ежедневной работе не могут вызывать у работающих заболевания или отклонения в состоянии здоровья. Такими концентрациями считаются, например, для аммиака - 20 мг/м, анилина - 3 мг/м, ацетона - 200 мг/м, бензола - 5 мг/м, бензина - 100 мг/м, серной кислоты - 1 мг/м.

При выполнении данной дипломной работы используются следующие элементы вычислительной техники:

персональный компьютер IBM PC 486DX;

струйный принтер Canon Bubble Jet.

Персональный компьютер питается напряжением 220В/50Гц, которое превышает безопасный предел 42 В. Следовательно возникает опасность поражения электрическим током.

Воздействие на человека электрического тока приводит к общим травмам (электроудары) и местным (ожоги, металлизация кожи, электрические знаки, электроофтальмия, механические повреждения).

Возникновение рентгеновского излучения обусловлено наличием на аноде электронно-лучевой трубки дисплея напряжения до 30 кВ (а при напряжении 3-500 кВ присутствует рентгеновское излучение различной жесткости). Пользователь попадает в зону мягкого рентгеновского излучения.

При воздействии рентгеновского излучения на организм человека происходит:

образование чужеродных соединений молекул белка, обладающих даже токсическими свойствами;

изменение внутренней структуры веществ в организме, приводящее к развитию малокровия, образованию злокачественных опухолей, катаракты глаз.

При работе за экраном электронно-лучевой трубки дисплея пользователь попадает под воздействие ультрафиолетового излучения с длинами волн < 320 нм. Также при образовании строчной и кадровой разверток дисплея возникает излучение электромагнитных полей частотой до 100 кГц. Это может являться причиной возникновения следующих заболеваний:

обострение некоторых заболеваний кожи (угревая сыпь, себорроидная экзема, розовый лишай, рак кожи и др.);

нарушение в протекании беременности;

увеличение в 2 раза вероятности выкидышей у беременных женщин;

нарушение репродуктивной функции и возникновение рака;

нарушение режима терморегуляции организма;

изменения в нервной системе (потеря порога чувствительности);

понижение/повышение артериального давления.

При работе на персональном компьютере человек попадает под воздействие статического электричества. Под действием статических электрических полей дисплея пыль помещения электризуется и переносится на лицо пользователя, что приводит к заболеваниям (раздражению) кожи (дерматит, угри).

При работе за персональным компьютером для вывода информации на бумажный носитель применяется принтер. Принтер Canon Bubble Jet имеет уровень звука на расстоянии 1 метр от корпуса 49 дБ (используется 1 час в течении смены), что соответствует норме. Следовательно, вредного воздействия по звуку на пользователя не оказывается.

Таким образом пользователь, работающий с персональным компьютером подвергается воздействию следующих опасных и вредных факторов:

поражение электрическим током;

воздействие рентгеновского излучения;

ультрафиолетовое излучение и излучение электромагнитных полей радиочастот;

воздействие статического электричества.
4.3. ТРЕБОВАНИЯ К ВИДЕОТЕРМИНАЛЬНЫМ УСТРОЙСТВАМ
Основными поражающими факторами, при работе с компьютером, являются вредные излучения видеотерминального устройства.

Видеотерминальное устройство должно соответствовать следующим требованиям:

яркость свечения экрана не менее 100 кд/м²;

минимальный размер светящейся точки не более 0,4 мм для монохромного дисплея и не более 0,6 мм для цветного;

контрастность изображения знака не менее 0,8;

частота регенерации изображения при работе с позитивным контрастом в режиме обработки текста не менее 72 Гц;

количество точек на экране не менее 640;

экран должен иметь антибликовое покрытие;

размер экрана должен быть не менее 31 см по диагонали, а высота символов не менее 3,8 мм, при этом расстояние от экрана до глаз оператора должно быть 40–80 см.

При работе с текстовой информацией наиболее предпочтительным является предъявление чёрных знаков на светлом (белом) фоне.

Максимальные значения напряженности магнитного поля, измеренные на расстоянии 50 см от экранов наиболее распространённых мониторов.
Полоса частот	Магнитное поле, А/м	Нормы BGA
5 - 1000 Гц	0,2	160 - 0,8
10 - 150 кГц	0,17	0,8 - 0,6
150 - 300 кГц	-	0,6 - 0,42
0,3 - 30 Мгц	0,00000066	0,42 - 0,73
30 - 300 Мгц	0,00000066	0,73

Максимальная напряженность электрического поля, допускаемая нормами BGA, равна 2,5 кВ/м. Это значение установлено из расчёта того, чтобы при прикосновении к заряженной проводящей поверхности электрический разряд не стал причиной шока.

Максимальные значения напряженности электрического поля, измеренные на расстоянии 50 см от экранов наиболее распространённых мониторов.
Полоса частот	Электрическое поле, В/м	Нормы BGA
5 - 1000 Гц	4,8	2500 - 177
10 - 150 кГц	4,8	87
150 - 300 кГц	0,48	87
0,3 - 30 Мгц	0,0024	87 - 27,5
30 - 300 Мгц	0,0024	27,5

Измерения BGA показывают, что напряженность электростатического поля около монитора может превысить 7 кВ/м. Согласно полученным SSI и SEMKO (Швеция) данным, эти значения для некоторых устройств достигают 50 кВ/м.

В России нормирование электромагнитных полей осуществляется в соответствии с ГОСТ 12.1.006-84 и санитарными нормами СНиП2963-84.

В зоне индукции нормируется напряженность электрического и магнитного поля в зависимости от частоты. В зоне излучения нормируется плотность потока энергии в зависимости от времени пребывания.

Нормир.	Частота f, МГц
велич.	0.06-1.5	1.5-3.0	3.0-30	30-50	50-300	300-3*10⁵
Е, В/м	50	50	20	10	5	нет
Н, В/м	5.0	–	–	0.3	–	нет
I, Вб/м²	–	–	–	–	–	I₀ = e/T

Электромагнитные поля нормируются следующим образом:

электрические: E = 6/ÖT; 1 £ T £ 9, где Т- время воздействия;

магнитные: Hn £ 8 кА/м в течение рабочего дня; e = 2 (Вт r/м²) – энергетическая нагрузка на организм.
4.4. РАСЧЕТ ВРЕДНЫХ ИЗЛУЧЕНИЙ ВИДЕОДИСПЛЕЯ
Время работы на персональном компьютере по санитарным нормам не должно превышать 4 часа.Большинство используемых в России мониторов не соответствуют шведскому стандарту защита пользователя от излучений и имеют на расстоянии 5 см от экрана дисплея имеют мощность дозы рентгеновского излучения 100 мкР/час. Рассчитаем, какую дозу рентгеновского излучения получит пользователь на различном расстоянии от экрана дисплея.

P_r = P₀e^-^m^r, где

P_r- мощность дозы рентгеновского излучения на расстоянии r, мкР/час;

P₀ - уровень мощности дозы рентгеновского излучения на расстоянии 5 см от экрана дисплея, мкР/ч.

m - линейный коэффициент ослабления рентгеновского излучения воздухом, см^-1;

r - расстояние от экрана дисплея, см;

Возьмем m = 3.14*10^-2 см^-1.

r, см	5	10	20	30	40	50	60	70	80	90	100
P, мкР/ч	100	73	53	39	28	21	15	11	8	6	4

Среднестатистический пользователь располагается на расстоянии 50 см от экрана дисплея. Рассчитаем дозу облучения, которую получит пользователь за смену, за неделю, за год.

За смену	4 часа	4*21	84 мкР/ч
За неделю	5 дней	5*84	420 мкР/ч
За год	44 рабочие недели	44*420	18480 мкР/ч

4.5. РАЦИОНАЛЬНАЯ ОРГАНИЗАЦИЯ РАБОЧЕГО МЕСТА
Для повышения производительности труда при работе за компьютером необходимо создать на рабочем месте наиболее благоприятные условия с точки зрения эргономики и эстетики.

Разработка мероприятий по рациональной организации рабочего места инженера-программиста и инженера-разработчика может идти в следующих направлениях:

устранение неблагоприятных факторов:

снижение шума в помещении;

правильный выбор источников освещения;

устранение запылённости и загазованности.

оптимизация условий труда на рабочем месте:

эргономические требования;

психологические требования.

создание комфортных условий отдыха в течение рабочего дня.

Производственные помещения вычислительного центра должны проектироваться в соответствии с требованиями СНиП 2.03.04-87 – “Административные и бытовые здания и помещения производственных предприятий”.

Площадь помещения следует принимать из расчёта 6 м² на одного работника. При оснащении рабочих мест терминалами ЭВМ, печатающими устройствами и пр. площади помещения допускается увеличивать в соответствии с техническими условиями на эксплуатацию оборудования. Кубатура должна быть не менее 19,5 м³ с учётом максимального числа одновременно работающих.

Минимальная ширина проходов с передней стороны пультов и панелей управления ЭВМ при однорядном расположении должна быть не менее 1 м, при 2-х рядном расположении не менее 1,2 м. Видеотерминалы должны располагаться при однорядном размещении на расстоянии не менее 1 м от стен. Рабочие места с дисплеями должны располагаться между собой на расстоянии не менее 1,5 м.

На постоянных рабочих местах и в кабинах операторов должны быть обеспечены микроклиматические параметры, уровни освещённости, шума и состояния воздушной среды, определённые действующими санитарными правилами и нормами.
4.6. РЕКОМЕНДАЦИИ ПО СНИЖЕНИЮ УТОМЛЯЕМОСТИ
Необходимо расположить экран дисплея немного выше уровня глаз. Это создаст разгрузку тех групп окологлазных мышц, которые наиболее напряжены при обычном направлении взгляда - вниз или вперёд.

Помещение, где находятся компьютеры и видеомониторы, должно быть достаточно просторным с постоянным обновлением микроатмосферы. Минимальная площадь на один видеомонитор - 9-10 м². Крайне нежелателен визуальный контакт работника с другими мониторами или телевизионными экранами. Необходимо исключить наличие всевозможных бликов на экране монитора, часто возникающих на стеклянных экранах. Следует также избегать большой контрастности между яркостью экрана и окружающего пространства - оптимальным считается выравнивание яркости экрана и компьютера. Запрещается работа с компьютером в тёмном или полутёмном помещении.

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

Для большего эргономического комфорта целесообразно расположить в кресле опору - в районе поясничного изгиба позвоночника (в виде продолговатой подушечки или валика).

Если работник имеет те или иные рефракционные отклонения (близорукость, дальнозоркость и др.), то последние должны быть полностью коррегированы очками. При более серьёзных отклонениях вопрос о возможности работы с видеотерминалами должен решаться с участием врача-офтальмолога.

Через каждые 40-45 минут необходимо проводить физкультурную микропаузу: вращение глаз по часовой стрелке и обратно, лёгкие гимнастические упражнения для всего тела, например поднимание и опускание рук.

Каждый час необходимо делать перерыв и выполнять несколько упражнений на расслабление, которые могут уменьшить напряжение, накапливающиеся в мышцах при длительной работе за компьютером.
4.7. ЗАЩИТА ОТ НАПРЯЖЕНИЯ ПРИКОСНОВЕНИЯ. ЗАНУЛЕНИЕ
Занулением называется преднамеренное соединение нетоковедущих частей с нулевым защитным проводником (НЗП). Оно применяется в трехфазных сетях с глухозаземленной нейтралью в установках до 1000 вольт и является основным средством обеспечения электробезопасности.

При попадании напряжения сети на корпус ПЭВМ возникает режим короткого замыкания. Для защиты электрической сети от короткого замыкания и перегрузок применяются автоматические выключатели или предохранители. При проектировании защитного устройства необходимо рассчитать его номинальный ток срабатывания - I_ном:

I_alarm ³ KI_ном, где

I_ном = I_alarm/K

I_ном- номинальный ток срабатывания защитного устройства, A;

K - коэффициент, учитывающий тип защитного устройства:

K = 3 - для автомата с электромагнитным расцепителем,

K = 1.4 - для теплового автомата,

I_alarm- ток короткого замыкания, A.

Рассчитаем величину тока короткого замыкания:

I_alarm= U_f/(R_n + R_m/3)

R_n = R_f + R₁ + jx₁

U_f = 220 В

R_m = 0,312W

R_f = 0,412W

jx₁= 0,6W

R₁ = r/S

r - удельное сопротивление НЗП, [Wmm²/m];

l - длина НЗП, m;

r_cu = 0,0175 W mm² /m,

l = 50 m,

S = 1,5 mm²

R₁ = 0,0175(50/15) = 0,58W

R_n = (0,412 + 0,58 + 0,6) = 1,59W

I_alarm = 130 A

I_ном = 43 A

Для того, чтобы в случае короткого замыкания или других причин ПЭВМ отключалась от электрической сети необходимо в цепь питания поставить автомат с электромагнитным расцепителем с I_ном = 43 A.
4.8. ПОЖАРНАЯ БЕЗОПАСНОСТЬ
В помещениях ВЦ существуют все три основные фактора, необходимые для возникновения пожара.

Горючими материалами на ВЦ являются: строительные материалы для акустической и эстетической отделки помещений, перегородки, двери, полы, изоляция силовых и сигнальных кабелей, шкафы, жидкости для очистки элементов и узлов ЭВМ и т.д.

Для отвода тепла от ЭВМ в производственных помещениях ВЦ постоянно действует система кондиционирования. Поэтому кислород, как окислитель процессов горения, имеется в любой точке помещений ВЦ.

Источниками зажигания на ВЦ могут оказаться электронные схемы ЭВМ, приборы, приборы, применяемые для технического обслуживания, устройства электропитания, кондиционеры воздуха.

По пожарной опасности ВЦ относятся к категории “В” (в производстве обращаются твердые сгораемые вещества и материалы). Исходя из этого ВЦ проектируется с II степенью огнеустойчивости.

Минимальные пределы огнеустойчивости в часах:

Cтепень огнестойкости зданий и сооружений	II
Основные строительные конструкции:
Несущие стены, стены лестничных клеток, колонны	2
Лестничные площадки	1
Наружние стены из навесных панелей	0,25
Внутренние несущие стены, перегородки	0,25
Несущие конструкции междуэтажных перекрытий	0,75
Плиты, настилы и др.	0,25

Для обнаружения начальной стадии загорания используют систему автоматической пожарной сигнализации (АПС). АПС состоит из пожарных извещателей, линий связи и приемных пультов (станций).

В помещениях ВЦ применят дымовые пожарные извещатели типа РИД-1.

Принцип действия РИД-1 основан на изменении величины электрического тока, протекающего через ионизационную камеру, при попадании в нее дыма.

Технические показатели для РИД-1:

чувствительный элемент	ионизационная камера
параметр срабатывания	тлеющий фитиль
инерционность, сек	10
диапазон температур, С	-30 ...+50
относительная влажность, %	80

Норма расстановки пожарных извещателей в помещениях с гладким полом:

Тип	Защищаемая площадь, i²	Расстояние между извещателями, м
		максимальное	в узких коридорах
РИД-1	100	12	15

Линии связи систем АПС с приемными станциями строятся по лучевому принципу. Приемные станции АПС устанавливаются в помещении дежурного по ВЦ, где организуется круглосуточное дежурство.

Приемные станции обеспечивают следующие функции:

прием сигналов от пожарных извещателей с индикацией номера луча;

непрерывный контроль состояния лучей по всей длине с автоматическим выявлением характера повреждения;

световая и звуковая сигнализация тревоги;

автоматическое переключение на резервный источник питания при сбоях сети с включением соответствующей сигнализации.

На ВЦ используется приемная станция РОУП-1.

Технические характеристики устройства РОУП-1:

извещателей РИД-1, шт	до 300
шлейфов блокировки, компл.	до 30
напряжение питания, В	220±10
потребляемая мощность, Вт	не более 180
диапазон температур, С	+5 ... +50
относительная влажность, %	до 80
срок службы, лет	8
дополнительные функции	может управлять устройствами пожаротушения

На ВЦ применяются установки газового тушения пожара, действие которых основано на быстром заполнении помещения газом с низким содержанием кислорода. Используется автоматическая установка газового пожаротушения (АУГП) с электрическим пуском.

Технические характеристики АУГП с электрическим пуском:

число пусковых баллонов, шт	2
число рабочих баллонов, шт	4
заряд пускового баллона	сжатый воздух
заряд рабочего баллона	фреон 114 Вч
вместимость пускового баллона, л	27
вместимость рабочего баллона, л	40
давление в пусковом баллоне, МПа	125
давление в рабочем баллоне, МПа	12,5
продолжительность пуска, с	65
масса батареи без заряда, кг	480

При использовании АУГП для предотвращения отравления персонала предусмотрена предупредительная звуковая и световая сигнализация, срабатывающая при ручном, дистанционном и автоматическом включении за 30 секунд до начала выпуска газа.

Расчет необходимого количества баллонов с сжатым воздухом и огнегасительным средством:

Количество огнегасительного вещества (фреона)

G_т = G_вW_пK_у, где G_т - количество огнегасительного вещества,

W_п - расчетный объем защищаемого помещения, м³,

G_в- огнегасительная концентрация газового состава, кг/м³,

K_у - коэффициент, учитывающий особенности процессов газообмена в защищаемом помещении.

Для ВЦ G_в= 0,25 кг/м³, K_у= 1,2.

W_п= SH, где S - площадь помещения, м².

H - высота помещения, м.

S = 100 м². H = 3 м. W_п= 300 м³.

G_т = 0,25*300*1,2 = 90 кг.

Необходимое количество баллонов

N_б = G_т/V_бra, где V_б - объем баллона, м³,

r - плотность, кг/л,

a - коэффициент наполнения баллона.

V_б = 40 л, r = 2,17 кг/л, a = 0,9.

N_б = 90/(40*2,17*0,9) = 2.

Объем воздушных баллонов

W_б= (Р^с_мин+1)(W_с+W_т)/(Р_макс-Р^б_мин), где Р^с_мин и Р^б_мин- конечное давление в воздушных баллонах и баллонах с огнегасительным средством, МПа,

Р_макс- начальное давление воздуха в баллоне, МПа,

W_сиW_т- объем баллонов с огнегасительным составом и трубопроводов, л.

Р^с_мин = Р^б_мин= 5 Мпа, W_с= 2*40 = 80 л, W_т= 20л, Р_макс = 125 МПа.

W_б= (5+1)(80+20)/(125-5) = 4,8 л.

5. СПИСОК ЛИТЕРАТУРЫ.

1. «Основы теории полета космических аппаратов» / Под ред. Г.С.Нариманова, М.К.Тихонравова. М., Машиностроение, 1972.

2. А.П.Разыграев «Основы управления полетом космических аппаратов». М., Машиностроение, 1990.

3. Г.Г.Бебенин, Б.С.Скребишевский, Г.А.Соколов «Системы управления полетом космических аппаратов». М., Машиностроение, 1978.

4. К.Б.Алексеев, Г.Г.Бебенин «Управление космическими летательными аппаратами». М., Машиностроение, 1974.

5. В.В.Солодовников, В.Н.Плотников, К.В.Яковлев «Теория автоматического управления технических систем». М., изд.МГТУ им.Баумана, 1993.

6. Б.Страуструп «Язык программирования С++». М., «Радио и связь», 1991.

7. А.В.Бошкин, П.Н.Дубнер «Работа с С++». М., «Юкис», 1991.

8. В.В.Арсеньев, Б.Ю.Сажин «Методические указания к выполнению организационно-экономической части дипломных проектов по созданию программной продукции», М., изд. МГТУ им.Баумана, 1994.

9. ГОСТ 2.103-68 НИР. М.: Изд-во стандартов, 1968.

10. В.К.Зелинский «НОТ в проектно-конструкторской организации». М.: «Экономика», 1969.

11. «Управление трудовым коллективом» / Г.П.Зайцев, Э.В.Минько, Н.В.Артамонова и др. Свердловск, Изд-во УГУ, 1989.

12. «Типовые нормы времени на программирование задач для ЭВМ», утвержденные постановлением Государственного комитета СССР по труду и социальным вопросам и Секретариата ВЦСПС от 27 июля 1987 г. №454/22-70

13. Ю.Г.Сибиров «Охрана труда в ВЦ». М., «Машиностроение», 1985.

14. Сибиров Ю.Г., «Основы инженерной психологии» / под ред. Б.Ф.Ломова. М., «Машиностроение», 1986.

15. СНиП 2.09.04-87 «Административные и бытовые здания и помещения производственных предприятий».

16. «Зрение» / под ред. Н.И.Кудряшовой, М., «Машиностроение», 1995.

17. «Временные рекомендации труда операторов за дисплеями». ГОСТ 12.1.006-84.

18. СНиП2963-84 «Нормирование электромагнитных полей».

19. «Современные нормы электростатического и электромагнитного излучения», «Computer World» №7, 1995.3

6. ПРИЛОЖЕНИЕ. ТЕКСТЫ ПРОГРАММ ДЛЯ BORLAND C++ И MATHLAB 4.0 FOR WINDOWS
6.1. ОСНОВНОЙ ПРОГРАММНЫЙ МОДУЛЬ MAIN.CPP
#include <fstream.h>

#include <iostream.h>

#include <conio.h>

#include <stdlib.h>

#include "rk5.h"

#include "sfun.h"

#include "init.h"

#include <math.h>
typedef long double real;
const float g_r = M_PI/180.;

const float r_g = 180./M_PI;
real t_beg;

real t_end;

real dt;

real toler;

int Np;

int Curp;

real dTp;

real mu_z;

real mu_s;

real mu_l;

real m;

real m_t;

real W;

real w_s;

real w_z;

real w_l;

real ww_l;

real xs,ys,zs;

real xl,yl,zl;

real Fz,Fs,Fl,Fa,U20;

real J1,J2,J3;

int nomin;

real par[8];

real parn[8];

real a_p,e_p,p_p,Om_p,i_p,om_p,Rp_p,Ra_p;

real y_main[6];

real prmt[5];

int Fl_u;

real u_last;

int Fl_ka;

int Fl_kp;

int Fl_ki;

int Fl_i;

int Fl_p;

int Fl_a;

int Fl_lu;

int Fl_pkT;

real dl;

real T_vd;

real dRa;

real dRp;

int Sig;

int Sig_a;

real Tkor;

real Tkore;

real Vkor[3];

real akor[3];

int Fl_l0;

int Fl_l1;

int Fl_pki;

real dV_ps;

real dV_as;

real dV_is;

real dV_ss;
ofstream m_y ("m_y.dat");

ofstream m_f ("m_f.dat");

ofstream m_s ("m_s.dat");

ofstream m_l ("m_l.dat");

ofstream m_par ("m_par.dat");

ofstream u_f ("u_f.dat");

ofstream u_par ("u_par.dat");

ofstream k_par ("k_par.dat");
void out_p(real,real *,real*,int,int,real*);
void main()

{

clrscr();

init_m();

real dery[]={ .167, .167, .167, .167, .166, .166};

int ihlf;

int ndim = 6;
Drkgs(prmt,y_main,dery,ndim,ihlf,fct,out_p);

clrscr();

if (ihlf<11)

{

cout << '\n' << "Успешное завершение моделирования" << '\n';

cout << " t0 = " << t_beg << " tk = " << t_end << " dt = " << dt;

cout << '\n' << "Число делений шага=" << ihlf;

}

else

{

cout << '\n' << "Ненормальное завершение моделирования" << '\n';

cout << " t0 = " << t_beg << " tk = " << t_end << " dt = " << dt;

cout << '\n' << "Число делений шага=" << ihlf;

}

getch();
m_y.close();

m_f.close();

m_s.close();

m_l.close();

m_par.close();

u_f.close();

u_par.close();

k_par.close();

}
void out_p(real x,real *y,real*,int,int,real*)

{

if (x >= (dTp*Curp))

{

Curp++;

gotoxy(1,20);

cout << "Процесс выполнения:" << float(Curp)*100./Np << " % " << '\n';
cout.precision(7);
m_y << x << '\t' << y[0] << '\t' << y[1] << '\t' << y[2] << '\t'

      << y[3] << '\t' << y[4] << '\t' << y[5] << '\n';

m_f << x << '\t' << Fz << '\t' << Fs << '\t' << Fl << '\t' << Fa

      << '\t' << U20 << '\n';

m_s << x << '\t' << xs << '\t' << ys << '\t' << zs << '\n';

m_l << x << '\t' << xl << '\t' << yl << '\t' << zl << '\n';

m_par << x << '\t' << par[0] << '\t' << par[1] << '\t' << par[2]

          << '\t' << par[3] << '\t' << par[4] << '\t' << par[5]

          << '\t' << par[6] << '\t' << par[7] << '\n';

}

if (Fl_u && (par[7] > parn[7]))

{

Fl_u = 0;

dl = -(w_z-w_s)*(par[6]-parn[6]);
u_par << x << '\t' << par[0] << '\t' << par[1] << '\t' << par[2]

          << '\t' << par[3] << '\t' << par[4] << '\t' << par[5]

          << '\t' << par[6] << '\t' << par[7] << '\n';

u_f << x << '\t' << Fz << '\t' << Fs << '\t' << Fl

      << '\t' << Fa << '\t' << U20 << '\n';

}

if ((x > 79000) && (x < 81000))

{

k_par << x << '\t' << par[5] << '\t' << par[7] << '\n';

}

}
6.2. ПОДПРОГРАММА РАСЧЕТА ВОЗМУЩАЮЩИХ УСКОРЕНИЙ, ПАРАМЕТРОВ ОРБИТЫ И КОРРЕКЦИИ SFUN.CPP
#include "sfun.h"
const real p = 4.64e-6;

const real sm_s = 8.;

const real A = 1.496e11;

const real Cx = 2.;

const real sm_a = 2.5;

const real ro = 5.098e-13;
void korr (real& t, real *f, real *dery);
void fct(real& t, real *f, real *dery)

{

real x = f[0];

real y = f[1];

real z = f[2];

real Vx = f[3];

real Vy = f[4];

real Vz = f[5];
real Tet_s = (28.1+60*g_r)+w_s*t;

real e_0 = 23.45*g_r;
xs = A*cos(Tet_s);

ys = A*sin(Tet_s)*cos(e_0);

zs = A*sin(Tet_s)*sin(e_0);
real Tet_l = 0+w_l*t;

real Om_l = 0-ww_l*t;

real i_l = acos(cos(e_0)*cos(5.15*g_r)-sin(e_0)*sin(5.15*g_r)*cos(Om_l));

real rsr_l = 3.8448e8;
xl = rsr_l*(cos(Tet_l)*cos(Om_l)-cos(i_l)*sin(Tet_l)*sin(Om_l));

yl = rsr_l*(cos(Tet_l)*sin(Om_l)+cos(i_l)*sin(Tet_l)*cos(Om_l));

zl = rsr_l*sin(i_l)*sin(Tet_l);
real R_ka = sqrt(x*x+y*y+z*z);
real Fz_x = -mu_z*x/pow(R_ka,3.);

real Fz_y = -mu_z*y/pow(R_ka,3.);

real Fz_z = -mu_z*z/pow(R_ka,3.);
real mu_sd = p*sm_s*A*A/m;
real R_s = sqrt((x-xs)*(x-xs)+(y-ys)*(y-ys)+(z-zs)*(z-zs));
real Fs_x = -(mu_s-mu_sd)*x/pow(R_s,3.);

real Fs_y = -(mu_s-mu_sd)*y/pow(R_s,3.);

real Fs_z = -(mu_s-mu_sd)*z/pow(R_s,3.);
real R_l = sqrt((x-xl)*(x-xl)+(y-yl)*(y-yl)+(z-zl)*(z-zl));
real Fl_x = -mu_l*x/pow(R_l,3.);

real Fl_y = -mu_l*y/pow(R_l,3.);

real Fl_z = -mu_l*z/pow(R_l,3.);
real V_ka = sqrt(Vx*Vx+Vy*Vy+Vz*Vz);
real Fa_x = (-Cx*sm_a/(2*m))*ro*V_ka*Vx;

real Fa_y = (-Cx*sm_a/(2*m))*ro*V_ka*Vy;

real Fa_z = (-Cx*sm_a/(2*m))*ro*V_ka*Vz;
const real c20 = -1.09808e-3;

const real c22 = 5.74e-6;

const real d22 = -1.58e-6;

const real r_e = 6378137.;
real cr = mu_z*r_e*r_e/pow(R_ka,5);

real lr = 2*atan(y/x);

real mr = 3*(c22*cos(lr)+d22*sin(lr));
real U20_x = cr*x*(c20*(1.5-7.5*z*z/pow(R_ka,2))+mr*(5*z*z/pow(R_ka,2)-3));

real U20_y = cr*y*(c20*(1.5-7.5*z*z/pow(R_ka,2))+mr*(5*z*z/pow(R_ka,2)-3));

real U20_z = cr*z*(c20*(4.5-7.5*z*z/pow(R_ka,2))+5*mr*(z*z/pow(R_ka,2)-1));
dery[0] = Vx;

dery[1] = Vy;

dery[2] = Vz;

dery[3] = (Fz_x+U20_x+Fs_x+Fl_x+Fa_x+akor[0]);

dery[4] = (Fz_y+U20_y+Fs_y+Fl_y+Fa_y+akor[1]);

dery[5] = (Fz_z+U20_z+Fs_z+Fl_z+Fa_z+akor[2]);
Fz = sqrt(Fz_x*Fz_x+Fz_y*Fz_y+Fz_z*Fz_z);

Fs = sqrt(Fs_x*Fs_x+Fs_y*Fs_y+Fs_z*Fs_z);

Fl = sqrt(Fl_x*Fl_x+Fl_y*Fl_y+Fl_z*Fl_z);

Fa = sqrt(Fa_x*Fa_x+Fa_y*Fa_y+Fa_z*Fa_z);

U20 = sqrt(U20_x*U20_x+U20_y*U20_y+U20_z*U20_z);
parn[3] = parn[3]+w_s*t;
par_or(f,par);
korr(t,f,dery);
if ((u_last-par[7]) > 300*g_r)

Fl_u = 1;
u_last = par[7];

}
void korr(real& t, real *f, real *)

{

if (t > (Tkor+172800.))

{

if ((fabs(dl) > 0.1*g_r) && (!Fl_ka) && (!Fl_kp) && (!Fl_ki))

    {

    Fl_kp = 1;

    Fl_ka = 0;

    Fl_ki = 0;

    cout << "Результат измерений накоплен" << '\n';

    cout << "Необходима коррекция периода. dl=" << dl*r_g << "град." << '\n';

    cout << "Период ном.=" << parn [6] << "Период тек.=" << par[6] << '\n';

    cout << "Параметры орбиты" << '\n';

    cout << " Rp = " << par[2]*(1-par[1]) << '\n';

    cout << " Ra = " << par[2]*(1+par[1]) << '\n';

    cout << " p = " << par[0] << '\n';

    cout << " a = " << par[2] << " e = " << par[1] << "\n T = "

          << par[6] << " w = " << par[5]*r_g << " u = " << par[7]*r_g

          << '\n';

    clrscr();

    }

}
Fl_a = 0;

Fl_p = 0;

Fl_lu = 0;
real da;

if (par[5] > par[7])

da = fabs(par[5]-par[7]-M_PI);

else

da = fabs(par[5]-par[7]+M_PI);
if (da < .1*g_r)

{

Fl_a = 1;

}
if (fabs(par[5] - par[7]) < .1*g_r)

{

Fl_p = 1;

}
if (par[7] < .1*g_r )

{

Fl_lu = 1;

}
real Vk;
if (T_vd)
if (t >= (T_vd +20))

{

T_vd = 0;

akor[0] = 0;

akor[1] = 0;

akor[2] = 0;

cout << "Выкл.дв. \n t = " << t;

}

if (((Fl_kp && Fl_a) || (Fl_ka && Fl_p) || (Fl_ki && Fl_lu)) && (!T_vd))

{

cout << " \n Коррекция \n";

cout << "\n Начало t=" << t << "сек \n";

int sim;
if ((t-Tkor) < 2500)

    {

    cout << "Не корректировать!";

    return;

    }

Tkor = t;
real R_t = sqrt(f[0]*f[0]+f[1]*f[1]+f[2]*f[2]);

real V_t = sqrt(f[3]*f[3]+f[4]*f[4]+f[5]*f[5]);

real R_n = parn[0];
if (Fl_a)

    {

    dRa = R_t-R_n;

    dRp = par[2]*(1-par[1])-R_n;

    cout << "Апоцентр dRp:" << dRp << "м \n";

    cout << "dRa:" << dRa << "м \n";

    cout << "w=" << par[5]*r_g << "u=" << par[7]*r_g << '\n';

    real l,ln;

    l = -(w_z-w_s)*par[6];

    ln = -(w_z-w_s)*parn[6];

    dl = -(w_z-w_s)*(par[6]-parn[6]);

    cout << "T=" << par[6] << "Тном=" << parn[6] << " T-Tном="

          << par[6]-parn[6] << '\n' << "l=" << l*r_g << "lном="

          << ln*r_g << "l-lном=" << (l-ln)*r_g << "dl=" << dl

          << '\n';
    if (dRp > 0)

      Sig_a = -1;

    else

      Sig_a = 1;
    cout << "Знак ускорения:" << Sig_a << '\n';

    clrscr();
    real Rp = par[2]*(1-par[1]);

    real Ra_p = par[2]*(1+par[1]);

    real Rp_p2 = Rp;

    real Ra_p2 = R_t;

    cout << "Rp=" << Rp_p2 << "Ra=" << Ra_p2 << '\n';

    cout << "Ra_p=" << Ra_p << "\n Rt=" << R_t << '\n';
    if (fabs(Rp - R_n) < 500)

      {

      Fl_kp = 0;

      Fl_ka = 1;

      cout << "Закончить коррекцию в апоцентре \n" << "dRp=" << Rp-R_n

             << "dRa=" << dRa << "t=" << t << '\n';

      cout << "Параметры орбиты: \n" << "Rp=" << par[2]*(1-par[1])

             << "Ra=" << par[2]*(1+par[1]) << "\n p=" << par[0]

             << "a=" << par[2] << "e=" << par[1] << "\n T="

             << par[6] << "w=" << par[5]*r_g << "u=" << par[7]*r_g

             << '\n';

      cout << "Суммарный импульс для коррекции перицентра=" << dV_ps << '\n';

      clrscr();

      }

    else

      {

      if (R_t > R_n)

          {

          Rp_p = R_n;

          Ra_p = R_t;

          a_p = (Ra_p+Rp_p)/2.;

          e_p = 1-Rp_p/a_p;

          p_p = a_p*(1-e_p*e_p);
          Vk = sqrt(mu_z/p_p)*(1-e_p);

          }

      else

          {

          Rp_p = R_t;

          Ra_p = R_n;

          a_p = (Ra_p+Rp_p)/2.;

          e_p = 1-Rp_p/a_p;

          p_p = a_p*(1-e_p*e_p);
          Vk = sqrt(mu_z/p_p)*(1+e_p);

          }

      real dV = Vk-V_t;

      real dVmax = 20*25./m;
      cout << "\n dVтреб=" << dV << "dVmax за 20 сек=" << dVmax;
      if (fabs(dV) > dVmax)

          {

          akor[0] = Sig_a*(25./m)*f[3]/V_t;

          akor[1] = Sig_a*(25./m)*f[4]/V_t;

          akor[2] = Sig_a*(25./m)*f[5]/V_t;

          cout << "\n dV=" << dV << "dVmax=" << dVmax;

          cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

               << '\t' << akor[2] << '\t' <<

               sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';
          dV_ps = dV_ps+dVmax;

          cout << "Суммарный импульс=" << dV_ps << '\n';

          }

      else

          {

          akor[0] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[3]/V_t;

          akor[1] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[4]/V_t;

          akor[2] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[5]/V_t;

          cout << "\n dV=" << dV << "dVmax=" << dVmax;

          cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

               << '\t' << akor[2] << '\t' <<

               sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';
          dV_ps = dV_ps+fabs(dV);

          cout << "Суммарный импульс=" << dV_ps << '\n';

          }
      if (dVmax > fabs(dV))

          {

          dVmax = fabs(dV);

          real Vk_r = Sig_a*dVmax+V_t;

          real Ra_r = R_t;

          real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;

          real a_r = Ra_r/(1+e_r);

          real p_r = a_r*(1-e_r*e_r);

          real Rp_r = a_r*(1-e_r);

          cout << "Параметры орбиты: \n" << " Rp_r = " << Rp_r

               << " Ra_r = " << Ra_r << "\n p_r = " << p_r << " a_r = "

               << a_r << " e_r = " << e_r << '\n';

          }

      else

          {

          real Vk_r = Sig_a*dVmax+V_t;

          real Ra_r = R_t;

          real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;

          real a_r = Ra_r/(1+e_r);

          real p_r = a_r*(1-e_r*e_r);

          real Rp_r = a_r*(1-e_r);

          cout << "Параметры орбиты: \n" << " Rp_r = " << Rp_r

               << " Ra_r = " << Ra_r << "\n p_r = " << p_r << " a_r = "

               << a_r << " e_r = " << e_r << '\n';

          }

      T_vd = t;

      cout << "Вкл.дв. t=" << T_vd << '\n';

      }

    }
if (Fl_p)

    {

    dRp = R_t-R_n;

    dRa = par[2]*(1+par[1])-R_n;

    cout << " Перицентра - dRp:" << dRp << "м \n";

    cout << "dRa:" << dRa << "м. \n";

    cout << "w=" << par[5]*r_g << "u=" << par[7]*r_g << '\n';

    real l,ln;

    l = -(w_z-w_s)*par[6];

    ln = -(w_z-w_s)*parn[6];

    dl = -(w_z-w_s)*(par[6]-parn[6]);

    cout << "T=" << par[6] << "Tном=" << parn[6] << "T-Tном="

          << par[6]-parn[6] << '\n' << "l=" << l*r_g << "lном="

          << ln*r_g << "l-lном=" << (l-ln)*r_g << "dl=" << dl << '\n';
    if (dRa > 0)

      Sig_a = -1;

    else

      Sig_a = 1;
    cout << "Знак ускорения:" << Sig_a << '\n';

    clrscr();
    real Ra = par[2]*(1+par[1]);

    real Rp_p1 = R_t;

    real Ra_p1 = Ra;

    cout << "Rp=" << Rp_p1 << "Ra=" << Ra_p1 << '\n';
    if ((fabs(Ra-R_n) < 500) || (fabs(dl*r_g) < .0001))

      {

      cout << "Закончить коррекцию в перицентре \n" << "dRa="

             << Ra-R_n << "dRp=" << dRp << "t=" << t << '\n';

      cout << "Параметры орбиты: \n " << "Rp="

             << par[2]*(1-par[1]) << "Ra=" << par[2]*(1+par[1])

             << " \n p=" << par[0] << "a=" << par[2] << "e="

             << par[1] << " \n T=" << par[6] << "w=" << par[5]*r_g

             << "u=" << par[7]*r_g << '\n';
      cout << "Суммарный импульс для коррекции перицентра=" << dV_as << '\n';
      clrscr();
      Fl_ka = 0;

      Fl_ki = 1;

      }

    else

      {

      if (R_t > R_n)

          {

          Rp_p = R_n;

          Ra_p = R_t;

          a_p = (Ra_p+Rp_p)/2.;

          e_p = 1-Rp_p/a_p;

          p_p = a_p*(1-e_p*e_p);
          Vk = sqrt(mu_z/p_p)*(1-e_p);

          }

      else

          {

          Rp_p = R_t;

          Ra_p = R_n;

          a_p = (Ra_p+Rp_p)/2.;

          e_p = 1-Rp_p/a_p;

          p_p = a_p*(1-e_p*e_p);
          Vk = sqrt(mu_z/p_p)*(1+e_p);

          }

      real dV = Vk-V_t;

      real dVmax = 20*25./m;
      cout << "\n dVнадо=" << dV << " dVmax за 20 сек=" << dVmax;
      if (fabs(dV) > dVmax)

          {

          akor[0] = Sig_a*(25./m)*f[3]/V_t;

          akor[1] = Sig_a*(25./m)*f[4]/V_t;

          akor[2] = Sig_a*(25./m)*f[5]/V_t;

          cout << "\n dV=" << dV << "dVmax=" << dVmax;

          cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

               << '\t' << akor[2] << '\t' <<

               sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';
          dV_as = dV_as+dVmax;

          cout << "Суммарный импульс=" << dV_as << '\n';

          }

      else

          {

          akor[0] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[3]/V_t;

          akor[1] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[4]/V_t;

          akor[2] = Sig_a*(fabs(dV)/dVmax)*(25./m)*f[5]/V_t;

          cout << "\n dV=" << dV << " dVmax=" << dVmax;

          cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

               << '\t' << akor[2] << '\t' <<

               sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';
          dV_as = dV_as+fabs(dV);

          cout << "Суммарный импульс=" << dV_as << '\n';

          }
      if (dVmax > fabs(dV))

          {

          dVmax = fabs(dV);

          real Vk_r = Sig_a*dVmax+V_t;

          real Ra_r = R_t;

          real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;

          real a_r = Ra_r/(1+e_r);

          real p_r = a_r*(1-e_r*e_r);

          real Rp_r = a_r*(1-e_r);

          cout << "Параметры орбиты: \n" << "Rp_r=" << Rp_r

               << "Ra_r=" << Ra_r << "\n p_r=" << p_r << "a_r="

               << a_r << "e_r=" << e_r << '\n';

          }

      else

          {

          real Vk_r = Sig_a*dVmax+V_t;

          real Ra_r = R_t;

          real e_r = -(Vk_r*Vk_r*Ra_r/mu_z)+1;

          real a_r = Ra_r/(1+e_r);

          real p_r = a_r*(1-e_r*e_r);

          real Rp_r = a_r*(1-e_r);

          cout << "Параметры орбиты: \n" << "Rp_r=" << Rp_r

               << "Ra_r=" << Ra_r << "\n p_r=" << p_r << "a_r="

               << a_r << "e_r=" << e_r << '\n';

          }

      T_vd = t;

      cout << "Вкл.дв. t=" << T_vd << '\n';

      }

    }
if (Fl_lu)

    {

    real di = par[4]-parn[4];

    cout << "Линия узлов - di: " << di*r_g << "град \n";

    cout << "w=" << par[5]*r_g << "u=" << par[7]*r_g << '\n';

    real l,ln;

    l = -(w_z-w_s)*par[6];

    ln = -(w_z-w_s)*parn[6];

    dl = -(w_z-w_s)*(par[6]-parn[6]);

  cout << "T=" << par[6] << "Tном=" << parn[6] << "T-Tном="

          << par[6]-parn[6] << '\n' << "l=" << l*r_g << "lном="

          << ln*r_g << "l-lном=" << (l-ln)*r_g << "dl=" << dl

          << "\n i=" << par[4]*r_g << "iном=" << parn[4]*r_g << '\n';

    cout << "Параметры орбиты: \n " << "Rp="

          << par[2]*(1-par[1]) << "Ra=" << par[2]*(1+par[1])

          << " \n p=" << par[0] << "a=" << par[2] << "e="

          << par[1] << " \n T=" << par[6] << "w=" << par[5]*r_g

          << "u=" << par[7]*r_g << " \n i=" << par[4]*r_g << '\n';

    clrscr();
  real Vk_x,Vk_y,Vk_z;

    if (fabs(di) < .0001*g_r)

      {

      Fl_ki = 0;

      cout << "Закончить коррекцию наклонения \n "

             << "di=" << (par[4]-parn[4])*r_g << "t=" << t << '\n';

      cout << "Параметры орбиты: \n " << "Rp="

             << par[2]*(1-par[1]) << "Ra=" << par[2]*(1+par[1])

             << " \n p=" << par[0] << "a=" << par[2] << "e="

             << par[1] << " \n T=" << par[6] << "w=" << par[5]*r_g

             << "u=" << par[7]*r_g << " \n i=" << par[4]*r_g << '\n';

      cout << "Суммарный импульс=" << dV_is

             << '\n';
      clrscr();

      }

    else

      {

      real teta;

      if (par[7] > par[5])

          teta = 2*M_PI+par[7]-par[5];

      else

          teta = par[7]-par[5];
      real Vr_i = sqrt(mu_z/par[0])*par[1]*sin(teta);

      real Vn_i = sqrt(mu_z/par[0])*(1+par[1]*cos(teta));

      V_t = sqrt(f[3]*f[3]+f[4]*f[4]+f[5]*f[5]);
      Vk_x = -Vn_i*cos(parn[4])*sin(par[3])+Vr_i*cos(par[3]);

      Vk_y = Vn_i*cos(parn[4])*cos(par[3])+Vr_i*sin(par[3]);

      Vk_z = Vn_i*sin(parn[4]);

      Vk = sqrt(Vk_x*Vk_x+Vk_y*Vk_y+Vk_z*Vk_z);
      real dV_x = Vk_x-f[3];

      real dV_y = Vk_y-f[4];

      real dV_z = Vk_z-f[5];

      real dV = sqrt(dV_x*dV_x+dV_y*dV_y+dV_z*dV_z);
      real dVmax = 20*25./m;

      cout << "Vнач=" << V_t << "Vк=" << Vk << "teta=" << teta*r_g

             << '\n';

      cout << "dV=" << dV << "dVmax за 20 сек=" << dVmax;
      if (dV > dVmax)

          {

          akor[0] = (25./m)*dV_x/dV;

          akor[1] = (25./m)*dV_y/dV;

          akor[2] = (25./m)*dV_z/dV;

          cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1] <<

                '\t' << akor[2] << '\t' <<

                sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';
          dV_is = dV_is+dVmax;

          cout << "Суммарный импульс=" << dV_is << '\n';

          }

      else

          {

          akor[0] = (fabs(dV)/dVmax)*(25./m)*dV_x/dV;

          akor[1] = (fabs(dV)/dVmax)*(25./m)*dV_y/dV;

          akor[2] = (fabs(dV)/dVmax)*(25./m)*dV_z/dV;

          cout << "\n Корректирующее ускорение:" << akor[0] << '\t' << akor[1]

               << '\t' << akor[2] << '\t'<<

               sqrt(akor[0]*akor[0]+akor[1]*akor[1]+akor[2]*akor[2]) << '\n';
          dV_is = dV_is+fabs(dV);

          cout << "Суммарный импульс=" << dV_is << '\n';

          }
      T_vd = t;

      cout << "Вкл.дв. t=" << T_vd << '\n';

      }

    }
if ((!Fl_ka) && (!Fl_kp) && (!Fl_ki))

    {

    cout << "Коррекция окончена!" << '\n';
    real m_t;

    dV_ss = dV_ps+dV_as+dV_is;

    m_t = m*(1-exp(-dV_ss/W));

    cout << "Потребный импульс: \n - перицентра dV_ps="

          << dV_ps << "\n апоцентра dV_as=" << dV_as

          << "\n Суммарный импульс=" << dV_ss << "Масса топлива=" << m_t

          << '\n';
    dV_ps = 0;

    dV_as = 0;

    dV_is = 0;

    dV_ss = 0;

    m_t = 0;

    }

}

}
void par_or(real *f, real *par)

{

real x = f[0];

real y = f[1];

real z = f[2];

real Vx = f[3];

real Vy = f[4];

real Vz = f[5];
real c1 = (y*Vz-z*Vy);

real c2 = (z*Vx-x*Vz);

real c3 = (x*Vy-y*Vx);

real C = sqrt(c1*c1+c2*c2+c3*c3);
par[0] = (C/mu_z)*C;
real R_ka = sqrt(x*x+y*y+z*z);

real V_ka = sqrt(Vx*Vx+Vy*Vy+Vz*Vz);
real f1 = (Vy*c3-Vz*c2)-(mu_z*x/R_ka);

real f2 = (Vz*c1-Vx*c3)-(mu_z*y/R_ka);

real f3 = (Vx*c2-Vy*c1)-(mu_z*z/R_ka);

real F = sqrt(f1*f1+f2*f2+f3*f3);
real h = V_ka*V_ka-(2*mu_z/R_ka);
if ((1+h*C*C/(mu_z*mu_z)) < 0)

{

cout << " Ошибка! \n";

getch();

}
par[1] = F/mu_z;
if ((1-par[1]*par[1]) < 0)

{

cout << " (1-e*e) < 0 Ошибка! \n";

getch();

}
par[2] = par[0]/(1-par[1]*par[1]);
par[4] = acos(c3/C);
real s_Om = c1/(C*sin(par[4]));

real c_Om = -c2/(C*sin(par[4]));
if (s_Om >= 0)

par[3] = acos(c_Om);

else

par[3] = 2*M_PI-acos(c_Om);
real c_om = (f1*cos(par[3])+f2*sin(par[3]))/F;

real s_om = f3/(F*sin(par[4]));
if (s_om > 0)

par[5] = acos(c_om);

else

par[5] = 2*M_PI - acos(c_om);
if (par[2] < 0)

{

cout << " Ошибка! \n";

getch();

}
par[6] = 2*M_PI*sqrt((par[2]/mu_z)*par[2]*par[2]);
real c_u = (x*cos(par[3])+y*sin(par[3]))/R_ka;

real s_u = z/(R_ka*sin(par[4]));
if (s_u > 0)

par[7] = acos(c_u);

else

par[7] = 2*M_PI - acos(c_u);

}
#include "rk5.h"

#include <iostream.h>
void Drkgs(real *prmt,real *y,real *dery,int ndim,int& ihlf,

    void (*fct)(real &,real*,real*),

    void (*out_p)(real,real*,real*,int,int,real*))

{

static real a[] = { 0.5, 0.292893218811345248, 1.70710678118665475,

                       0.16666666666666667 };

static real b[] = { 2.0, 1.0, 1.0, 2.0 };

static real c[] = { 0.5, 0.292893218811345248, 1.70710678118665475, 0.5 };

real *aux[8];

int i,j,imod,itest,irec,istep,iend;

real delt,aj,bj,cj,r,r1,r2,x,xend,h;

for (i=0; i<8; i++) aux[i] = new real[ndim];

for (i=0; i<ndim; i++) aux[7][i] = (1./15.)*dery[i];

x = prmt[0];

xend = prmt[1];

h = prmt[2];

prmt[4] = 0.0;
fct(x,y,dery);

r = h*(xend-x);

if (r <= 0.0)

{

ihlf = 13;

if (r == 0.0) ihlf = 12;

goto l39;

}

for(i=0; i<ndim; i++)

{

aux[0][i] = y[i];

aux[1][i] = dery[i];

aux[2][i] = 0.0;

aux[5][i] = 0.0;

}

irec = 0;

h = h+h;

ihlf = -1;

istep = 0;

iend = 0;

l4: r = (x+h-xend)*h;

if (r >= 0.0)

{

iend = 1;

if (r > 0.0) h = xend-x;

}

out_p(x,y,dery,irec,ndim,prmt);

if (prmt[4] != 0.0) goto l40;

itest = 0;

l9: istep++;

j = 0;

l10: aj = a[j];

bj =b[j];

cj = c[j];

for (i=0; i<ndim; i++)

{

r1 = h*dery[i];

r2 = aj*(r1-bj*aux[5][i]);

y[i] = y[i]+r2;

r2 = r2+r2+r2;

aux[5][i] += r2-cj*r1;

}

if (j-3 < 0)

{

j++;

if (j-2 != 0) x = x+0.5*h;

fct(x,y,dery);

goto l10;

}

if (itest <= 0)

{

for (i=0; i<ndim; i++) aux[3][i] = y[i];

itest = 1;

istep = istep+istep-2;

l18: ihlf++;

x = x-h;

h = 0.5*h;

for (i=0; i<ndim; i++)

    {

    y[i] = aux[0][i];

    dery[i] = aux[1][i];

    aux[5][i] = aux[2][i];

    }

goto l9;

}

imod = istep/2;

if (istep-imod-imod != 0)

{

fct(x,y,dery);

for (i=0; i<ndim; i++)

    {

   aux[4][i] = y[i];

    aux[6][i] = dery[i];

    }

goto l9;

}

delt = 0.0;

for (i=0; i<ndim; i++)

delt += aux[7][i]*fabs(aux[3][i]-y[i]);

if (delt-prmt[3] > 0.0)

{

if (ihlf-10 >= 0)

    {

    ihlf = 11;

    fct(x,y,dery);

    goto l39;

    }

for (i=0; i<ndim; i++) aux[3][i] = aux[4][i];

istep = istep+istep-4;

x = x-h;

iend = 0;

goto l18;

}

fct(x,y,dery);

for (i=0; i<ndim; i++)

{

aux[0][i] = y[i];

aux[1][i] = dery[i];

aux[2][i] = aux[5][i];

y[i] = aux[4][i];

dery[i] = aux[6][i];

}

out_p(x-h,y,dery,ihlf,ndim,prmt);

if (prmt[4] != 0) goto l40;

for (i=0; i<ndim; i++)

{

y[i] = aux[0][i];

dery[i] = aux[1][i];

}

irec = ihlf;

if (iend > 0) goto l39;

ihlf--;

istep = istep/2;

h = h+h;

if (ihlf < 0) goto l4;

imod = istep/2;

if ((istep-2*imod != 0) || (delt-0.02*prmt[3] > 0.0)) goto l4;

ihlf--;

istep = istep/2;

h = h+h;

goto l4;

l39: out_p(x,y,dery,ihlf,ndim,prmt);

l40: for (i=0; i<ndim; i++) delete aux[i];

return;

}
6.3. ФАЙЛ НАЧАЛЬНОЙ ИНИЦИАЛИЗАЦИИ INIT.H
ifndef _INIT

#define _INIT
#include "def.h"

#include <stdlib.h>

#include <fstream.h>
ifstream if_init;
void nex_ln (void);
void init_m()

{

Np = 150;

t_beg = 0;

t_end = 8000000;

dt = 2;

toler = .05;
dTp = (t_end-t_beg)/float(Np);

Curp = 0;
J1 = 532;

J2 = 563;

J3 = 697;

m = 597.;

W = 2200;
mu_z = 3.9858e14;

mu_s = 1.3249e20;

mu_l = 4.9027e12;
w_s = 2*M_PI/(365.2422*24*3600);

w_z = 2*M_PI/(24*3600);

w_l = 2*M_PI/(27.32*24*3600);

ww_l = 2*M_PI/(18.6*365.2422*24*3600);
parn[0] = 6952137.;

parn[1] = 0;

parn[2] = 6952137;

parn[3] = 28.1*g_r;

parn[4] = 97.6*g_r;

parn[5] = 63.1968*g_r;

parn[6] = 5769.;

parn[7] = 5.751*g_r;
Fl_u = 1;

u_last = parn[7];
Fl_ka = 0;

Fl_kp = 0;

Fl_ki = 0;

Fl_p = 0;

Fl_a = 0;

Fl_i = 0;

Fl_pkT = 0;

Tkor = 0;

T_vd = 0;

akor[0] = 0;

akor[1] = 0;

akor[2] = 0;

dV_ps = 0;

dV_as = 0;

dV_is = 0;

dV_ss = 0;

Fl_l0 = 0;

Fl_l1 = 0;

Fl_pki = 0;
real x0 = 6137262.9+7000;

real y0 = 3171846.1+7000;

real z0 = 689506.95+7000;

real Vx0 = -201.288+5;

real Vy0 = -1247.027+5;

real Vz0 = 7472.65+5;
prmt[0] = t_beg;

prmt[1] = t_end;

prmt[2] = dt;

prmt[3] = toler;

prmt[4] = 0.0;
y_main[0] = x0;

y_main[1] = y0;

y_main[2] = z0;

y_main[3] = Vx0;

y_main[4] = Vy0;

y_main[5] = Vz0;

}
void nex_ln (void)

{

char ch;

if_init.get(ch);

while (ch != '\n')

if_init.get(ch);

}

#endif
6.4 ФАЙЛ ОПИСАНИЯ ПЕРЕМЕННЫХ DEF.H
#ifndef _DEFH

#define _DEFH
#include <math.h>
typedef long double real;
extern const float g_r;

extern const float r_g;
extern int Np;

extern int Curp;

extern real dTp;

extern real t_beg;

extern real t_end;

extern real dt;

extern real toler;

extern real J1,J2,J3;

extern real mu_z;

extern real mu_s;

extern real mu_l;

extern real m;

extern real m_t;

extern real W;

extern real w_s;

extern real w_z;

extern real w_l;

extern real ww_l;

extern real xs,ys,zs;

extern real xl,yl,zl;

extern real Fz,Fs,Fl,Fa,U20;

extern int nomin;

extern real par[8];

extern real parn[8];

extern real a_p,e_p,p_p,Om_p,i_p,om_p,Rp_p,Ra_p;

extern real y_main[6];

extern real prmt[5];

extern int Fl_u;

extern real u_last;

extern int Fl_ka;

extern int Fl_kp;

extern int Fl_ki;

extern int Fl_i;

extern int Fl_p;

extern int Fl_a;

extern int Fl_lu;

extern int Fl_pkT;

extern real dl;

extern real T_vd;

extern real dRa;

extern real dRp;

extern int Sig;

extern int Sig_a;

extern real Vkor[3];

extern real akor[3];

extern real Tkor;

extern real Tkore;

extern real dV_ps;

extern real dV_as;

extern real dV_is;

extern real dV_ss;

extern int Fl_l0;

extern int Fl_l1;

extern int Fl_pki;
#endif
6.5 ФАЙЛ SFUN.H
#ifndef _SFUN

#define _SFUN
#include "def.h"

#include <iostream.h>

#include <conio.h>

#include <math.h>
void out_p(real x,real *y,real*,int,int,real *);

real interpl(real*,real*,int,real);

void fct(real& ,real *y,real *dery);

void par_or(real *,real *);
#endif
6.5 ФАЙЛ RK5.H
#ifndef _RK5

#define _RK5

#include "def.h"

#include <iostream.h>

#include <conio.h>

#include "sfun.h"
void Drkgs(real *prmt,real *y,real *dery,int ndim,int& ihlf,

void (*fct)(real&,real*,real*),

void (*out_p)(real,real*,real*,int,int,real*));
#endif
6.6 ПРОГРАММА ПОСТРОЕНИЯ ВРЕМЕННЫХ ДИАГРАММ
clc

g_r = pi/180;

r_g = 180/pi;

load m_y.dat

t = m_y(:,1);

x = m_y(:,2);

y = m_y(:,3);

z = m_y(:,4);

Vx = m_y(:,5);

Vy = m_y(:,6);

Vz = m_y(:,7);

clear m_y;

s_tmp = size(t);

s_m = s_tmp(1);

clear s_tmp;

load m_f.dat

Fz = m_f(:,2);

Fs = m_f(:,3);

Fl = m_f(:,4);

Fa = m_f(:,5);

U20 = m_f(:,6);

clear m_f;

load m_s.dat

xs = m_s(:,2);

ys = m_s(:,3);

zs = m_s(:,4);

clear m_s;

load m_par.dat

p = m_par(:,2);

e = m_par(:,3);

a = m_par(:,4);

Om = m_par(:,5);

i = m_par(:,6);

omg = m_par(:,7);

T = m_par(:,8);

u = m_par(:,9);

clear m_par;

p_n = 6952137.;

e_n = 0;

a_n = 6952137.;

Om_n0 = 28.1*g_r;

i_n = 97.6*g_r;

omg_n = 346.725*g_r;

T_n = 5765;

ws = 2*pi/(365.2422*24*3600);

for j = 1:s_m, tmp(j) = Om_n0+ws*t(j);

end

Om_n = tmp';

clear tmp;

map = [1,1,1];

colormap(map);

plot(t,p,'y-',[min(t) max(t)],[p_n p_n],'r-'), title (' Фокальный параметр '), grid on;

print -dwin;

pause;

plot(t,p-p_n,'y-'), title (' dp '), grid on;

print -dwin;

pause;

plot(t,e,'y-',[min(t) max(t)],[e_n e_n],'r-'), title (' Эксцентриситет '), grid on;

print -dwin;

pause;

plot(t,e-e_n,'y-'), title (' de '), grid on;

print -dwin;

pause;

plot(t,a,'y-',[min(t) max(t)],[a_n a_n],'r-'), title (' Большая полуось орбиты '), grid on;

print -dwin;

pause;

plot(t,a-a_n,'y-'), title (' da '), grid on;

print -dwin;

pause;

plot(t,Om*r_g,'y-',t,Om_n*r_g,'r-'), title (' Долгота восходящего узла '), grid on;

print -dwin;

pause;

plot(t,Om*r_g-Om_n*r_g,'y-'), title (' dOm '), grid on;

print -dwin;

pause;

plot(t,i*r_g,'y-',[min(t) max(t)],[i_n*r_g i_n*r_g],'r-'), title (' Наклонение '), grid on;

print -dwin;

pause;

plot(t,i*r_g-i_n*r_g,'y-'), title (' di '), grid on;

print -dwin;

pause;

plot(t,T,'y-',[min(t) max(t)],[T_n T_n], 'r-'), title (' Период '), grid on;

print -dwin;

pause;

plot(t,T-T_n,'y-'), title (' dT '), grid on;

print -dwin;

pause;

plot3(x,y,z,'b')

axis([min(x) max(x) min(y) max(y) min(z) max(z)])

set(gca,'box','on')

title (' Положение МКА ')

hold on

plt = plot3(0,0,0,'.','erasemode','xor','markersize',24);

dk = ceil(length(y)/2500);

for k = 1:dk:length(y)

set(plt,'xdata',x(k),'ydata',y(k),'zdata',z(k))

drawnow

end

hold off

pause;

plot(t,Fz,'y-'), title (' Гравитация Земли ' ), grid on;

print -dwin;

pause;

plot(t,Fs,'y-'), title (' Гравитация Солнца и солнечное давление '), grid on;

print -dwin;

pause;

plot(t,Fl,'y-'), title (' Гравитация Луны '), grid on;

print -dwin;

pause;

plot(t,Fa,'y-'), title (' Сопротивление атмосферы '), grid on;

print -dwin;

pause;

plot(t,U20,'y-'), title (' Нецентральность гравитационного поля Земли '), grid on;

print -dwin;

pause;

plot(t,Fz+Fs+Fl+Fa+U20,'y-'), title (' Суммарное возмущающее ускорение '), grid on;

print -dwin;

pause;

clear all
clc

g_r = pi/180;

r_g = 180/pi;

p_n = 6952137.;

e_n = 0;

a_n = 6952137.;

Om_n0 = 28.1*g_r;

i_n = 97.6*g_r;

omg_n = 346.725*g_r;

T_n = 5765;

load u_par.dat

t_u = u_par(:,1);

p_u = u_par(:,2);

e_u = u_par(:,3);

a_u = u_par(:,4);

Om_u = u_par(:,5);

i_u = u_par(:,6);

omg_u = u_par(:,7);

T_u = u_par(:,8);

u_u = u_par(:,9);

clear u_par;

load u_f.dat;

Fz_u = u_f(:,2);

Fs_u = u_f(:,3);

Fl_u = u_f(:,4);

Fa_u = u_f(:,5);

U20_u = u_f(:,6);

clear u_f;

s_tmp = size(t_u);

s_m_u = s_tmp(1);

clear s_tmp;

ws = 2*pi/(365.2422*24*3600);

for j = 1:s_m_u, tmp(j) = Om_n0+ws*t_u(j);

end

Om_n_u = tmp';

clear tmp;

plot(t_u,p_u,'y-',[min(t_u) max(t_u)],[p_n p_n],'r-'), title (' Фокальный параметр '), grid on;

print -dwin;

pause;

plot(t_u,p_u-p_n,'y-'), title (' dp '), grid on;

print -dwin;

pause;

plot(t_u,e_u,'y-',[min(t_u) max(t_u)],[e_n e_n],'r-'), title (' Эксцентриситет '), grid on;

print -dwin;

pause;

plot(t_u,e_u-e_n,'y-'), title (' de '), grid on;

print -dwin;

pause;

plot(t_u,a_u,'y-',[min(t_u) max(t_u)],[a_n a_n],'r-'), title (' Большая полуось орбиты '), grid on;

print -dwin;

pause;

plot(t_u,a_u-a_n,'y-'), title (' da '), grid on;

print -dwin;

pause;

plot(t_u,Om_u*r_g,'y-',t_u,Om_n_u*r_g,'r-'), title (' Долгота восходящего узла '), grid on;

print -dwin;

pause;

plot(t_u,Om_u*r_g-Om_n_u*r_g,'y-'), title (' dOm '), grid on;

print -dwin;

pause;

plot(t_u,i_u*r_g,'y-',[min(t_u) max(t_u)],[i_n*r_g i_n*r_g],'r-'), title (' Наклонение '), grid on;

print -dwin;

pause;

plot(t_u,i_u*r_g-i_n*r_g,'y-'), title (' di '), grid on;

print -dwin;

pause;

plot(t_u,T_u,'y-',[min(t_u) max(t_u)],[T_n T_n], 'r-'), title (' Период '), grid on;

print -dwin;

pause;

plot(t_u,T_u-T_n,'y-'), title (' dT '), grid on;

print -dwin;

pause;

clear all

Bukvasha