Моделирование фотонных кристаллов в программной среде MEEP. Знакомство со средой программирования

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Оглавление

Оглавление

1. Вступление

2.Теоретическая часть

2.1 Теория фотонных кристаллов

2.2 Уравнения Максвелла

2.3 Finite Difference Time Domain(FTDT)

2.4 PML-слой

3. Описание программных продуктов

3.1 Установка дополнений для гостевой операционной системы

3.2 Папки с общим доступом на локальном компьютере

3.3 Установка программы

3.4 Описание программных продуктов

3.5 Описание библиотек и утилит, используемых MEEP и MPB

4. Описание языка программирование и разработанных программ

4.1 Первая программа

4.2 Запуск программы на выполнение и представление данных

4.3 Программа, описывающая циклы и условные операторы

4.4 Пример программы , анализирующей структуру кристалла

4.5 Расчет двумерных структур размером 3х3

5. Перспективы и применение фотонных кристаллов

Выводы

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

1. Вступление

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

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

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

Рис1. Примеры фотонных кристалов, встречающиеся в природе.

Фотонные кристаллы можно встретить и в природе. Крылья бабочек переливаются всеми цветами радуги. Драгоценный камень опал тоже является природным фотонным кристаллом.

Фотонные кристаллы обладают двумя важными свойствами:

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

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

Второе свойство является главной чертой, отличающее фотонные кристаллы от обычных дифракционных решеток.

Целью моего курсового проекта является изучение свойств фотонных кристаллов и их моделирование с помощью программных продуктов MEEP и MPB.

2.Теоретическая часть

2.1 Теория фотонных кристаллов

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

Одномерные (плёночные структуры типа «брэгговского зеркала»), в которых коэффициент преломления периодически изменяется в одном пространственном направлении. Такие фотонные кристаллы состоят из параллельных друг другу слоев различных материалов с разными коэффициентами преломления и могут проявлять свои специфические оптические свойства только в одном пространственном направлении, перпендикулярном слоям.

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

Трёхмерные (структуры в виде сверхрешёток), в которых коэффициент преломления периодически изменяется в трёх пространственных направлениях.

Как и электрические среды в зависимости от ширины запрещённых и разрешённых зон, фотонные кристаллы можно разделить на проводники (света), диэлектрики (зеркала), полупроводники (выборочное отражение) и сверхпроводники.

Схематично показано явление дифракции лучей света на периодических структурах различной размерности. При рассеянии фотонов на 1D- и 2D-структурах всегда находятся такие направления распространения дифрагировавших лучей, для которых условие максимума интерференции выполнено. Для одномерного кристалла - нити (2.1.1 а), такие направления образуют конические поверхности, а в двумерном случае (2.1.1 б) - совокупность отдельных, изолированных друг от друга лучей.

а) б) в)

Рис2,1,1. Примеры 1-но, 2-х, и 3-х мерных структур кристаллов

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

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

показанный на рис. 2.1.3 точечной прямой линией; c - скорость света в вакууме, Q - волновой вектор и щ(Q)=щ(?Q). Введение модуляции диэлектрической проницаемости с периодом d приводит к появлению зоны Бриллюэна размером 1/d. На краях зоны Бриллюэна закон дисперсии фотонов отклоняется от (1), как показано на рис. 2.1.3

Рис2.1.2 Периодическая двумерная рис2.1.3 Закон дисперсии электромагнитных волн в структура фотонном кристалле (схема расширенных зон).

Для такой волны среда квазиоднородна, причем ее эффективная диэлектрическая проницаемость соответствует нулевому волновому вектору (Q>0). Как следствие, для предельно длинных электромагнитных волн в бесконечном фотонном кристалле справедлив закон дисперсии (1) с диэлектрической проницаемостью

(2)

где f обозначает относительную долю объема, занимаемого средой е1. В приближении почти

свободных фотонов, когда разность |е1 - е2| мала, выражение (1) можно экстраполировать в область длин волн, где Qd~1. Тогда частота фотонов

(3)

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

(4)

которое является условием брэгговского рассеяния фотонов на решетке с периодом d.

2.2 Уравнения Максвелла

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

Уравнения Максвелла:

 (5)

А также:

(6)

При расчетах предполагается, что материал абсолютный диэлектрик, а значит µ=1.

2.3 Finite Difference Time Domain(FTDT)

Finite Difference Time Domain - Метод конечных разностей во временной области.

Этот метод используется для решения численных задач в электродинамике при расчете уравнений Максвелла.

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

Такая постановка задачи позволяет делать расчеты в одно-, двух- и трехмерных пространствах. Для упрощения расчетов сетки полей H и E сдвигаются друг относительно друга на h/2.

Рис 2.3.1 Значение и направление полей E и H во времени

Такая схема построения расчетов называется схемой Yee. Она позволяет формировать пространственно временную сетку. В каждой ячейке сетки заданы Значения магнитных и электрических полей, значения диэлектрических и магнитной проницаемости, проводимости.

Рис 2.3.2 Пространственная сетка Yee.

Рассмотрим двумерный случай. При постановке задачи предполагается, что у нас задана фотонная кристаллическая структура с известным периодом решетки, формой и размером стержней. Также известна магнитная проницаемость среды е1 и е2.

