Здравствуйте, уважаемые хабровчане! Расскажу о своём небольшом опыте работы в проекте, который был посвящён астрономии, и о математике, с которой пришлось повозиться.
Дело было так - учился я в универе, изучал электронику, и после окончания курса по азам С++ , мне предложили подработку в одном проекте, где нужно было написать небольшую программу для моделирования рассеивания света звёздной пылью. Пришёл я на собеседование, которое проводили два доктора астрофизика, поспрашивали про математику и астрономию, поняли что не математик и уж тем более не астроном, но на работу взяли.
Первым моим заданием было подтянуть математику и боль-мень освоить фортран, так как на фортране был софт одного известного учёного М.И.Мищенко (Michael I. Mishchenko). Этот софт определял оптические свойства пыли в зависимости от её концентрации, формы пыли, длины волны и индекса преломления света.
Далее мне нужно было написать свою программу, которая моделирует рассеивание света. Модель простая - звезда, вокруг неё пылевое облако (параметры оптических свойств пылевого облака можно было рассчитать программой М.И.Мищенко), вдали от облака находится оптический телескоп. Моя программа моделирует полученное изображение на CCD матрице телескопа, а именно значение четырёх векторов Стокса, для каждого пикселя CCD матрицы. Метод, который использовался для моделирования рассеивания света - Монте-Карло. Я уговорил начальника - доктора астрофизики использовать C++, т.к. я не хотел писать её на фортране.
Задача была не сложной, начальник выводил формулы и уравнения, я решал их численными методами и ваял программу.
Упрощённая блок-схема программы, к которой я пришёл интуитивно, а после сравнил с одной публикацией (Three Monte Carlo programs of polarized light transport into scattering media: part I (Jessica C. Ramella-Roman, Scott A. Prahl, and Steve L. Jacques) ) о разных методах Монте-Карло для рассеивания света. И понял что моя совпадает с классическим методом Монте-Карло для поляризованного излучения.
![Блок-схема моей программы Блок-схема моей программы](https://habrastorage.org/getpro/habr/upload_files/87a/514/91b/87a51491b7aebe131d8d7f3ad6c6b7ed.png)
Из звезды вылетают пакеты фотонов (каждый пакет фотонов состоит из 1000 фотонов), которые летят в случайном направлении, после часть энергии рассеивается, в зависимости от направления пакета фотонов и оптических свойств пыли рассчитывается дальнейшее направление пакета фотонов. Каждое рассеивание фиксируется на CCD матрицу телескопа. И так много миллионов раз. Чем больше пакетов фотонов рассеяно, тем точнее изображение.
Часть математики я вывел сам - это моделирование излучения пакетов фотонов из сферического источника света и моделирование проекции пылевого облака на плоскость телескопа.
Моделирование сферического источника света:
![Звезда и испускание пакета фотонов из точки её поверхности Звезда и испускание пакета фотонов из точки её поверхности](https://habrastorage.org/getpro/habr/upload_files/6ba/eff/cde/6baeffcde4c484a52d8555afa69ec0c3.png)
- Декартовы координаты центра звезды.
- сферические координаты точки, на поверхности звезды, из которой будет испускаться пакет фотонов.
- угол между нормалью к поверхности звезды и вектором эмиссии (испускания) пакета фотонов.
- полярный и азимутальный углы эмиссии пакета фотонов.
Программа генерирует углы и
, радиус звезды
- константа. Декартовы координаты точки эмиссии пакета фотонов рассчитываются по формуле (1.1)
![](https://habrastorage.org/getpro/habr/upload_files/804/47c/36b/80447c36be784e50ce6c8f5ba25c9555.png)
Единичный вектор рассчитывается по формуле (1.2)
![](https://habrastorage.org/getpro/habr/upload_files/87e/9de/d01/87e9ded01fc2d5399bea1a19ac0fdc66.png)
Угол между нормалью к поверхности звезды и вектором эмиссии
.
Программа рассчитывает вектор эмиссии по формуле (1.3).
![](https://habrastorage.org/getpro/habr/upload_files/698/6a6/2c1/6986a62c1c96df59e7f7ee49059ed36c.png)
Поворачиваем вектор эмиссии вокруг вектора нормали к поверхности звезды
на угол
.
Используем матрицу поворота , чтобы рассчитать вектор эмиссии
после поворота.
![](https://habrastorage.org/getpro/habr/upload_files/265/7df/440/2657df44066be39028f985a70574c0dc.png)
![](https://habrastorage.org/getpro/habr/upload_files/22f/c54/968/22fc54968e35857b6d6db2590ee43ade.png)
![](https://habrastorage.org/getpro/habr/upload_files/2d4/cc2/cdb/2d4cc2cdb4f99b067d881104b3a3bc47.png)
Когда вектор получен,
рассчитываются полярный и азимутальный углы эмиссии
![](https://habrastorage.org/getpro/habr/upload_files/bd6/49c/cce/bd649ccce7de11658c52a5ffa5ddf74e.png)
В результате мы получили пакет фотонов, который испускается из точки поверхности звезды с координатами и летит в направлении которое задают полярный и азимутальный углы
.
-
Моделирование проекции пылевого облака на плоскость телескопа:
Сферическое облако и позиция телескопа до поворотов в программно заданную точку.
Центр системы координат , центр сферического облака
, радиус сферического облака 39. (Это описание математической модели, в программе соответственно используются астрономические расстояния).
Координаты телескопа до поворотов:
![](https://habrastorage.org/getpro/habr/upload_files/7ca/f12/aab/7caf12aab725c5d242bc96d247f26f0c.png)
Первый поворот осуществляется вокруг оси которая параллельна Y оси ,
Телескоп поворачивается на угол (полярный угол позиции телескопа):
![](https://habrastorage.org/getpro/habr/upload_files/bd8/7cf/20c/bd87cf20c5658ecc081d84e8c4848fd8.png)
![](https://habrastorage.org/getpro/habr/upload_files/6bf/005/444/6bf00544448998d4823522ff19608970.png)
![](https://habrastorage.org/getpro/habr/upload_files/3a5/a06/91d/3a5a0691da17beff58eee26dda2032a3.png)
Первый поворот телескопа на угол совершён.
![Позиция телескопа после первого поворота Позиция телескопа после первого поворота](https://habrastorage.org/getpro/habr/upload_files/1e5/e8a/4dd/1e5e8a4ddafebc9ce15c6256aea4e7ad.png)
Далее будет произведён поворот вокруг оси параллельной Z оси , на угол
(азимутальный угол позиции телескопа).
![](https://habrastorage.org/getpro/habr/upload_files/72e/50b/746/72e50b74621a2a099792fd9a1fce66e0.png)
![](https://habrastorage.org/getpro/habr/upload_files/646/ee2/baa/646ee2baafb3762d0803e6ef50cabb97.png)
![](https://habrastorage.org/getpro/habr/upload_files/e0b/f3a/6ca/e0bf3a6ca046fe6fe38ae8d10102de93.png)
Второй поворот телескопа на угол совершён.
![Позиция телескопа после второго поворота Позиция телескопа после второго поворота](https://habrastorage.org/getpro/habr/upload_files/c65/011/d5e/c65011d5e2666baf8fdb9b256bbc64ba.png)
Далее необходимо рассчитать единичный вектор для плоскости телескопа (для проекции облака на плоскость телескопа):
![](https://habrastorage.org/getpro/habr/upload_files/0d5/6d8/f2f/0d56d8f2f3c04ef7a2375f8e9510a2dc.png)
Необходимо рассчитать константу для уравнения плоскости телескопа:
![](https://habrastorage.org/getpro/habr/upload_files/56c/809/517/56c809517e91db41d372821b71aed498.png)
Уравнение плоскости телескопа:
![](https://habrastorage.org/getpro/habr/upload_files/e5a/4e9/362/e5a4e9362ed0251646be23f83ae4f6e6.png)
Необходимо рассчитать рассчитать расстояния от точки облака до плоскости телескопа, - точки облака
![](https://habrastorage.org/getpro/habr/upload_files/5be/64e/ef4/5be64eef411b375c99c83798e7cb116a.png)
Далее получаем точки проекции облака на плоскость телескопа:
![](https://habrastorage.org/getpro/habr/upload_files/447/aa7/d1e/447aa7d1eb48eb72c339bfa3ab5ebec0.png)
Для того чтобы правильно отобразить все точки облака и векторы Стокса на CCD матрице телескопе, необходимо выполнить обратные повороты.
![](https://habrastorage.org/getpro/habr/upload_files/38c/dcb/108/38cdcb1080d5b12fafb2357bcfb6e0de.png)
![](https://habrastorage.org/getpro/habr/upload_files/180/8f5/cb3/1808f5cb3d67cfeacf41c3a3b2e05649.png)
![](https://habrastorage.org/getpro/habr/upload_files/009/efd/688/009efd6888a3921da9796a6d4fe9623c.png)
![](https://habrastorage.org/getpro/habr/upload_files/2e6/22b/443/2e622b4431ae7b2ce90afc93e9b336a8.png)
![](https://habrastorage.org/getpro/habr/upload_files/c14/2f3/549/c142f35492d3e51db57edb10f3abaef5.png)
![](https://habrastorage.org/getpro/habr/upload_files/23a/01a/4ff/23a01a4ffdb86f6c6853c2fd9fbb8b01.png)
С математикой, которую выводил самостоятельно усё.
Написал я программу , которая рассчитана под линукс, используется C++11 concurency, распараллеливание на ядра процессора (сколько ядер, столько одновременно пакетов фотонов рассеиваются), запускал на кластере суперкомпьютера, у которого было 24 ядра. Получил результаты, но нужно было их как-то проверить.
Начальник принял решение - написать другую программу, которая решит систему интегральных уравнений Фредгольма для простой модели со звездой, сферическим облаком вокруг неё и Рэлеевским рассеиванием света. Для этой задачи я решил использовать питон, из-за удобства библиотек (да и интересно было что-то на питоне написать).
Далее была необходима нормализация параметров, т.к. разные программы давали результаты в разных величинах измерения. Результаты совпали, но для большей убедительности результаты сравнили также с программой известного учёного Cornelis Dullemond RADMC3D. Там результаты тоже совпали.
Были конференции, прикладываю презентацию одной.
А также видео другой (выступает мой бывший начальник) :
И также написана научная публикация (это уже после окончания моей работы, её к сожалению выложить не могу).
Ну вот - как-то так, поработал, написал свой первый софт, подтянул математику, мне понравилось :)
vkams
Навеяло трогательные воспоминания о начале моей работы в теоротделе ИВТ РАН: тоже писал на С программу, которая рассчитывала рассеяние методом Монте-Карло, и даже уравнения Фредгольма изучил, хотя не пригодились. Компьютеры тогда были - не чета нынешним: 386DX 25 МГц, а больше 286 16 МГц 1 Мб оперативки.