Предполагается, что в системе нет свободных зарядов и нет токов. И все пространство разбито на сетку с шагом ?x, ?y и ?t по осям X, Y и времени.

Уравнение Максвелла в общем виде, расписанное для каждой из 3-х координат будет выглядеть следующим образом:

(7)

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

(8)

Набор синих уравнений (8) имеет аббревиатуру TMz(traveling-magnetic mode - вектор E находится в плоскости перпендикулярной оси Z), а набор красных уравнений (8) имеет название TEz(traveling-electric mode - вектор E находится в плоскости параллельной оси Z). Как видно из системы уравнений, красные и синие уравнения не связаны между собой, а значит могут быть рассчитаны независимо друг от друга.

Теперь, воспользовавшись схемой Yee, мы можем получить следующие расчетные формулы:



(9)

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

Начальное условие: Начальные значения для E и H известны во всех точках для момента времени t=0. Для расчета значений на полушаге используется линейная аппроксимация вида:

(10)

Приведенная выше система уравнений (8) и расчетные формулы к ним (9, 10), являются основой математической базы программных продуктов MEEP и MPB. Реализация MEEP и MPB скрывает эти расчеты от пользователя и весь порядок расчетов берет на себя. Программисту требуется только задать структуру, дефекты в ней, расположение и вид источника, временной интервал в котором будут проводиться измерения. То есть от программиста требуется только задать начальные условия, воспользовавшись программными средствами, и просмотреть конечный результат работы.

2.4 PML-слой

Рис 2.4.1 Расчетная область, окруження PML-слоем.

Так как расчеты фотонных кристаллов происходят на ограниченной области, то необходимо задавать специальные условия для E и H. От правильности их задания зависит результат расчета поля на каждом шаге на любом промежутке времени. Существует несколько способов их задания. Одним из них являются условия PML. Именно этим способом я пользовался при написании программ.

PML (Perfect Matched Layer) - идеально согласованный слой. Основным достоинством PML является наличие маленького коэффициента отражения, который практически не зависит от угла падения луча. Расчетная область окружается слоем заданной толщины. Все лучи, попавшие в него, со временем поглощаются. Идеальным считается задавать толщину PML--слоя порядка 10% от размера расчетной области.

При этом уравнения Максвелла заменяются следующими 12-ю уравнениями:

(11)

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

(12)

3. Описание программных продуктов

3.1 Установка дополнений для гостевой операционной системы

Первой с чего мы должны начать после установки системы на виртуальную машину - это с установки дополнений для гостевой операционной системы. Эти дополнения нужны для более удобной работы программиста на виртуальной машине(ВМ). Особенно важные операции это использование общего буфера как для хоста, так и для гостевой операционной системы(ОС). Еще одной важной функцией дополнений гостевой ОС является: более оптимальные настройки ВМ под ресурсы компьютера. С ними гостевая ОС будет работать быстрее, что очень важно при моделировании прохождения световых волн через фотонные кристаллы. Сам процесс моделирования занимает много времени(порядка 30 секунд при моделировании простейших структур кристаллов и нескольких минут при 3D-моделировании). Размер, создаваемых временных файлов в процессе компиляции может достигать нескольких сотен мегабайт.

Опишем процесс подготовки виртуальной машины. В первую очередь монтируем файл VBoxGuestAdditions.iso через меню “Установка””Добавить виртуальный образ”. Затем, находясь в смонтированной папке вбиваем команду sudo sh ./VBoxLinuxAdditions-x86.run. И после установки всех драйверов - перезагружаем гостевую операционную систему. Если перед установкой была обнаружена старая версия “Дополнений гостевой операционной системы, то вначале будет выполнено удаление старой версии. Система перезагрузится, но без установки новой версии. После того как будет удалена старая версия процесс установки необходимо повторить.

В случае установки на мою ОС возникли проблемы. Программа установки не смогла определить версию X11. По указанным путям в скрипте необходимых папок не оказалось, а структура самих папок, которые были найдены в других директориях, отличались от предложенных. Эту проблему удалось решить путем обновления самого Virtual Box. При обновлении был также закачен новый образ гостевой ОС. После полного обновления VirtualBox и дополнений гостевой операционной системы, все настройки были спешно закончены.

3.2 Папки с общим доступом на локальном компьютере

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

Из гостевой ОС наш хост будет рассматриваться как удаленный компьютер. Обращение к хосту из гостевой ОС осуществляется по протоколу smb. Обращение к хосту идет не по сетевому имени, а по ip адресу. IP адрес хоста можно узнать с помощью команды ipconfig. Введя в адресную строку:

smb://<ip-адрес хоста>/ + <Enter>

можно увидеть все доступные для общего доступа папки

Настройка VirtualBox и дополнений к ней закончена. Теперь можно приступить к установке программ.

3.3. Установка программы

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

Второй способ гораздо проще и не требует абсолютно ни каких затрат сил и времени. Так уж случилось, что программы Meep и MPB уже входят в пакет обновлений для Linux версии Ubuntu 10.0 и выше.

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

Воспользовавшись вторым способом, я подгрузил программы Meep и MPB в виде обновлений прямо с официального сайта Ubuntu. При этом я воспользовался командой apt-get. Последующие компиляции примеров программ и разработка собственных, доказали полную работоспособность установленных программ Meep и MPB.

Стоит также отметить, что хоть программы и работают под Linux, но рекомендуемым дистрибутивом является Ubuntu. Доверяя рекомендациям разработчиков, я выполнял разработку программ для своей курсовой работы именно под этим дистрибутивом(Ubuntu v11). Установка Meep и MPB на другие версии Linux не производилась. На сайте разработчиков говорится, что остальные дистрибутивы поддерживают не все функции программ. Удостовериться в данном утверждении мне не получилось, так как это не входит в рамки курсовой и вообще не обязательно.

3.4 Описание программных продуктов

Основной целью написания этого и последующих курсовых проектов является моделирование в таких программных продуктах, как MEEP и MPB. Эти программные продукты предназначены для расчета наноструктур. MEEP и MPB разрабатывались исследовательской группой ученых из Массачусетского технологического университета(MIT), в которую входят Steven G. Johnson, John D. Joannopoulos и Marin Soljaиiж.

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

Б)Meep - программа для моделирования электромагнитных по, которая в своей основе использует “Метод конечных разностей во временной области”(FTDT). С помощью этого метода происходит расчет электромагнитных полей в различных структурах. При этом все пространство разбивается на дискретную сетку, а затем вычисления происходят через равные промежутки времени. Чем меньше временной и пространственный шаг, тем выше точность вычислений, а значит можно проводить моделирование практических задач с высокой точностью. Позволяет проводить расчеты в одно-, двух и трехмерных пространствах. Позволяет делать расчеты электрической и магнитной проницаемостей, эффекты дисперсии в электромагнитных полях, а также позволяет проводить расчеты нелинейных диэлектрических и магнитных материалов.

В процессе написания своего курсового проекта я учился программировать в среде MEEP, строя и изучая параллельно различные диэлектрические структуры. В связи с тем, что большая часть используемых структур будет использоваться и в MPB, изучение основ программирования на Scheme позволит, в последствии, перейти к моделированию реальных трехмерных фотонных структур. Изучение основ владения программными средствами проходило на моделирование двумерных периодических структур с различными видами дефектов. Проведем обзор языка программирования Scheme и моделей, которые удалось смоделировать.

Как отмечено в документации к MEEP значения е0, µ0 и с скрыты от программиста. Их значения уже занесены в расчетные формулы и задавать их значения не нужно. Значения этих констант уже занесено в некоторую одну единицу.

В связи с тем, что уравнение Максвелла инвариантно(увеличение всех величин в уравнении в 10 приводит к увеличению частоты в 10 раз), то это позволяет использовать безразмерные величины. Это может быть очень удобно, так как присваивая самостоятельно одной из величин размерность мы можем проводить расчеты в очень широком диапазоне частот и периодов решетки.

В MEEP скорость света принято приравнивать c=1, а - период решетки. Частота в MEEP определяется как f=c/a (или для щ=2рc/a). Эти условия позволяют определять f как f=a/л.

3.5 Описание библиотек и утилит, используемых MEEP и MPB

Meep и MPB и библиотеки к ним написаны на С и С++, но также вызываются библиотеки BLAS и LAPACK, которые обычно компилируются на фортране. Поэтому при компиляции могут возникнуть проблемы связанные с различием версий между компилятором и библиотеками. По умолчанию используются компилятор фортрана f77 и g77

BLAS и LAPACK библиотек необходимы для того, чтобы установить Harminv. Harminv не требуется для MEEP, но настоятельно рекомендуется при расчетах резонанса. Но они могут пригодиться при работе с MPB.

1)BLAS - расшифровывается как "Basic Linear Algebra Subroutines"(Базовые подпрограммы линейной алгебры). Используется при перемножении матриц и работе с векторами . В ней находятся все базовые операции линейной алгебры, которые могут использоваться при написании более сложных приложений, а также при работе с LAPACK. Используемые алгоритмы отличаются высокой оптимальностью.

Состоит из 3-х уровней:

Выполнение скалярных, векторных и вектор-вектор операций;

Выполнение матрица-вектор операций;

Выполнение матрица-матрица операций.

2)LAPACK - Linear Algebra PACKage - стандартная библиотека процедур на основе BLAS, для более сложных операций линейной алгебры таких как: инверсия матриц, диагонализация матриц и т.д. Использует, по большей части, 3-й уровень BLAS.

3)GNU Guile - является расширенным скриптовым языком программирования, основанным на Scheme. Используется как, основанная Libctl библиотеках, оболочка для Meep(по умолчанию используется С++ интерфейс).

4)Libctl - бесплатная, основанная на Guile(язык программирования для создания приложений в Linux, позволяющий включать их в другие программы в качестве плагинов) библиотека, позволяющая создавать управляющие файлы для научного моделирования.

5)MPI - Если необходимо, то MEEP может работать на параллельных машинах с распределенной памятью, а для этого используется стандартный интерфейс передачи сообщений MPI. Существует две разновидности стандарта MPI: MPICH и LAM. MPI не рекомендуется устанавливать на однопроцессорные системы.

MPI, конкретно к используемым продуктам, реализован так, что несколько программ могут одновременно обращаться к одному control файлу. C++ интерфейс не поддерживает данную функцию. Чтобы MPI работало, необходимо также установить HDF5.

6)H5utils программы - отвечает за визуализацию научных данных в бесплатный и легко переносимый формат HDF5. H5utils предоставляет простые инструменты по визуализации данных в PNG формате(2d моделирование), plain text, Vis5d, VTK.

В h5topng входят следующие утилиты:

· H5totxt и h5fromtxt - переводит данные из HDF5 файла в/из ASCII текст

· H5topng - переводит 2-d расчеты из HDF5 файла в PNG, с возможность задания опций при выводе в цвете.

· H5tov5d - переводит данные из HDF5 файла в 3d+ проект для дальнейшей визуализации в программе Vis5d. Чтобы воспользоваться этой утилитой Vis5d уже должна быть установлена на компьютере

· H5tovtk - переводит HDF5 в VTK формат свободной программы Visualization ToolKit

· H5math - создание и комбинация данных из HDF5 файла с простыми математическими выражениями.

· H5read.oct - плагин для GNU Octave для чтения 2-d из HDF5 файлов

· H5fromh4 и H5fromh5 - конвертирует HDF(версия 4) в HDF5.

· Hdf5-to-cdf - программа для конвертирования HDF5 в CDF формат

7)Libctl - бесплатная, основанная на Guile(язык программирования для создания приложений в Linux, позволяющий включать их в другие программы в качестве плагинов) библиотека, позволяющая создавать управляющие файлы для научного моделирования.

8)NLopt - бесплатная open-source библиотека для нелинейной оптимизации.

· простой интерфейс, позволяющий вызывать различные алгоритмы оптимизации, меняя только 1 параметр.

· Поддерживает оптимизацию с большим количеством параметров (некоторые алгоритмы поддерживают размеры в миллион параметров и тысячи констант)

· Позволяет использовать локальные и глобальные алгоритмы оптимизации

· Имеет алгоритмы безусловной оптимизации, связано условной оптимизации

4. Описание языка программирование и разработанных программ

4.1 Первая программа

фотонный кристалл операционный программа

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

Пример 1(Листинг программы test2.ctl):

;1)

(set! geometry-lattice (make lattice (size 10 8 no-size)))

;2)

(set! geometry (list

(make block (center 0 0) (size 1.5 1.5 infinity) (material (make dielectric (epsilon 12))) )

(make block (center -1.5 0) (size 1 1 infinity) (material (make dielectric (epsilon 12))) )

(make block (center 1.5 0) (size 1 1 infinity) (material (make dielectric (epsilon 12))) )

(make block (center 0 -1.5) (size 1 1 infinity) (material (make dielectric (epsilon 12)) )

(make block (center -1.5 -1.5) (size 1 1 infinity) (material (make dielectric (epsilon 12))) )

(make block (center 1.5 -1.5) (size 1 1 infinity) (material (make dielectric (epsilon 12))) )

(make block (center 0 1.5) (size 1 1 infinity) (material (make dielectric (epsilon 12))) )

(make block (center -1.5 1.5) (size 1 1 infinity) (material (make dielectric (epsilon 12))) )

(make block (center 1.5 1.5) (size 1 1 infinity) (material (make dielectric (epsilon 12)))) )

)

;3)

(set! sources (list (make source (src (make continuous-src (wavelength 1.5) (width 20))) (component Ez) (center -4 0)(size 0 2.5))))

;4)

(set! pml-layers (list (make pml (thickness 1.0))));5)

(set! resolution 20)

;6)

(run-until 200 (at-beginning output-epsilon) (to-appended "ez" (at-every 0.6 output-efield-z)))

Начнем описание программы команда за командой

1) Прежде чем объявить структуру кристалла, расположение и частоту источника и др. параметры, необходимо задать размеры расчетной области:

(set! geometry-lattice (make lattice (size 10 8 no-size)))

· geometry-lattice - задает область в которой будут проводиться расчеты

· size 10 8 no-size - задаем размеры расчетной области. Общий формат команды size x y z, где x y z - размер в каждой из 3-х координат соответственно. С помощью этой команды мы задали двухмерную область размером 10х8. Эта команда позволяет задавать трехмерные области, но в нашем случае размер вдоль оси Z не указан. Т.к мы не задали свойства расчетной области, то по умолчанию,е=1.

2) Теперь можно разместить внутри какие-нибудь объекты. Следующая структура:

(set! geometry (list

(make block (center 0 0) (size 1.5 1.5 infinity) (material (make dielectric (epsilon 12))) )

(make block (center -1.5 0) (size 1 1 infinity) (material (make dielectric (epsilon 12))) )

---------------

)

описывает массив элементов. Распишем ее более подробно:

· (set! geometry … )- задает объекты, которые будут участвовать в симуляции.

· (list … ) - создает список элементов, которые будут участвовать в симуляции. Используется с целью создания массивов и других форм объединения структур между собой. Например в If для разделения результатов true и else. Используется практически во всех командах. В нашем случае list объединяет несколько блоков в единый массив блоков.

· make block - задает блок(квадрат, прямоугольник - 2-у мерный случай; прямоугольный параллелепипед - 3-х мерный). Кроме этого поддерживаются следующие геометрические фигуры:

· sphere radius [number] - сфера;

· cylinder radius [number] height [number] axis [vector3] - цилиндр

· cone radius2 [number] - конус

· ellipsoid - элипс

· Рассмотрим параметры make block более подробно:

· (center 0 0) - указываем где будет располагаться центр блока(в нашем случае в начале координат). По умолчанию координате Z присваивается значение 0.

· (size 1.5 1.5 infinity) - задает размеры блока. В нашем случае описывается стержень с основанием 1.5*1.5 и бесконечной длинной.

· (material (make dielectric (epsilon 12))) - задает свойства материала. В нашем случае: Задаем диэлектрик с е=12.

· C помощью описанной выше команды, создается структура из 9 стержней, отстоящих друг от друга на расстоянии а=1.5. Все стержни, кроме центрального, имеют основание размером 1*1. Центральный имеет размер 1.5*1.5. Каждый из стержней имеет диэлектрическую проницаемость е1=12. Диэлектрическая проницаемость окружающей среды е2=1.

3) Теперь опишем параметры источника:

(set! sources (list (make source (src (make continuous-src (wavelength 1.5) (width 20))) (component Ez) (center -4 0)(size 0 2.5))))

· (set! sources … )- команда задающая вид излучения, размеры и его частоту, а также положение в пространстве.

· (make source (src … ) - создать отдельный источник. Meep поддерживает 3 вида источников:

· continuous-src - синусоида

· gaussian-src - гаусиановый импульсный источник пропорциональный ехр (- iщt - (Т - Т0) 2 / 2W2).

· custom-src - заданная пользователем функция

В нашем случае, задается синусоидальная функция с длинной волны л=1.5. Поле может быть описано также с помощью частоты (frequency х) и периода (period x)Параметр (width 20) указывает на то, что функция будет введена в моделируемую область не сразу, а постепенно в течении 20-ти временных интервалов.

· (component Ez) - источник поля является компонентой Jz.

· (center -4 0)(size 0 2.5) - задает координаты центра и размер источника излучения. В нашем случае, источник вытянут вдоль оси Y на 2.5 и не имеет размера в X и Z направлениях

4) Опишем толщину PML слоя. Известно, что при расчетах волна может пройти сквозь слой PML отразиться от границы расчетов и вернуться обратно в расчетную область. Чтобы этого избежать рекомендуется задавать толщину слой PML порядка 10% от размера расчетной области. Также следует помнить, что PML находится внутри расчетной области.

(set! pml-layers (list (make pml (thickness 1.0)))) - толщина PML=1.

5) (set! resolution 20) - Следующая команда задает количество пикселей на единицу длинны. Получаем, что размер расчетной области равен 200х160 pixels.

6) Программа готова, теперь ее нужно запустить:

(run-until 200 (at-beginning output-epsilon) (to-appended "ez" (at-every 0.6 output-efield-z)))

· (run-until 200 - запуск расчетов(на 200 временных интервалов).

· (at-beginning output-epsilon)(to-appended "ez" (at-every 0.6 output-efield-z))) - С помощью этих команд указывается какие компоненты поля нам нужно вывести. В теле программы можно проводить любые вычисления, а результат выводить с помощью этих команд. В нашем случае выводится две компоненты поле - это е и Ez.

(at-every 0.6 - значение компоненты Ez выводится каждые 0.6 временных интервала. Это впоследствии пригодится для создания GIF анимации.

(to-appended "ez" - причем выходной файл будет иметь приставку ez в названии.

4.2 Запуск программы на выполнение и представление данных

1) meep first2 - запуск программы

После запуска программы будут созданы два файла: first2-eps-000000.00.h5 и first2-ez.h5. Имя файла формируется из имени исходного ctl файла и величины, которую необходимо вывести.

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

2) Существует несколько способов представления данных. Один из них:

h5topng -S3 first2-eps-000000.00.h5 - из имеющегося *.h5 файла создается черно белый PNG с изображением расчетной области



Рис 4.2.1. Расчетная область first2.ctl

3) Это изображение не показательно и не дает представления о каких-либо качественных характеристиках фотонных кристаллов. Более интересен результат в виде анимированного изображения.

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

$ h5ls first2-ez.h5

ez Dataset {200, 160, 333/Inf} - 333 изображения можно построить.

б) Пользуясь средствами h5utils, создадим изображения для создания GIF анимации.

$ h5topng -t 0:332 -R -Zc dkbluered -a yarg -A first2-eps-000000.00.h5 first2-ez.h5

-t 0:332 Создает 333 PNG файла, изображающих состояние поля каждые 0,6 единиц времени

-Zc dkbluered - переводит стандартную черно-белую схему в сине(отриц.)-бело(ноль)-красную(полож.)

-a yarg -A - обводит диэлектрик сероватым контуром

в) Для создания GIF анимации можно воспользоваться любым доступным конвертором например:

convert first2-ez.t*.png first2-ez.gif

4.3 Программа, описывающая циклы и условные операторы

Следующая программа более функциональна и использует больше возможностей языка Scheme. Эта программа создает кристалл mat c диэлектрическою проницаемостью е=12. Затем создает в нем массив дырок с радиусом r=0.25, заполненных воздухом. В кристалл также внесен дефект в виде креста.

Рис 4.3.1 Расчетная область test4_42.ctl

;Листинг test4_42

;1

(define lx 20) ;Задаем размер расчетной области по оси X и Y

(define ly 20)

(define r 0.25) ;радиус отверстий

(define eps 12) ;материал

(define-param dmin -5) ;Задаем параметры цикла

(define-param dmax 5)

(define-param ds 1) ; величина шага

(define-param dpml 1) ;величина области PML;

(define-param l_place (+ (/ (+ dmin dmax) 2) 1)) ;задаем координаты линий

;--------------------------------------------

;2

(define mat (make dielectric (epsilon eps))) ;создаем диэлектрик mat с e=12

(set! geometry-lattice (make lattice (size lx ly no-size))) ;Задаем размеры расчетной области

;3

(set! geometry (list(make block (center 0 0) (size infinity infinity infinity) (material mat))))

; Заполняем расчетную область диэлектриком mat

;4

;Цикл по созданию отверстий в диэлектрике(материал - воздух)

(do ((x dmin (+ x ds))) ((> x dmax))

(do ((y dmin (+ y ds))) ((> y dmax))

;5

(if (and (not (eqv? x l_place)) (not (eqv? y l_place)))

(set! geometry (append geometry(list

(make cylinder (center x y) (material air) (radius r) (height infinity))

)))

)

)

)

(set! pml-layers (list (make pml (thickness dpml))))

(set! resolution 10)

;6

(set! sources (list(make source(src (make continuous-src(wavelength 2.5)(width 20)))

(component Ez)(center -9 0)(size 0 3))))

(run-until 200(at-beginning output-epsilon)(to-appended "ez" (at-every 0.6 output-efield-z)))

1) Комментарии отмечаются знаком “;”. Все что стоит за ним до конца строки компилятором не воспринимается.

2) define-param - используется если нам нужно изменить, в случае необходимости, значение параметра не открывая файл ctl. Прежде чем запустить компиляцию программы можно изменить один из параметров:

meep dmax=8 test4_42.ctl.

3) Вместо того, чтобы каждый раз прописывать свойства материала при описании каждого участка структуры, можно задать характеристики материала в виде отдельной переменной:

(define eps 12)

(define mat (make dielectric (epsilon eps))) - зная синтаксис предыдущего примера, легко определить, что мы определили диэлектрик с е=12.

4) Теперь переменную mat можно использовать при описании структур:

(set! geometry (list(make block (center 0 0) (size infinity infinity infinity) (material mat))))

· make block - создаем block бесконечных размеров посреди р

· (size infinity infinity infinity) (material mat)- бесконечных размеров из материала mat.

· (center 0 0) - посреди расчетной области.

5) Следующим шагом опишем реализацию циклов, формирующих дырки в кристаллической структуре.

Циклы реализуют одну из схем do-loop языка Scheme.

(do ((x dmin (+ x ds))) ((> x dmax))

(do ((y dmin (+ y ds))) ((> y dmax))

<тело цикла>

)

)

· dmin, dmax, ds - константы, задающие минимальное и максимальное значение счетчика цикла, ds - шаг счетчика цикла

· x dmin - цикл выполняется, начиная с x=dmin(dmin=-5)

· (+ x ds) - с шагом ds(ds=1)

· (> x dmax) - пока x не станет больше dmax (dmax=5).

Аналогично и для координаты y. В итоге мы описали структуру из дырок размером 11*11, с шагом ds=1.

6) Следующий оператор описывает расположение дефекта в кристалле

(if (and (not (eqv? x l_place)) (not (eqv? y l_place)))

(set! geometry (append geometry(list

(make cylinder (center x y) (material air) (radius r) (height infinity))

)))

)

· Математические операции в Scheme описываются несколько непривычным способом: вначале идет знак операции, а затем только числа.

Например, (2+3-1) можно представить как (+ 2 3 -1). “-“ здесь используется не как математическая операция, а как знак числа. Аналогичным образом формируются и более сложнее выражения.

· (not (eqv? x l_place)) - eqv? Одна из нескольких операций сравнения. Eqv? Сравнивает два числа между собой и возвращает true если они равны друг другу. Not - меняет результат логической операции на противоположный.

Благодаря структурам (5) и (6) происходит построение дырок внутри кристалла. При вычислении новых значений X и Y, происходит их сравнение с l_place(в нашем случае l_place=1). Дырка в кристалле делается только тогда, когда X или Y не равны l_place. Тем самым формируется дефект в виде креста.

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

4.4 Пример программы , анализирующей структуру кристалла

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

;flux1

;1

(define fcen 0.4) ;цент частот в пакете волн

(define df 0.1) ;ширина частот

;если no-bend=true, то не рисуем кристаллическую структуру

(define-param no-bend? false)

(define lx 20) ;Задаем размер расченой области

(define ly 20)

(define r 0.25) ;радиус отверстий

(define eps 12) ;материал

(define-param dmin -5) ;Задаем параметры цикла

(define-param dmax 5)

(define-param ds 1) ; величина шага

(define-param dpml 1) ;величина области PML;

(define-param l_place (+ (/ (+ dmin dmax) 2) 1)) ;задаем координаты линий

(define mat (make dielectric (epsilon eps))) ;создаем диэлектрик mat с e=12

;----------------------------------------------------------------------------

(set! geometry-lattice (make lattice (size lx ly no-size))) ;Задаем размеры расчетной области

(set! geometry (list(make block (center 0 0) (size infinity infinity infinity) (material mat))))

;2

; Заполняем расчетную область диэлектриком mat

(if (not no-bend?)

;Цикл по созданию отверстий в диэлектрике(материал - воздух)

(do ((x dmin (+ x ds))) ((> x dmax))

(do ((y dmin (+ y ds))) ((> y dmax))

(if (and (not (eqv? x l_place)) (not (eqv? y l_place))

(set! geometry (append geometry(list

(make cylinder (center x y) (material air) (radius r) (height infinity))

)))))))

(set! pml-layers (list (make pml (thickness dpml))))

(set! resolution 10)

;3

(set! sources (list(make source(src (make gaussian-src (frequency fcen) (fwidth df)))

(component Ez)(center -9 0)(size 0 3))))

;4

(define-param nfreq 100) ; количество частот в диапазон[fcen-df/2; fcen+df/2]

(define trans ; количество прошедшего потока

(add-flux fcen df nfreq

(if no-bend?

(make flux-region

(center (- dmin 2) l_place) (size 0 (* 2 ds)))

(make flux-region

(center (+ dmax 2) l_place) (size 0 (* 2 ds))))))

;5

(run-sources+

(stop-when-fields-decayed 50 Ez

(if no-bend?

(vector3 (- dmin 2) l_place)

(vector3 (+ dmax 2) l_place))

1e-3))

;6

(display-fluxes trans)

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

Для этих целей используется параметр no-bend.

No-bend=true - располагаем в расчетной области только источник излучения и проводим первый расчет

No-bend=false - заново располагаем в расчетной области источник излучения и кристалл. Производим второй расчет.

Значение no-bend=false - задается в программе по умолчанию

2) (if (not no-bend?) ;<текст> ) - Здесь и далее идет проверка параметра No-bend. Как видно из структуры, запуск программы нужно произвести два раза. Первый раз с no-bend=true, а второй раз с no-bend=false.

3) (set! sources

(list

(make source(src

(make gaussian-src (frequency fcen) (fwidth df)))

(component Ez)

(center -9 0)(size 0 3))))

· (make gaussian-src - задаем гаусиановый импульсный источник

· (frequency fcen) - с основной частотой fcen

· (fwidth df) - и шириной импульса df

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

(define trans ; количество прошедшего потока

(add-flux fcen df nfreq

(if no-bend?

(make flux-region

(center (- dmin 2) l_place) (size 0 (* 2 ds)))

(make flux-region

(center (+ dmax 2) l_place) (size 0 (* 2 ds))))))

· (define trans - структура, в которую будет заноситься значение энергий волны для каждой из частот.

· (add-flux fcen df nfreq - задает в каком диапазоне частот нужно производить расчеты. Nfreq - описывает количество частот из этого диапазона.

· (make flux-region - задает область пространства в которой нужно произвести расчет.

5) Затем задается в течении какого промежутка времени нужно производить расчеты. Существует два способа. Первый, использованный в программе:

(run-sources+

(stop-when-fields-decayed 50 Ez

(if no-bend?

(vector3 (- dmin 2) l_place)

(vector3 (+ dmax 2) l_place))

1e-3))

· (stop-when-fields-decayed - расчет производится до тех пор, пока величина поля не станет меньше определенной величины. В общем виде:

(stop-when-fields-decayed dT component pt decay-by

· 50 - dT - через какие промежутки времени нужно проводить измерение величины…

· Ez - component

· (if no-bend?

(vector3 (- dmin 2) l_place)

(vector3 (+ dmax 2) l_place)) - pt -… чтобы ее значение в одной из расчетных точек(при наличии или отсутствии решетки) …

· 1e-3)) - decay-by - … было бы меньше заданной величины(в нашем случае Ez<0.001)

Во втором способе, можно указать сколько раз нужно произвести расчеты поля после выключения источника излучения. Пример:

(run-sources+ 500) - расчет прекратится через 500 временных интервалов после выключения источника.

6) Все расчеты произведены и можно выполнить вывод результатов.

· (display-fluxes trans)

Результат:

flux1:, 0.35, 5.58516844105163e-5

flux1:, 0.351010101010101, 8.45812001141722e-5

flux1:, 0.352020202020202, 1.27883318656805e-4

flux1:, 0.353030303030303, 1.90004998820637e-4

flux1:, 0.354040404040404, 2.75747528407135e-4

flux1:, 0.355050505050505, 3.92351676416799e-4

flux1:, 0.356060606060606, 5.52562540438976e-4

flux1:, 0.357070707070707, 7.76885293373528e-4

flux1:, 0.358080808080808, 0.00109244029356672

flux1:, 0.359090909090909, 0.00152809295946778

…

В первой колонке задается частота излучения, во второй - количество прошедшей энергии.

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

7) Произведем анализ полученных результатов:

а) Запуск программы: Запуск программы производится два раза. В первый раз без кристаллической структуры:

meep no-bend?=true flux1.ctl | tee bend0.out

Во второй раз - со структурой:

meep flux1.ctl | tee bend.out

В обоих случаях результаты перенаправляются в файл. Bend0.out - хранит значение энергий падающего излучения; bend.out - хранит значения энергий прошедшего излучения.

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

б) Далее переведем результаты программы в формат dat, который воспринимается такими приложениями как Excel и Matlab:

grep flux1: bend0.out >bend0.dat

grep flux1: bend.out > bend.dat

в) Представление результатов:

Для представления результатов я пользовался средствами Excell:

Частота

Величина прошедшего излучения бралась как отношение количества падающего излучения(2 колонка в bend0.dat) к количеству прошедшего сквозь структуру излучения(вторая колонка в bend.dat).

Как видно из диаграммы количество прошедшего через дефект излучения при частоте 0,36 очень мало (~30%) , когда при остальных частотах величина прошедшего излучения выше 60%.

4.5 Рассчет двумерных структур размером 3х3

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

Листинг test2_2.ctl

;Создадим 2-х мерную область размером 10х8)

(set! geometry-lattice (make lattice (size 10 8 no-size)))

;Расположим в этой области стержни бесконечной длины и имеющие размер в ;плоскости ;1х1. Зададим е1=12.

(set! geometry (list

(make block (center 0 0) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center -1.5 0) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center 1.5 0) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center 0 -1.5) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center -1.5 -1.5) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center 1.5 -1.5) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center 0 1.5) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center -1.5 1.5) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))

(make block (center 1.5 1.5) (size 1 1 infinity)

(material (make dielectric (epsilon 12))))))

;Зададим характеристики источника излучения: синусоида с л=1,5; центр источника ;расположен в ;координатах (-4 0); размер источника по осям ХY (0 2,5); .

(set! sources (list

(make source

(src (make continuous-src (wavelength 1.5) (width 20)))

(component Ez)

(center -4 0)(size 0 2.5))))

;Зададим толщину PML слоя( толщина =1)

(set! pml-layers (list (make pml (thickness 1.0))))

;Зададим число пикселей на единицу расстояний

(set! resolution 20)

;Запускаем выполнение программы на 200 временных интервалов. Измерения будем ;проводить каждые 0,6 от временного интервала: 200/0.6=333 - количество записанных ;временных интервалов.

(run-until 200

(at-beginning output-epsilon)

(to-appended "ez" (at-every 0.6 output-efield-z)))

Рис 4.5.1 Распределение волны в среде с показателями: a=1.5; D=1; е1=1; е2=12; л=1.5;

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

Рис 4.5.2 Структура таже, но размер центрального Dцентр=1.5;

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

Рассмотрим теперь структуры более большие по размеру структуры. С числом ячеек в решетке порядка 11х11.

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

Рис 4.5.3 Структура 11х11 с показателями: f=0.5; е1=1; е2=12; D=1; a=1.5;

Как видно из рисунков, волна попадает в структуру кристалла с большим трудом. Большая часть лучей пытается обогнуть созданную решетку. Те же лучи, которые попали в кристаллическую сетку лучше распространяются в диагональном направлении. Попадая в соседнюю по диагонали клетку, лучи также начинают расходится по диагонали. Только такое распространение лучей позволяет волне проникать внутрь кристалла. Это подтверждает предположение о том, что нам удалось описать структуру с запрещенной зоной при f=0.5.

Рис. 4.5.4 Структура таже что и на рис 3.2.3, но источник расположен внутри кристалла

В следующей программе я построил такую же структуру, но разметил источник волн прямо по центру кристалла и задал размер 0.5х0.5. Этот размер меньше, чем размер стержня в 2 раза. На представленных картинках видно, что свет не может выйти из кристалла вдоль направлений X и Y. В диагональном направление часть энергии все-таки уходит и возбуждает соседние с ними ячейки. На рисунках четко видно, как свет выходит через углы структуры, но принимается сильно интерферировать внутри решетки. Что приводит к резким скачкам максимумов и минимумов интенсивности, которые в итоге суммируются и глушат друг друга.

5. Перспективы и применение фотонных кристаллов

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

рис 4.1 Футуристический компьютер частями кристалла), можно формировать

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

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

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

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

Выводы

- В ходе курсового проекта удалось научиться пользоваться программным продуктом MEEP и изучить язык программирования Scheme; писать программы для Meep.

- Удалось провести моделирование двухмерных структур с разными характеристиками решетки, видами дефектов и длинами волн излучения, проходящего через решетку

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

- Удалось изучить основные структуры языков программирования

- Удалось провести теоретический обзор, изученной литературы

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

1.www.ab-initio.mit.edu

2. “Photonic Crystals. Molding The Flow of Light” Second Edition. Автор: John J. Joannopulus

3. А.В Закиров, В.Д. Левченко, “Эффективный алгоритм для трехмерного моделирования распространения электро магнитных волн в фотонных кристаллах” Москва 2008г.