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

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

Федеральное государственное автономное

образовательное учреждение

высшего профессионального образования

"Сибирский федеральный университет"

Институт инженерной физики и радиоэлектроники

Кафедра физики конденсированного состояния вещества

.65, физика конденсированного состояния вещества

Дипломная работа

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

Выпускник

А.Е. Ануфриев

Научный руководитель

д.т.н., проф. В.В. Тюрнев

Красноярск 2011

Содержание

Введение

. Расчёт однородной линии передачи

. Обзор литературы и вывод основных формул

.1 Коаксиальная линия

.2 Линия квадратного сечения с проводником в центре симметрии

. Постановка задачи расчёта однородной линии

. Проведение теоретических исследований однородной линии

.1 Линия квадратного сечения

.2 Прямоугольная линия с двумя симметричными проводниками

. Расчёт неоднородной двухпроводной линии передачи. Первая модель

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

. Постановка задачи для первой модели неоднородной линии

. Проведение теоретических исследований первой модели неоднородной линии передачи

.1 Вычисление коэффициентов разложения потенциала

.2 Проверка граничных условий

.3 Вычисление  и коэффициента связи

. Расчёт неоднородной двухпроводной линии передачи. Вторая модель

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

. Постановка задачи для второй модели неоднородной линии

. Проведение теоретических исследований второй модели неоднородной линии передачи

.1 Вычисление коэффициентов разложения потенциала

.2 Проверка граничных условий

.3 Вычисление  и коэффициента связи

Заключение

Список использованных источников

Приложение

Введение

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

1. Расчёт однородной линии передачи

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

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

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

Добротность резонатора или линии Q можно выразить через следующее соотношение [1]:

, (1)

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

. (2)

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

Таблица 1 − Параметры материалов, широко применяемых в СВЧ технике.

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

. (3)

В этом выражении h - толщина подложки, D − толщина скин-слоя. Данная формула (3) является точной для плоского волновода. Применять её для микрополосковой линии можно при условии, что её ширина гораздо больше h.

Для меди на частоте 1 ГГц D=2.09 мкм то лщина h имеет значения порядка 1 мм, при этих условиях Q_c»500. Заметим, что Q_c<<Q_d, и учитывая соотношение (1), Q_c вносит больший вклад в значение внутренней добротности Q.

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

Линия, рассмотренная в данной работе отличается от микрополосковой тем, что на её параметры значительно влияют как поперечные, так и продольные координаты проводников в поперечном сечении, а электромагнитная энергия более равномерно распределена по объёму.

Исследование проводилось на основе готового программного модуля, написанного в среде Compaq Visual Fortran 6.5, использующей язык Fortran-95, цель работы которого - получать набор Z_m - волновых сопротивлений для каждой моды электромагнитных колебаний. Входные параметры модуля - количество проводников N, их радиусы r_i, координаты центров x_i ,y_i, высота и ширина оболочки линии, а также относительная диэлектрическая проницаемость. Данная среда программирования выбрана по той причине, что имеет в своём составе большое число библиотек процедур и функций. В числе этих библиотек есть библиотека IMSL (Internal Mathematical and Statical Library), позволяющая работать с комплексными числами, облегчающая решение матричных и векторных задач.

В процессе работы в программу были внесены дополнения: возможность получения матрицы погонной ёмкости C в Ф/м , матрицы погонной индуктивности L в Гн/м, нормированных амплитуд токов на i-м проводнике от m-й волны I_im.

Визуальная оболочка программы была написана с использованием среды Delphi 7, так как язык Object Pascal удобен для объектно-ориентированного программирования. Оболочка упрощает работу с координатами, позволяет автоматизировать процесс исследования, просматривать графики в ходе эксперимента, сохранять результаты исследования в формате Microsoft Excel. Вычисление добротности также происходит в теле оболочки.

2. Обзор литературы и вывод основных формул

Метод, реализованный в алгоритме, основан на аналитическом решении телеграфных уравнений:

 (4)

Продифференцируем второе уравнение системы (4) по координате, первое по времени и подставим первое во второе:

 (5)

Уравнение (5) называется уравнением Гельмгольца для линии. Решение этого уравнения ищем в виде вектора токов гармонической Т-волны, бегущей вдоль линии (в данном случае координатная ось z направлена вдоль линии), где I_m - вектор амплитуд m-й нормальной волны:

.(6)

В результате решения получается формула:

.(7)

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

.(8)

Формула (7) при подстановке выражения (8) преобразуется к следующему виду:

.(9)

Так как напряжение меняется по тому же закону, что и ток

, (10)

то при подстановке (10) в первое телеграфное уравнение получаем формулу, связывающую амплитуды напряжений и токов:

 (11)

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

Формулу (11), используя уравнение (9), можно выразить через обратную матрицу ёмкости.

 (12)

 (13)

Для расчёта матрицы погонной ёмкости можно перейти к погонным зарядам Q_j:

 (14)

Дальнейшим шагом в данном методе является применение конформного отображения с помощью эллиптических функций, в результате которого прямоугольная область с круглыми контурами проводников преобразуются в полуплоскость с гладкими контурами над ней. Конформным называется такое отображение, которое переводит одну область в другую, причём если в первой области существовало решение уравнения Лапласа DF =0, то решение будет существовать и во второй области.

Если x, y - координатные оси в старой системе координат, то u, v - в новой

.

Конформное отображение w(z) приводится в [2]:

 (15)

где k и  - дополнительные модули, K(k) и  - связанные полные нормальные эллиптические интегралы Лежандра первого рода, dn(z,k) - эллиптическая функция Якоби.

Значение k находится из условия [2]

 (16)

Погонный заряд Q_j на j-м проводнике связан с потенциалом формулой

 (17)

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

В итоге решение сводится к нахождению потенциала j - го проводника

, (18)

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

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

, (19)

где  - линейно независимые векторы единичных зарядов,

. (20)

Затем обратная матрица погонной ёмкости подставляется в формулу (13) ,что даёт нам амплитуды напряжений U_jm на j - м проводнике от m -го проводника, по которому течёт единичный ток.

Согласно (12) и (13), формулы для C _, L выглядят таким образом:

, (21)

, (22),

где U- матрица, составленная из элементов U_jm.

Учитывая соотношения (11), (21) и связь между амплитудами напряжений и токов через волновое сопротивление получаем уравнение:

 или  (23),

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

Добротность линии Q будет равна добротности резонатора, сделанного из отрезка такой линии [1]:

, (24)

W - запасённая энергия, P - усреднённая по времени мощность потерь,

, (25)

. (26)

Согласно закону приращения индуктивности [1],

, (27)

формула (24) преобразуется в

, (28)

где ∆Z- приращение волнового сопротивления при переходе от идеального проводника к реальному. В данном методе это переход от исходной конфигурации проводников к такой, когда ширина и высота прямоугольной линии увеличиваются на ∆, а диаметры проводников уменьшаются на ∆, где ∆ − толщина скин-слоя [1],

 (29)

Если учесть малость D по сравнению с размерами линии и проводников, то для аналитического определения ∆Z многопроводной линии применимо следующее выражение [1]:

, (30)

где n_i - нормаль к поверхности i-го проводника. Подставляя (30) в (28) получаем выражение для Q:

 (31)

Формула (31) применяется если есть аналитическое выражение для Z. Если же Z рассчитывается численными или приближёнными методами, то используется (28).

2.1 Коаксиальная линия

Коаксиальная линия является наиболее простой однородной однопроводной линией в расчёте благодаря цилиндрической симметрии. Рассчитывается аналитически.

Погонная ёмкость[5]:

 (32)

Погонная индуктивность, согласно (12):

 (33)

Волновое сопротивление:

 (34)

График зависимости волнового сопротивления от с точностью до :

Добротность, согласно (34) и (31):

 (35)

Ниже приведён график зависимости  с точностью до множителя :

Очевидно, что существует некий оптимальный параметр ,при котором добротность максимальна. Найдём p:

Тогда  (пФ/м);  (мкГн/м);  (Ом); .

При =1,=0,5 см,  , ∆=2,09 мкм (медь, частота 1 ГГц) добротность .

.2 Линия квадратного сечения с проводником в центре симметрии

Метод расчёта такой линии взят из источника [6]. Согласно этому методу, волновое сопротивление равно

 (Ом) (36)

Формула (36) является приближённой и применима при .

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

Согласно (31) и (36) получаем выражение для добротности коаксиальной линии:

 (37)

По аналогии с коаксиальной линией, оптимальное отношение

;  (Ом); .

При =1,=1 см, , ∆=2,09 мкм (медь, частота 1 ГГц) добротность .

3. Постановка задачи расчёта однородной линии

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

4. Проведение теоретических исследований однородной линии

В данной главе приведены результаты работы программы, основанной на описанном методе, найдены оптимальные параметры исследуемых линий.

.1 Линия квадратного сечения

Новый метод расчёта позволяет более точно находить параметры линии. Ниже представлены результаты работы программы для такой же линии, как на рисунке 7, чтобы удостовериться в соответствии результатов нового метода с ранее известным методом. (W=H=b=1 см, D=2.09 мкм, e=1, радиус проводника меняется от 5 мм до 0,25 мм с шагом 0,01 мм).

Оптимальное отношение

, ()

при максимальной добротности Q=1413; Z=77.5 Ом.

Отличие от предыдущих результатов составляет не более 1%.

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

Двухпроводная линия является основой для изготовления резонаторов. В линии такого типа будет уже 2 типа нормальных колебаний (мод). 1-я мода − синфазные колебания (чётная, even), 2-я мода - колебания в противофазе тока и напряжения на разных проводниках (нечётная, odd).

Рассматриваемая линия составлена из двух квадратных линий, между которыми убрана проводящая перегородка. W=10 мм, H=20 мм, D=2.09 мкм. Радиусы обоих проводников изменяются от 5 мм до 0,25 мм с шагом 0,01 мм, e=1.

Как видно, волновые сопротивления для разных мод различаются. Ниже представлен график отношения .

Добротность для обоих типов колебаний:

Графики добротности для двухпроводной линии интересны тем, что имеется два экстремальных значения Q, а также имеет место точка пересечения этих графиков. Добротность связана с поглощением энергии в линии соотношением (25), а сигнал, передающийся в данной линии является суперпозицией двух типов нормальных колебаний. Разное поглощение для разных мод может исказить форму сигнала, поэтому оптимальная конфигурация линии - когда добротности мод равны:

Q_e= Q_o=1405 ; Z_e=109,7 Ом ; Z_o=100,4 Ом ; W/d=4,95.

Пик чётной волны:

Q_e=1478 ; Z_e=86,2 Ом; Z_o=77,1 Ом; W/d=3,36.

Пик нечётной волны:

Q_o=1500,5; Z_e=83,0 Ом; Z_o=74,0 Ом; W/d=3,18.

В рамках этой задачи можно рассмотреть влияние расстояния между проводниками на параметры линии. В том случае, когда Q_e= Q_o ,при тех же размерах поперечного сечения диаметр проводников d=1,01 мм. В ходе данного исследования два проводника перемещаются от противоположных стенок навстречу друг-другу до соприкосновения.

Ниже приведены графики зависимости Z(R) и Q(R), где R - расстояние между центрами проводников в миллиметрах.

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

линия передача многопроводный резонатор

5. Расчёт неоднородной двухпроводной линии передачи. Первая модель

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

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

6. Обзор литературы и вывод основных формул для первой модели неоднородной линии

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

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

 (38)

Продифференцируем второе уравнение системы (4) по координате, первое по времени и подставим первое во второе:

 (39)

Уравнение (39),которое уже было представлено (5), называется уравнением Гельмгольца для линии. Решение этого уравнения ищем в виде вектора токов гармонической Т-волны, бегущей вдоль линии (в данном случае координатная ось z направлена вдоль линии), где I_m - вектор амплитуд токов m-й нормальной волны, U_m - вектор амплитуд напряжений m-й нормальной волны:

, (40)

 (41)

В результате решения получается формула:

 (42)

Волновой вектор зависит в данной задаче только от частоты и дилектрической проницаемости для m-й нормальной волны:

 (43)

Формула (42) при подстановке выражения (43) преобразуется к следующему виду:

 (44)

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

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

, (45)

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

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

Рассмотрим случай индуктивной связи. Будем считать, что колебательные контуры имеют вырожденные резонансные частоты (как в фильтре, сделанном на основе двухпроводной линии передачи):

 (46)

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

 (47)

Подставляя выражение (47) в (45) получим:

 (48)

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

, (49)

 (50)

Случай комбинированной связи даёт нам следующее:

, (51)

 (52)

Следует заметить, что формула (52) при =0 либо =0 с точностью до знака преобразуется соответственно в формулы (50) и (48). Положим, что знаки  и  определены, тогда подстановка (50) и (48) в (52) даст следующее выражение для коэффициента комбинированной связи:

 (53)

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

7. Постановка задачи для первой модели неоднородной линии

Для решения уравнения (44) относительно e_m необходимо знать матрицу погонной ёмкости . Для определения  необходимо найти функцию двумерного потенциала F(x, y). Для этого необходимо получить общее решение уравнения Лапласа

 (54)

Общее решение задачи Дирихле (54) можно строго представить в виде

 (55)

Функция  - потенциал произвольно распределенных зарядов с погонной плотностью Q_i на поверхности i-го цилиндрического проводника при отсутствии щелей на экране. Этот потенциал обращается в нуль на поверхности диэлектрика и вне его.

Функция  - поправка к потенциалу , всех цилиндрических проводников на s-й щели в экране, при отсутствии остальных щелей.

Общее выражение для потенциала  получено в работе методом конформных отображений. Оно имеет вид

, (56)

Где  (57)

 (58)

- неопределенные вещественные коэффициенты, характеризующие мультипольные моменты распределения заряда Q_i, l - порядок мультиполя, - комплексное число, сопоставляемое координатам x и y, z_i - комплексное число, сопоставляемое центру i-го проводника.

Функция комплексного аргумента wz(z) в (58) определяется формулой

, (59)

где  - эллиптическая функция Якоби (дельта амплитуды), k и  - дополнительные модули эллиптической функции,  и  - связанные полные эллиптические интегралы Лежандра первого рода [8].

Модуль k в формуле (59) является корнем уравнения

, (60)

где W - ширина, а H - полувысота экранированной линии передачи.

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

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

Пусть поправка к потенциалу  на поверхности s-й копланарной линии описывается функциями

на wh, oh: , (61)

на hw, ow:  (62)

Искомые функции представим рядами

 (63)

где . Использование функций Чебышева второго рода  обеспечивает быструю сходимость рядов, так как каждый член ряда удовлетворяет условию Мейкснера на краях проводников копланарной линии [1].

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

1)      являться гармоническими (условие (54)),

2)      удовлетворять условию эквипотенциальности цилиндрических проводников,

3)    удовлетворять граничным условиям на щели

, где  ,.

Здесь и далее индексом int обозначаются функции потенциала, определённые внутри экранированного блока (), индексом ext , соответственно, функции, определённые вне экранированного блока.

Чтобы найти , разложим в ряд Фурье функции (61) и (62):

на wh, oh: , (64)

на hw, ow: . (65)

Коэффициенты  находятся с использованием выражений (63) (вычислено на основе справочных данных из [9]):

на wh, oh: , (66)

на hw, ow:  (67)

Согласно (64)-(67), и, учитывая условие гармоничности,  определяется следующим образом:

на wh: , (68)

на oh: , (69)

на hw: , (70)

на ow:  (71)

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

Для варианта а) (Рисунок 23) функция  выбирается аналогично (64)-(71) в виде разложения в ряд Фурье:

на wh: ,  (72)

на oh: , (73)

на hw: , (74)

на ow:  (75)

Для варианта б) (Рисунок 23) функции ,  необходимо представить иначе, чем (61) - (62), а именно:

на wh, oh: , (76)

на hw, ow:  (77)

Функции (76) и (77), содержащие бесконечные пределы, целесообразно представить в виде интеграла Фурье:

на wh, oh: , (78)

на hw, ow: . (79)

Функции  определяется следующим образом (на основе справочных данных из [10]):

на wh, oh:  (80)

на hw, ow: , (81)

Используя выражения (78)-(81), мы можем записать формулы для  по аналогии с (72)-(75):

на wh: , (82)

на oh: , (83)

на hw: , (84)

на ow:  (85)

Вариант в) (Рисунок 23) допускает ещё более высокую точность нахождения , но в данной работе не будет рассматриваться.

На данный момент мы знаем функцию потенциала (55), которая пригодится для вычисления матрицы погонной ёмкости С и матрицы погонной индуктивности L, необходимых для решения уравнения (44) относительно  и .

Элементы обратной матрицы погонной ёмкости могут быть вычислены следующим образом:

, (86)

где  - линейно независимые векторы единичных зарядов,

 (87)

Прямая матрица С может быть получена путём обращения . Известно, что матрица погонной ёмкости является функцией диэлектрической проницаемости среды С=С(e_r), в то время, как матрица погонной индуктивности является функцией магнитной проницаемости среды L=L(m_r). Так как мы рассматриваем немагнитное заполнение линии (m_r=1), то нет необходимости в прямом вычислении матрицы L . Как было показано в [1] , матрицы С и L в случае однородной воздушной линии передачи являются зависимыми:

 (88)

Таким образом, вычисление L сводится к вычислению С для воздушной однородной линии передачи. Подстановка (88) в (44) даст следующий результат:

 (89)

Решение уравнения (89) относительно  в случае двухпроводной линии передачи даёт нам два значения эффективных диэлектрических проницаемостей для двух типов нормальных колебаний:  - для нечётной моды,  - для чётной моды, обычно .

Для резонатора СВЧ, сделанного из отрезка данной линии, формулу (45) можно записать как

 (90)

Рассмотрим резонансное взаимодействие проводников отрезка двухпроводной линии передачи, взаимодействующих по всей длине l_c

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

, (91)

где  - электрическая длина резонатора. Для чётных и нечётных колебаний электрические длины равны соответственно

 (92)

Соотношения (91) и (92) в совокупности приводят к выражению, связывающему частоты нормальных колебаний и эффективные диэлектрические проницаемости:

 (93)

В результате подстановки (93) в (90) получаем окончательное выражение для коэффициента связи, которое и будет использоваться далее:

 (94)

8. Проведение теоретических исследований первой модели неоднородной линии передачи

На основе вышеприведённых выкладок и формул (55) - (94) была написана программа на языке программирования Fortran в среде Compaq Visual Fortran 6.6 (см. Приложение А), вычисляющая 2NL коэффициентов  и Ns·Nj коэффициентов , где Ns - число щелей, Nj - максимальное число коэффициентов j, а также количество точек внутри щели, для которых записывается система уравнений . Программно реализованы метод а) и б) (Рисунок 23). Числа L, Nj допустимо брать не более 10, а число коэффициентов k требуется взять как можно больше

.

.1 Вычисление коэффициентов разложения потенциала

Ниже приведён пример работы программы для однопроводной (r=1 мм) воздушной линии квадратного сечения 10Ч10 мм с одной щелью в экране S=2 мм (Рисунок 25). L=10; Nj=10, Nk=1000. Размерность коэффициентов разложения функции потенциала

,

где  - заряд центрального проводника.

Таблица 2 −Результаты вычислений для однопроводной линии методом а)

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

Далее приведены результаты работы программы для двухпроводной (r₁= r₂=1 мм) воздушной линии прямоугольного сечения 20Ч10 мм с одной щелью в экране S=2 мм (Рисунок 26). L=20; Nj=10, Nk=1000. Размерность коэффициентов -

,

где  - заряд m-го проводника, при отсутствии зарядов на других проводниках.

Таблица 3 −Результаты вычислений методом а) для двухпроводной линии при m=1

Таблица 4 −Результаты вычислений методом а) для двухпроводной линии при m=2

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

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

8.2 Проверка граничных условий

Ниже приведены графики зависимости потенциала на проводниках двухпроводной воздушной линии прямоугольного сечения 20Ч10 мм с одной щелью в экране S=1 мм (Рисунок 27) в единицах

в зависимости от угла j. Число L варьируется, Nj=5, Nk=1000.

Судя по данным графикам, точность вычисления потенциала заметно растёт с ростом L. При L>5 погрешность вычисления не превышает 1%.

На следующих графиках представлено поведение электростатического потенциала на стенке с частичным отсутствием металлизации. Число Nk варьируется, L=5 Nj=5, S=1 мм.

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

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

8.3 Вычисление  и коэффициента связи

Представим зависимости  и  от относительной диэлектрической проницаемости материала, заполняющего отрезок линии (r₁= r₂=1 мм ,W=20 мм, 2H=10 мм, S=2 мм (Рисунок 24), L=5 , Nj=5, Nk=100).

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

Получим зависимость  от ширины щели S, которая будет варьироваться в пределах от 1 до 20 мм при = 100 ,Nk=500,1000.

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

Можно получить зависимость коэффициента связи от радиуса проводников при S=2 мм. Радиусы равны между собой и варьируются в пределах от 0,1 мм до 5 мм, = 100 ,Nk=500.

Далее приведены зависимости коэффициента связи от расстояния между проводниками при S=2 мм. Расстояние варьируется в пределах от 2 мм до 18 мм, проводники перемещаются вдоль оси x, сохраняя симметрию линии и координаты y_1,2=5 мм, = 100 ,Nk=500.

Следующим шагом будет получение зависимости относительной величины расщепления резонансных частот от координаты щели, ширина которой будет постоянна S=2 мм, = 100 ,Nk=500. Координата варьируется от 1 мм до 19 мм. Данная величина коэффициентом связи не является, хотя формально считается по формуле (94), поскольку изменение координаты щели нарушает условие равенства резонансных частот для связанных резонаторов.

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

9. Расчёт неоднородной двухпроводной линии передачи. Вторая модель

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

Данная линия при условии e_r >1 является неоднородной линией передачи, у которой различаются скорости распространения чётных и нечётных волн. Практическое значение такой линии заключается в том, что на основе таких линий могут быть созданы узкополосные монолитные диэлектрические фильтры.

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

На этом этапе работы были использованы результаты работ [7], [2]. Основные формулы для второй модели неоднородной линии такие же, как и для первой модели, они кратко приведены ниже.

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

 (95)

Продифференцируем первое уравнение системы (95) по времени, второе по координате и подставим первое во второе;

 (96)

Уравнение (96) называется уравнением Гельмгольца для линии. Решение этого уравнения ищем в виде вектора токов гармонической Т-волны, бегущей вдоль линии (в данном случае координатная ось z направлена вдоль линии), где I_m - вектор амплитуд токов m-й нормальной волны, U_m - вектор амплитуд напряжений m-й нормальной волны:

, (97)

. (98)

В результате решения получается формула:

 (99)

Волновой вектор зависит в данной задаче только от частоты и дилектрической проницаемости для m-й нормальной волны:

 (100)

Формула (99) при подстановке выражения (100) преобразуется к следующему виду:

 (101)

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

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

, (102)

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

11. Постановка задачи для второй модели неоднородной линии

Для решения уравнения (101) относительно e_m необходимо знать матрицу погонной ёмкости . Для определения  необходимо найти функцию двумерного потенциала F(x, y). Для этого необходимо получить общее решение уравнения Лапласа

 (103)

Общее решение задачи Дирихле (103) можно строго представить в виде

 (104)

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

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

Общее выражение для потенциала  получено в работе [2] методом конформных отображений. Оно имеет вид

, (105)

Где  (106)

 (107)

- неопределенные вещественные коэффициенты, характеризующие мультипольные моменты распределения заряда Q_i, l - порядок мультиполя, - комплексное число, сопоставляемое координатам x и y, z_i - комплексное число, сопоставляемое центру i-го проводника.

Функция комплексного аргумента wz(z) в (107) определяется формулой

, (108)

где  - эллиптическая функция Якоби (дельта амплитуды), k и  - дополнительные модули эллиптической функции,  и  - связанные полные эллиптические интегралы Лежандра первого рода [8].

Модуль k в формуле (108) является корнем уравнения

 (109)

где W - ширина, а H - полувысота экранированной линии передачи.

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

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

, (110)

 (111)

Условие (110) означает непрерывность потенциала на границе раздела диэлектрического заполнения и диэлектрического отверстия, где - относительная диэлектрическая проницаемость материала диэлектрического отверстия. Следующее условие (111) означает непрерывность нормальной составляющей вектора электрической индукции на данной границе раздела.

Используя соотношения

, (112)

, (113)

можно представить условие (111) в виде

, (114)

где  - производная по нормали к границе раздела диэлектриков заполнения и отверстия.

Верхний индекс ext, который уже упоминался выше, означает функцию, определённую вне диэлектрического отверстия, тогда как верхний индекс int означает функцию, определённую внутри диэлектрического отверстия

Пусть поправка к потенциалу  на границе раздела диэлектрика и s-го диэлектрического отверстия описывается гармонической функцией

 (115)

, (116)

которая аналогична функции (106) для цилиндрических проводников и отличается от него отсутствием логарифмического слагаемого, что связано с отсутствием свободных зарядов в диэлектриках. Эта поправка определена исключительно вне s-го диэлектрического отверстия, поскольку имеет точку разрыва при , - неопределенные вещественные коэффициенты, l - порядок мультиполя, - комплексное число, сопоставляемое координатам x и y, z_s - комплексное число, сопоставляемое центру s-го диэлектрического отверстия.

Пусть поправка , определенная исключительно внутри s-го диэлектрического отверстия, описывается гармонической функцией

 (117)

которая непрерывна внутри данного диэлектрического отверстия.

Различие в выборе функций (115) и (117) связано с требованиями конечности и непрерывности электрического потенциала.

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

·        - условию непрерывности потенциала на границе диэлектриков (16),

·         - условию непрерывности нормальной составляющей вектора электрической индукции на границе диэлектриков (111) (114).

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

Элементы обратной матрицы погонной ёмкости могут быть вычислены следующим образом:

, (118)

где  - линейно независимые векторы единичных зарядов,

 (119)

Прямая матрица С может быть получена путём обращения . Известно, что матрица погонной ёмкости является функцией диэлектрических проницаемостей среды С=С(e_r, e_s), в то время, как матрица погонной индуктивности является функцией магнитной проницаемости среды L=L(m_r). Так как мы рассматриваем немагнитное заполнение линии (m_r=1), то нет необходимости в прямом вычислении матрицы L . Как было показано в [1], матрицы С и L в случае однородной воздушной линии передачи являются зависимыми:

 (120)

Таким образом, вычисление L сводится к вычислению С для воздушной однородной линии передачи. Подстановка (120) в (101) даст следующий результат:

 (121)

Решение уравнения (121) относительно  в случае двухпроводной линии передачи даёт нам два значения эффективных диэлектрических проницаемостей для двух типов нормальных колебаний:  - для нечётной моды,  - для чётной моды.

Для резонатора СВЧ, сделанного из отрезка данной линии, формулу (102) можно записать как

 (122)

Рассмотрим резонансное взаимодействие проводников отрезка двухпроводной линии передачи, взаимодействующих по всей длине l_c

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

, (123)

где  - электрическая длина резонатора. Для чётных и нечётных колебаний электрические длины равны соответственно

 (124)

Соотношения (123) и (124) в совокупности приводят к выражению, связывающему частоты нормальных колебаний и эффективные диэлектрические проницаемости:

 (125)

В результате подстановки (125) в (122) получаем окончательное выражение для коэффициента связи, которое и будет использоваться далее:

 (126)

12. Проведение теоретических исследований второй модели неоднородной линии передачи

На основе вышеприведённых выкладок и формул (104) - (126) была написана программа на языке программирования Fortran в среде Compaq Visual Fortran 6.6 (см Приложение), вычисляющая 2NL коэффициентов , 2NsL коэффициентов  и столько же коэффициентов  где Ns - число диэлектрических отверстий, L - максимальный порядок мультипольного разложения. Для достаточной точности вычислений L допустимо брать не более 10.

.1 Вычисление коэффициентов разложения потенциала

Ниже приведены результаты работы программы для двухпроводной (r₁= r₂=1 мм) симметричной линии прямоугольного сечения 20Ч10 мм с одним воздушным отверстием r_s=1 мм (Рисунок 39), , L=10.

Размерность коэффициентов кратна

,

где  - заряд m-го проводника, при отсутствии зарядов на других проводниках.

Таблица 5 − Результаты вычислений для двухпроводной линии при m=1

Таблица 6 − Результаты вычислений для двухпроводной линии при m=2

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

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

.2 Проверка граничных условий

Ниже приведены графики зависимости потенциала на проводниках двухпроводной линии прямоугольного сечения 20Ч10 мм с одним воздушным отверстием r_s=1 мм (Рисунок 40) в единицах  в зависимости от угла j. Число L варьируется, .

Судя по данным графикам, точность вычисления потенциала заметно растёт с ростом L. При L>5 погрешность вычисления не превышает 1%.

На следующих картах представлено распределение электростатического потенциала по площади поперечного сечения неоднородной двухпроводной линии при двух разных  и , при этом L=10, r_s=1 мм.

Рисунок 43 − Распределение потенциала по площади поперечного сечения при  в единицах . Заряженный проводник слева.

Рисунок 44 − Распределение потенциала по площади поперечного сечения при  в единицах . Заряженный проводник слева.

На рисунке 44 изображена зелёная линия, пересекающая воздушное отверстие, вдоль которой построен линейный график зависимости потенциала от координаты вдоль линии. Он представлен на рисунке 45.

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

.3 Вычисление  и коэффициента связи

Представим зависимость  от относительной диэлектрической проницаемости материала, заполняющего резонатор, изготовленный на основе отрезка линии (r₁= r₂= r_s=1 мм ,W=20 мм, 2H=10 мм, (Рисунок 3), L=10).

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

Получим зависимость  от радиуса воздушного отверстия, который будет варьироваться в пределах от 0,1 мм до 1 мм при = 100.

Следующим шагом будет получение зависимости модуля коэффициента связи от вертикальной координаты воздушного отверстия при постоянной горизонтальной координате =10 мм. Координата варьируется от 1 мм до 9 мм, , r₁= r₂= r_s=1 мм ,W=20 мм, 2H=10 мм, L=10.

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

Заключение

В ходе первого этапа работы была получена теоретическая и практическая подготовка для проведения расчётов и исследования свойств экранированных многопроводных линий, включающая в себя ознакомление с телеграфными уравнениями, с методом решения граничного уравнения Лапласа с помощью конформных отображений. Для выполнения численных расчётов в течении работы были изучены основы языка программирования Fortran 95, получены навыки использования dll-библиотек, написанных на языке Fortran к программам, написанным в среде Delphi 7.

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

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

Результатом второго этапа работы явилась программа написанная и отлаженная в среде Maple 5.4, а затем переведённая с языка среды Maple 5.4 на язык Fortran 95 в среде Compaq Visual Fortran 6.0, решающая систему уравнений, являющих собой граничные условия, относительно коэффициентов разложения потенциала. С помощью получившегося инструмента можно восстановить потенциал в любой точке внутри экранированного блока, вычислить такие параметры линии, как коэффициент связи, погонную ёмкость, погонную индуктивность, эффективные диэлектрические проницаемости для разных типов нормальных колебаний.

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

Результатом третьего этапа работы явилась программа на языке Fortran 95 в среде Compaq Visual Fortran 6.0, решающая систему уравнений, являющих собой граничные условия, относительно коэффициентов разложения потенциала. С помощью получившегося алгоритма можно восстановить потенциал в любой точке внутри экранированного блока, вычислить такие параметры линии, как коэффициент связи, погонную ёмкость, погонную индуктивность, эффективные диэлектрические проницаемости для разных типов нормальных колебаний.

Список использованных источников

1. Тюрнев В.В. Теория цепей СВЧ / В.В. Тюрнев. − Красноярск. ИПЦ КГТУ. 2003. − 199 с.

. Тюрнев В.В. Расчет поперечных волн в экранированной линии передачи, содержащей круговые цилиндрические проводники / В.В. Тюрнев // Радиотехника и электроника − 2006. − Т. 51, Вып. 7. − С. 839-842.

. Бахарев С.И. Справочник по расчету и конструированию СВЧ полосковых устройств / С.И. Бахарев, В.И. Вольман, Ю.Н. Либ и др. − М.: Радио и связь, 1982. − С. 328.

. Лифшиц Ю.А. Перспективы применения сред с повышенной диэл. проницаемостью в технике СВЧ ГИС / Ю.А. Лифшиц, Р.П. Сейсян // Микроэлектроника. − 1981. −Т. 10, вып. 5. − С. 433-439.

. Иродов И.Е. Электромагнетизм. Основные законы / И.Е. Иродов  − 4-е изд., испр. - М.: Лаборатория Базовых Знаний, 2002 − 320 с.: ил., С. 60.

. Фельдштейн А.П. Справочник по элементам волноводной техники / А.П. Фельдштейн − С. 94, таб. 3.1

. Тюрнев B.B. Анализ экранированной линии передачи, содержащей круговые цилиндрические проводники и копланарные линии на заземленном экране / B.B Тюрнев // Радиотехника и электроника. − 2007. − том 52, №11. − С. 1353­1359.

. Корн Г. Справочник по математике для научных работников и инженеров / Г. Корн, Т. Корн − М.: Наука, 1968.

. Бейтмен, Г. Высшие трансцендентные функции / Г. Бейтмен, А. Эрдейи − M.:Наука, 1974. − Т. 2. −С. 82-92.

10. Анго А. Математика для электро- и радиоинженеров./ А. Анго. − М.: Наука, 1965. − С. 478-485.

Приложение А

Листинг основной программы, вычисляющей волновое сопротивление однородной линии, на языке FORTRAN

SUBROUTINE anrod(er,hb,wb,l,n,error,outdata)       DFIMSLNONE, INTENT(IN)::N,L ! N - количество проводников/ L+1 высший порядок мультиполя

REAL,INTENT(IN)::Er,Hb,Wb ! Er - диэл. проницаемость/ 2Hb половина высоты линии/ Wb ширина линии

COMPLEX(8)Zi(N),Zm(N),Iim(N,N) ! Zi - координаты центров проводников на комплексной плоскости

! Zm - Волновые сопротивления мод / Iim - амплитуды токов от m-й моды на i-м проводнике(8) Ri(N)! Ri - радиусы проводников(8),INTENT(INout)::outdata(4*N+4*n*n) ! outdata - линейный массив для ввода и вывода данных из данной процедуры

REAL(8) Uim(N,N),Lim(N,N),Cim(N,N) ! Uim - напряжение от m-й моды на i-м проводнике / Lim и Cim - взаимные погонные индуктивности и ёмкости i -го и m -го проводников(8)::AM ! AM=m=k^2 параметр эллиптического интеграла(8)::K1 ! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода,PARAMETER::Pi=3.1415926536, Zc=376.73035,vellight=2.99792E+8         ! константы в единицах СИERROR,I,J,I1,L1,I2,L2,Nbig,M ! вспомогательные переменные

REAL:: Zr(8),ALLOCATABLE::Wi(:),Wil(:,:) ! Wi - отображённые координаты центров проводников на комплексной плоскости

! Wil - отображённые координаты точек на проводниках на комплексной плоскости

REAL(8),ALLOCATABLE:: A(:,:),B(:),X(:)(8)::dummyr(8) dummyc

!-----------------------------------------I=1,n ! извлечение входных данных(I)=(1,0)*outdata(I)+(0,1)*outdata(I+n)(i)=outdata(I+2*n)(Er<=0) THEN; ERROR=1; RETURN; ENDIF ! ERROR=1(Hb<=0) THEN; ERROR=2; RETURN; ENDIF       ! ERROR=2(Wb<=0) THEN; ERROR=3; RETURN; ENDIF ! ERROR=3(ANY(Ri<=0)) THEN; ERROR=4; RETURN; ENDIF ! ERROR=4( N<=0 ) THEN; ERROR=5; RETURN; ENDIF         ! ERROR=5( L<=0 ) THEN; ERROR=6; RETURN; ENDIF ! ERROR=6=Zc/(2*Pi*SQRT(Er)) !

AM=AMELL(K1) ! AM - параметр эллиптического интеграла m=k^2; K1 - связанный интеграл K'.

IF( ALLOCATED(Wi) .OR. ALLOCATED(Wil)) DEALLOCATE(Wi,Wil)(Wi(N),Wil(N,0:2*L),STAT=ERROR)( ERROR/=0 ) THEN; ERROR=7; RETURN; ENDIF                  ! ERROR=7I1=1,N ! Заготовки w-координат:

Wi(I1) =WMAP( Zi(I1) ) ! w-координаты отображений центра I1-окружностиL1=0,2*L(I1,L1)=WMAP( Zi(I1)+Ri(I1)*CEXP((0.,2.)*L1*Pi/(2*L+1)) ) ! w-координаты отображений 2L равноотстоящих точек на I1-окружности

ENDDO=N*2*L( ALLOCATED(A) .OR. ALLOCATED(B) .OR. ALLOCATED(X)) DEALLOCATE(A,B,X)(A(Nbig,Nbig),B(Nbig),X(Nbig),STAT=ERROR)( ERROR/=0 ) THEN; ERROR=8; RETURN; ENDIF ! ERROR=8

! ВСЕ МОДЫ ----------------------------------------------------------------------M=1,N ! M - все моды (т.е. заряды поочередно расположены по всем проводникам)

A=0; B=0; X=0I1=1,N ! I1 - проводник №1 (потенциал)L1=1,2*L ! L1 - точки на проводнике №1(IND1(I1,L1))=DLOG(CDABS( (Wil(I1,L1)-Wi(M))/(Wil(I1,L1)-CONJG(Wi(M)))*(Wil(I1,0)-CONJG(Wi(M)))/(Wil(I1,0)-Wi(M)) ))I2=1,N ! I2 - проводник №2 (заряды)L2=1,L ! L2 - порядок мультиполя на проводнике №2

dummyc= (Wil(I1,L1)-Wi(I2))**(-L2) - (DCONJG(Wil(I1,L1))-Wi(I2))**(-L2) &

& - (Wil(I1, 0)-Wi(I2))**(-L2) + (DCONJG(Wil(I1, 0))-Wi(I2))**(-L2)(IND1(I1,L1),IND2(I2,L2,1))= DREAL(dummyc)(IND1(I1,L1),IND2(I2,L2,2))=-DIMAG(dummyc)

ENDDO ! L2 - порядок мультиполя на проводнике №2! I2 - проводник №2 (заряды)! L1 - точки на проводнике №1! I1 - проводник №1   (потенциал)ERSET(3,0,0) ! Для ошибок уровня=3 в процедурах IMSL отключить (0) PRINT и STOPERSET(4,0,0) ! Для ошибок уровня=4 в процедурах IMSL отключить (0) PRINT и STOP

CALL DLSARG(Nbig,A,Nbig,B,1,X) ! IMSL-Solve a real general system of linear equations with iterative refinement: AX=B(IERCD()==1) THEN; ERROR= 9; RETURN; ENDIF ! ERROR= 9 Poor matrix A in DLSARG(IERCD()==2) THEN; ERROR=10; RETURN; ENDIF ! ERROR=10 Singular matrix A in DLSARG

CALL ERSET(3,2,2) ! Для ошибок уровня=3 в процедурах IMSL установить по умолчанию (2) PRINT и STOPERSET(4,2,2) ! Для ошибок уровня=4 в процедурах IMSL установить по умолчанию (2) PRINT и STOPI1=1,N ! I1 - проводник №1L1=1,2*L ! L1 - точки на проводнике №1= DLOG(CDABS( (Wil(I1,L1)-DCONJG(Wi(M))) / (Wil(I1,L1)-Wi(M)) ))I2=1,N ! I2 - проводник №2 (заряды)L2=1,L ! L2 - порядок мультиполя на проводнике №2

dummyc = (Wil(I1,L1)-Wi(I2))**(-L2) - (DCONJG(Wil(I1,L1))-Wi(I2))**(-L2)= dummyr + DREAL(dummyc)*X(IND2(I2,L2,1)) - DIMAG(dummyc)*X(IND2(I2,L2,2))

ENDDO ! L2 - порядок мультиполя на проводнике №2! I2 - проводник №2 (заряды)! L1 - точки на проводнике №1(I1,M)=dummyr ! Напряжение на I1-ом проводнике, когда единичный ток течёт по M-му проводнику! I1 - проводник №1! M - моды

Uim=Zr*Uim(A,B,X,STAT=ERROR)( ERROR/=0 ) THEN; ERROR=11; RETURN; ENDIF ! ERROR=11(Wi,Wil,STAT=ERROR)( ERROR/=0 ) THEN; ERROR=12; RETURN; ENDIF ! ERROR=12=(1/vellight*sqrt(Er))*UIM !вычисление погонной индуктивностиDLINRG(n,uim,n,cim,n) !вычисление погонной ёмкости (обратная матрица)=(sqrt(Er)/vellight)*cim !вычисление погонной ёмкостиDEVCRG(n,Uim,n,Zm,Iim,n) ! вычисление амплитуд токов и волновых сопротивлений (Zm-[U])*[I]=0j=1,n ! вывод полученных данных в массив

do I=1,n(I+n*(j-1)+3*n)=uim(i,j)(I+n*(j-1)+3*n+n*n)=Cim(i,j)(I+n*(j-1)+3*n+2*n*n)=Lim(i,j)(I+n*(j-1)+4*n+3*n*n)=real(Iim(i,j))j=1,n(3*n*n+3*n+j)=real(zm(j))

! конец___________________________________________      !вложенные процедурыFUNCTION IND1(I1,L1),INTENT(IN):: I1,L1 ! I1 - номер проводника; L1 - порядок мультиполя

IND1=(I1-1)*2*L+L1FUNCTION IND1FUNCTION IND2(I2,L2,M2),INTENT(IN):: I2,L2,M2 ! I2 - номер проводника; L2 - порядок мультиполя; M2=(1,2) - поляризация мультиполя (ReФ,ImФ)

IND2=(I2-1)*2*L+2*(L2-1)+M2FUNCTION IND2

COMPLEX(8) FUNCTION WMAP(Z) ! Отображение на верхнюю полуплоскость (w=u+iv), v>=0(8),INTENT(IN):: Z ! z=x+iy - комплексная координата в прямоугольнике (0<x<Wb); (0<y<2Hb)

COMPLEX(8) ZEJDNZEJDN ! IMSL: dn(z,m)=-ZEJDN(Z*K1/Hb,AM)/DSQRT(1-AM)FUNCTION(8) FUNCTION AMELL(K1)         ! Корень K(m)/K(1-m)=Wb/Hb, где m=AMELL, K(m) - полный эллиптический интеграл Лежандра I рода(8),INTENT(OUT):: K1! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода, m=k^2 - параметр интеграла

REAL(8) DELK,R,X1,X2,X,Y1,Y2,YDELK ! IMSL: K(m)=Wb/Hb=0; X2=1; X=0.5=DELK(X)/DELK(1-X)(Y==R) THEN(Y>R) THEN=X; Y2=Y(Y<R) THEN=X; Y1=Y(X1==0 .OR. X2==1) THEN=(X1+X2)/2=X1 + (X2-X1)/(Y2-Y1)*(R-Y1)(X<X1+(X2-X1)/10) X=X1+(X2-X1)/10(X>X2-(X2-X1)/10)         X=X2-(X2-X1)/10((X2-X1)/X <= 4*EPSILON(X)) EXIT=DELK(1-X)=XFUNCTION AMELLSUBROUTINE ANrod

Фрагмент кода на object pascal, вычисляющий добротность экранированной линии передачи

…// визуализация, объявление переменных, вычисление параметров эксперимента…j:=0 to numpoints do // цикл по всем точкам эксперимента

… // визуализация:= form6.SpinEdit2.Value; // L извлекается из текстового поля

form9.StaticText2.Caption:=inttostr(j);i:=0 to numrods-1 do begin // зануление массивов ввода-вывода[i,j]:=0;

newoutputd[i,j]:=0;;//поиск значения L, при котором не происходит ошибки начиная с максимально заданного

application.ProcessMessages;.StaticText4.caption:=inttostr(l);form9.tag=1 then break;i:=0 to numrods-1 do begin //цикл заполнения массивов ввода-вывода[i]:=rodsgeom[0,i,j]; // rodsgeom - хранилище данных о геометрии проводников[i+numrods]:=rodsgeom[1,i,j];[i+2*numrods]:=rodsgeom[2,i,j];

… //визуализация;

// первый вызов процедуры ANrod

unit5.anrod(shieldgeom[2,j],shieldgeom[1,j],shieldgeom[0,j],l,numrods,error1,outdata);

if error1=0 then // если не произошло ошибки, то сохраняем полученные волновые сопротивления в outputd

for i:=0 to numrods-1 do begin[i,j]:=(outdata[i+3*numrods+3*numrods*numrods]);;i:=0 to numrods-1 do begin // изменяем геометрию проводников[i]:= outdata[I+3*numrods+3*numrods*numrods]; //outdataprev - старые волновые сопротивления,outdata[...] -новые[i+2*numrods]:=rodsgeom[2,i,j]-0.5*shieldgeom[3,j]; // в shieldgeom[3,j] хранится толщина скин слоя

end;

// второй вызов процедуры ANrod с изменённой геометрией

unit5.anrod(shieldgeom[2,j],shieldgeom[1,j]+0.5*shieldgeom[3,j],shieldgeom[0,j]+shieldgeom[3,j],l,numrods,error2,outdata);

finally end;(l); // уменьшение L на случай, если произошла ошибка((error1<>9) and (error2<>9) ) or (l=0) ; // конец цикла поиска минимального L, не дающего ошибки

for i:=0 to numrods-1 do begin:=(outdata[i+3*numrods+3*numrods*numrods]-outdataprev[i]); // изменение волнового сопротивления для каждой модыdz<>0 then[i,j]:=outdataprev[i]/dz ; // в newoutputd помещаются добротности

…. // отрисовка графиков волнового сопротивления и добротности; // конец эксперимента

… // визуализация, отрисовка графиков, заполнение таблиц, и т.д

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

program ANROD41PORTLIBdfIMSL(4),parameter:: N=2,L=5,Ns=1,Nk=500,Nj=5,Nbig=N*2*L(8):: h2,Ri(N),Eps,Sb ,ERRABS,ERRREL,ERREST(1):: draws(8):: Zi(N),Wi(N),Wil(N,0:2*L),Zil(N,0:2*L),dummyc,Iim(n,n),Zm(n),zm2(n),Iim2(n,n),PARAMETER:: ow=101, hw=102,wh=103, oh=104 , metod=1, metodnk=1(8):: K1,AM,dx,BJCS(Ns,Nj,Nk),A(Nbig+Nj*Ns,Nbig+Nj*Ns),B(Nbig+Nj*Ns),X(Nbig+Nj*Ns),MX(Nbig+Nj*Ns,N),slots(Ns,4),aa,bb,dummyr,PsiInt,psiext,psislot,addd,BFF,CSSS,ksv(8):: Zslot(Ns,Nj),Wslot(Ns,Nj),zex(Ns,Nj),Wext(Ns,Nj),Zint(Ns,Nj),Wint(Ns,Nj),argl(4) k,i1,i2,i3,s,s1,s2,k3,ii,r1,mm,ARG1(8),PARAMETER::Pi=3.14159265358979323846264338327950288419716939937510d0, Zc=376.73035d0         ,vellight=2.99792E+8(8),PARAMETER::Wb=20,Hb=5,eps0=1,parameter::PHIn=100,WALLN=500(8),parameter::PHIREs=2*Pi/PHIn,WALLresW=Wb/WALLn,WALLresH=Wb/WALLn(8)::Zrad(N,0:PHIn),Wrad(N,0:PHIn),argl1,argl2,Zwall(Ns,0:WALLN),Wwall(Ns,0:WALLN)(8)::potrad(N,N,0:PHIn,0:1),psirad(N,N,0:PHIn,0:1),UU ,sumrad(N,N,0:PHIn,0:1),potwall(N,Ns,0:walln,0:1),psiwall(N,Ns,0:walln,0:1),sumwall(N,Ns,0:walln,0:1)(8)::SINNN1,SINNN2,EXP1,EXP2,EXP3,EXP4,EXP5,Aim(n,n),Aim1(n,n),Aim2(n,n),em(n),K0(8) ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10ggggg,f1,f2ARG1ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10(3, FILE = 'output.TXT')UMACH (-2, 3)PGOPT (-1, 255)ERSET (0, 0, 0)(3, *) CPSEC(),' SECONDS '

!WRITE (3, '(A50)') 'START***OK'=.false.=2*Hb(1)=(5,5)(2)=(15,5)=100(1)=1(2)=1=2=0=1e-10=0=1=0.001mm=20,180(1,1)=Sb(1,2)=wb/2(1,3)=0(1,4)=ow

! WRITE (3, '(A50)') 'LINE PARAMETERS *** OK'=AMELL(K1)I1 =1, N(I1)=WMAP( Zi(I1) )         L1=0,2*L(I1,L1)=Zi(I1)+Ri(I1)*exp((0,2)*L1*Pi/(2*L+1))(I1,L1)=WMAP(Zil(I1,L1))

!WRITE (3, '(A50)') 'CONFORM IMAGES OF RODS POINT *** OK's =1, Ns( slots(s,4)==wh ) thenK3 =1, Nj(s,K3)=slots(s,2)+(0,1)*(slots(s,3)-slots(s,1)/2+K3*slots(s,1)/(Nj+1))(s,K3)=WMAP(Zslot(s,K3))(s,K3)=Zslot(s,K3)-dx(s,K3)=WMAP(Zint(s,K3))(s,K3)=Zslot(s,K3)+dx(s,K3)=WMAP(zex(s,K3))( slots(s,4)==oh ) thenK3 =1, Nj(s,K3)=slots(s,2)+(0,1)*(slots(s,3)-slots(s,1)/2+K3*slots(s,1)/(Nj+1))(s,K3)=WMAP(Zslot(s,K3))(s,K3)=Zslot(s,K3)+dx(s,K3)=WMAP(Zint(s,K3))(s,K3)=Zslot(s,K3)-dx(s,K3)=WMAP(zex(s,K3))( slots(s,4)==ow ) thenK3 =1, Nj(s,K3)=slots(s,2)-slots(s,1)/2+K3*slots(s,1)/(Nj+1)+(0,1)*(slots(s,3))(s,K3)=WMAP(Zslot(s,K3))(s,K3)=Zslot(s,K3)+(0,1)*dx(s,K3)=WMAP(Zint(s,K3))(s,K3)=Zslot(s,K3)-(0,1)*dx(s,K3)=WMAP(zex(s,K3))( slots(s,4)==hw ) thenK3 =1, Nj(s,K3)=slots(s,2)-slots(s,1)/2+K3*slots(s,1)/(Nj+1)+(0,1)*(slots(s,3))(s,K3)=WMAP(Zslot(s,K3))(s,K3)=Zslot(s,K3)-(0,1)*dx(s,K3)=WMAP(Zint(s,K3))(s,K3)=Zslot(s,K3)+(0,1)*dx(s,K3)=WMAP(zex(s,K3))

! WRITE (3, '(A50)') 'SLOT POINTS (BORDER) *** OK's=1,Ns=getparamA(slots(s,4))=getparamB(slots(s,4))((slots(s,4)==oh) .or. (slots(s,4)==wh) ) thenk =1,Nki3 =1, Nj=DBESJN(i3,aa*k*slots(s,1)/2)=dcos(aa*k*slots(s,3)+Pi*i3/2)(s,i3,k)=-i3*2*BFF*CSSS/kk =1,Nki3 =1, Nj=DBESJN(i3,aa*k*slots(s,1)/2)=dcos(aa*k*slots(s,2)+Pi*i3/2)(s,i3,k)=-i3*2*BFF*CSSS/k

! WRITE (3, '(A50)') 'BESSEL*COSINE COEEFICIENTS *** DONE'

! WRITE (3, *) 'START FLOODING MATRICES A,B [AX=B] FOR RODS POINTS'getaim(eps)DEVCRG(n,Aim,n,Zm,Iim2,n)DLINRG(n,aim,n,aim1,n)getaim(eps0)=matmul(aim1,aim)DEVCRG(n,Aim2,n,Zm2,Iim,n)=real(zm2)*eps=(em(2)-em(1))/(em(2)+em(1))(3, *) dabs(ksv)!mm(draws==.true.) theni1=1,Nm=1,N(3, *)(3, *)(3, *)'M= ',M(3, *)'N=',i1

!CALL DWRRRN('PHI (angle)', PHIN+1, 2, potrad(m,i1,0:PhiN,0:1),phiN+1, 0)

!CALL DWRRRN('PSI (angle)', PHIN+1, 2, psirad(m,i1,0:PhiN,0:1),phiN+1, 0)DWRRRN('SUMMARY RODS POTENTIAL (angle)', PHIN+1, 2, sumrad(m,i1,0:PhiN,0:1),phiN+1, 0)

! read(3,*)s = 1 , Nsr1 = 0 , WALLn(s,r1)=WALLresW*r1(s,r1)=WMAP(Zwall(s,r1))m = 1 , N(1:Nbig+Nj*Ns)=MX(1:Nbig+Nj*Ns,M)r1 = 0 , WALLns1 = 1 , Ns=(Wwall(s1,r1)-Wi(m))=dlog(cdabs((Wwall(s1,r1)-dconjg(Wi(m)))/argl))I2 = 1 , NL2 = 1 , L=(Wwall(s1,r1)-Wi(I2))=dconjg(Wwall(s1,r1))-Wi(I2)=(argl1)**(-L2) - (argl2)**(-L2)=dummyr + dreal(dummyc)*X(IND2(I2,L2,1)) -dimag(dummyc)*X(IND2(I2,L2,2))         (m,s1,r1,1)=dummyr(m,s1,r1,0)=r1*wallREsW=0s = 1 , Ns=getparamA(slots(s,4))=getparamB(slots(s,4))(slots(s,4)==wh) thenii = 1 , Njk = Nk,1,-1=UU+X(IND3(ii,s))*dsin(k*aa*dimag(Zwall(s,r1)))*dsinh(k*aa*dreal(Zwall(s,r1)))/dsinh(k*aa*bb)*BJCS(s,ii,k)(slots(s,4)==oh) thenii = 1 , Njk = Nk,1,-1=UU+X(IND3(ii,s))*dsin(k*aa*dimag(Zwall(s,r1)))*dsinh(k*aa*(bb-dreal(Zwall(s,r1))))/dsinh(k*aa*bb)*BJCS(s,ii,k)(slots(s,4)==hw) thenii = 1 , Njk = Nk,1,-1=UU+X(IND3(ii,s))*dsin(k*aa*dreal(Zwall(s,r1)))*dsinh(k*aa*dimag(Zwall(s,r1)))/dsinh(k*aa*bb)*BJCS(s,ii,k)(slots(s,4)==ow) thenii = 1 , Nj((metod==1) ) thenk = Nk,1,-1=X(IND3(ii,s))*dsin(k*aa*dreal(Zwall(s,r1)))*dexp(-k*aa*dImag(Zwall(s,r1)))*((1-dexp(-2*k*aa*(bb-dImag(Zwall(s,r1)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,ii,k)=UU+addd((metodnk==2)) then=Nk+1=aa*dreal(Zwall(s,r1))=-aa*dImag(Zwall(s,r1))=-2*aa*(bb-dImag(Zwall(s,r1)))=-2*aa*bb=aa*slots(s,1)/2=-pI*Ii/2-pI/4dQdAGI(f2 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)=-ADDD*Ii*2.8284271247461900976033774484194D0/sqrt(PI)=uu+addd(metod==2) then=0=ii

!ARG3==dImag(Zwall(s,r1))=slots(s,1)/2=pI*Ii/2=slots(s,2)-dreal(Zwall(s,r1))DQDAGI(ggggg,K0,INTERVAL,ERRABS,ERRREL,addd,ERREST)=uu+X(IND3(ii,s))*addd*arg7/2(m,s,r1,1)=UU(m,s,r1,0)=r1*wallREsW(1:N,1:Ns,0:wallN,1)=potwall(1:N,1:Ns,0:wallN,1)+psiwall(1:N,1:Ns,0:wallN,1)(1:N,1:Ns,0:wallN,0)=potwall(1:N,1:Ns,0:wallN,0)s1=1,Nsm=1,N(3, *)(3, *)(3, *)'S= ',s1(3, *)'M=',M

!CALL DWRRRN('PHI (angle)', PHIN+1, 2, potrad(m,i1,0:PhiN,0:1),phiN+1, 0)

!CALL DWRRRN('PSI (angle)', PHIN+1, 2, psirad(m,i1,0:PhiN,0:1),phiN+1, 0)

!CALL DWRRRN('SUMMARY wall POTENTIAL (mm)', walln+1, 2, sumwall(m,s1,0:wallN,0:1),wallN+1, 0)

! CALL DWRRRN(' wall POTENTIAL (mm)', walln+1, 2, psiwall(m,s1,0:wallN,0:1),wallN+1, 0)(3, *) CPSEC(),' SECONDS '(3)getaim(er)(8) erM =1, NI1 =1, NL1 =1, 2*L(IND1(I1,L1))=dlog(cdabs((Wil(I1,L1)-Wi(M))/(Wil(I1,L1)-conjg(Wi(M)))*(Wil(I1,0)-conjg(Wi(M)))/(Wil(I1,0)-Wi(M))))

!WRITE (3, *) 'B(',IND1(I1,L1),')'I2 =1, N        L2 =1, L =(Wil(I1,L1)-Wi(I2))**(-L2)-(conjg(Wil(I1,L1))-Wi(I2))**(-L2)-(Wil(I1,0)-Wi(I2))**(-L2)+(conjg(Wil(I1,0))-Wi(I2))**(-L2)(IND1(I1,L1),IND2(I2,L2,1))= dReal(dummyc)

! WRITE (3, *) 'A(',IND1(I1,L1),IND2(I2,L2,1),')'(IND1(I1,L1),IND2(I2,L2,2))=-dImag(dummyc)

! WRITE (3, *) 'A(',IND1(I1,L1),IND2(I2,L2,2),')'      s =1, Ns=getparamA(slots(s,4))=getparamB(slots(s,4))( slots(s,4)==wh ) theni3 =1, Nj=0k =Nk,1,-1=(dsin(k*aa*dImag(Zil(I1,L1)))*dexp(k*aa*(dreal(Zil(I1,L1))-bb))*((1-dexp(-2*k*aa*(dreal(Zil(I1,L1)))))/(1-dexp(-2*k*aa*bb)))-dsin(k*aa*dImag(Zil(I1,0)))*dexp(k*aa*(dreal(Zil(I1,0))-bb))*((1-dexp(-2*k*aa*(dreal(Zil(I1,0)))))/(1-dexp(-2*k*aa*bb))))*BJCS(s,i3,k)=PsiInt+addd(IND1(I1,L1),IND3(i3,s))=PsiInt

! WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'( slots(s,4)==oh ) theni3 =1, Nj=0k =Nk,1,-1=(dsin(k*aa*dImag(Zil(I1,L1)))*dexp(-k*aa*dreal(Zil(I1,L1)))*((1-dexp(-2*k*aa*(bb-dreal(Zil(I1,L1)))))/(1-dexp(-2*k*aa*bb)))-dsin(k*aa*dImag(Zil(I1,0)))*dexp(-k*aa*dreal(Zil(I1,0)))*((1-dexp(-2*k*aa*(bb-dreal(Zil(I1,0)))))/(1-dexp(-2*k*aa*bb))))*BJCS(s,i3,k)=PsiInt+addd(IND1(I1,L1),IND3(i3,s))=PsiInt

! WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'( slots(s,4)==hw ) theni3 =1, Nj=0

! if ((metod==1) .or. (metod==2)) thenk =Nk,1,-1=(dsin(k*aa*dReal(Zil(I1,L1)))*dexp(k*aa*(dImag(Zil(I1,L1))-bb))*((1-dexp(-2*k*aa*(dImag(Zil(I1,L1)))))/(1-dexp(-2*k*aa*bb)))-dsin(k*aa*dReal(Zil(I1,0)))*dexp(k*aa*(dImag(Zil(I1,0))-bb))*((1-dexp(-2*k*aa*(dImag(Zil(I1,0)))))/(1-dexp(-2*k*aa*bb))))*BJCS(s,i3,k)=PsiInt+addd

! endif

! if (metod==2) then

!ENDIF(IND1(I1,L1),IND3(i3,s))=PsiInt

!WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'( slots(s,4)==ow ) theni3 =1, Nj=0

!if (metod==1) thenk =Nk,1,-1=dsin(k*aa*dReal(Zil(I1,L1)))=dexp(-k*aa*dImag(Zil(I1,L1)))=1-dexp(-2*k*aa*(bb-dImag(Zil(I1,L1))))=1-dexp(-2*k*aa*bb)=dsin(k*aa*dReal(Zil(I1,0)))=dexp(-k*aa*dImag(Zil(I1,0)))=1-dexp(-2*k*aa*(bb-dImag(Zil(I1,0))))

!EXP6=1-dexp(-2*k*aa*bb)=(SINNN1*EXP1*EXP2-SINNN2*EXP4*EXP5)*BJCS(s,i3,k)/EXP3=PsiInt+addd((metodnk==2)) then=Nk+1=aa*dReal(Zil(I1,L1))=-aa*dImag(Zil(I1,L1))=-2*aa*(bb-dImag(Zil(I1,L1)))=aa*dReal(Zil(I1,0))=-aa*dImag(Zil(I1,0))=-2*aa*(bb-dImag(Zil(I1,0)))=aa*slots(s,1)/2=-2*aa*bb=-pI*I3/2-pI/4dQdAGI(f1 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)=-ADDD*I3*2.8284271247461900976033774484194D0/sqrt(PI)=psiint+addd

! endif(IND1(I1,L1),IND3(i3,s))=PsiInt

! WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'            s1 =1, Nss2 =1, Nsi1 =1, Nji2 =1, Nj(IND3(i1,s1),IND3(i2,s2))=0

! WRITE (3, *) 'START FLOODING MATRICES A,B [AX=B] FOR SLOTS POINTS's =1, Ns=getparamA(slots(s,4))=getparamB(slots(s,4))i =1, Nj=Wint(s,i)-Wi(M)(IND3(i,s))=dlog(cdabs((Wint(s,i)-dconjg(Wi(M)))/argl))

! WRITE (3, *) dummyrI2 =1, NL2 =1, L=-((Wint(s,i)-Wi(I2))**(-L2)-(dconjg(Wint(s,i))-Wi(I2))**(-L2))

! WRITE (3, *) dummyc(IND3(i,s),IND2(I2,L2,1))= dreal(dummyc)

! WRITE (3, *) 'A(',IND3(i,s),IND2(I2,L2,1),')'(IND3(i,s),IND2(I2,L2,2))=-dimag(dummyc)

! WRITE (3, *) 'A(',IND3(i,s),IND2(I2,L2,2),')'  ( slots(s,4)==wh ) theni3 =1, Nj=0=0=0=dsin(i3*dacos( 2*(dImag(Zslot(s,i))-slots(s,3))/slots(s,1)) )k =Nk,1,-1=dsin(k*aa*dimag(Zint(s,i)))*dexp(k*aa*(dreal(Zint(s,i))-bb))*((1-dexp(-2*k*aa*(dreal(Zint(s,i)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,i3,k)=PsiInt+addd=dsin(k*aa*dimag(zex(s,i)))*dexp(k*aa*(bb-dreal(zex(s,i))))*BJCS(s,i3,k)=PsiExt+addd(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

! WRITE (3, *) 'A(',IND3(i,s),IND3(i3,s),')'( slots(s,4)==oh ) theni3 =1, Nj=0=0=0=dsin(i3*dacos( 2*(dimag(Zslot(s,i))-slots(s,3))/slots(s,1)) )k =Nk,1,-1=dsin(k*aa*dimag(Zint(s,i)))*dexp(-k*aa*dreal(Zint(s,i)))*((1-dexp(-2*k*aa*(bb-dreal(Zint(s,i)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,i3,k)=PsiInt+addd=dsin(k*aa*dimag(zex(s,i)))*dexp(k*aa*dreal(zex(s,i)))*BJCS(s,i3,k)=PsiExt+addd(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

! WRITE (3, *) 'A(',IND3(i,s),IND3(i3,s),')'( slots(s,4)==hw ) theni3 =1, Nj=0=0=0=dsin(i3*dacos( 2*(dreal(Zslot(s,i))-slots(s,2))/slots(s,1)) )k =Nk,1,-1=dsin(k*aa*dreal(Zint(s,i)))*dexp(k*aa*(dImag(Zint(s,i))-bb))*((1-dexp(-2*k*aa*(dImag(Zint(s,i)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,i3,k)=PsiInt+addd=dsin(k*aa*dreal(zex(s,i)))*dexp(k*aa*(bb-dimag(zex(s,i))))*BJCS(s,i3,k)=PsiExt+addd(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

! WRITE (3, *) 'A(',IND3(i,s),IND3(i3,s),')'( slots(s,4)==ow ) theni3 =1, Nj=0=0=0=dsin(i3*dacos( 2*(dreal(Zslot(s,i))-slots(s,2))/slots(s,1)) )k =Nk,1,-1=dsin(k*aa*dreal(Zint(s,i)))=dexp(-k*aa*dImag(Zint(s,i)))=1-dexp(-2*k*aa*(bb-dImag(Zint(s,i))))=1-dexp(-2*k*aa*bb)=SINNN1*EXP1*EXP2*BJCS(s,i3,k)/EXP3=PsiInt+addd((metodnk==2)) then=Nk+1=aa*dReal(Zint(s,i))=-aa*dImag(Zint(s,i))=-2*aa*(bb-dImag(Zint(s,i)))=0=0=0=aa*slots(s,1)/2=-2*aa*bb=-pI*I3/2-pI/4dQdAGI(f1 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)=-ADDD*I3*2.8284271247461900976033774484194D0/sqrt(PI)=psiint+addd((metod==1)) thenk =Nk,1,-1=dsin(k*aa*dreal(zex(s,i)))=dexp(k*aa*dimag(zex(s,i)))=SINNN2*EXP4*BJCS(s,i3,k)=PsiExt+addd

!WRITE (*, *) i3,psiext

!psiext=0((metodnk==2)) then=Nk+1=aa*dreal(zex(s,i))=-aa*dImag(zex(s,i))=-2*aa*(bb-dImag(zex(s,i)))=-2*aa*bb=aa*slots(s,1)/2=-pI*I3/2-pI/4dQdAGI(f2 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)=-ADDD*I3*2.8284271247461900976033774484194D0/sqrt(PI)=psiext+addd(metod==2) then=0=i3

!ARG3=dreal(zex(s,i))=dimag(zex(s,i))=slots(s,1)/2=pI*I3/2=slots(s,2)-dreal(zex(s,i))DQDAGI(ggggg,K0,INTERVAL,ERRABS,ERRREL,addd,ERREST)=addd*arg7/2(*, *) i3,psiext(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

! WRITE (3, *) 'A(',IND3(i,s),IND3(i3,s),')'

! WRITE (3, *) 'solving linear equation******'DLSARG(Nbig+Nj*Ns,A,Nbig+Nj*Ns,B,1,X)

! WRITE (3, *) 'Ok'WROPT (-6, 5,1)

! CALL DWRRRN('A', Nbig+Nj*Ns, Nbig+Nj*Ns, A,Nbig+Nj*Ns, 0)

! CALL DWRRRN('B', Nbig+Nj*Ns, 1, B,Nbig+Nj*Ns, 0)

!WRITE (3, *)'M= ',M

!CALL DWRRRN('COFFICIENTS (X)', Nbig+Nj*Ns, 1, X,Nbig+Nj*Ns, 0)

!READ (3,*)

!MA(1:Nbig+Nj*Ns,1:Nbig+Nj*Ns,M)=A(1:Nbig+Nj*Ns,1:Nbig+Nj*Ns)

!MB(1:Nbig+Nj*Ns,M)=B(1:Nbig+Nj*Ns)(1:Nbig+Nj*Ns,M)=X(1:Nbig+Nj*Ns)m=1,Nr1=0,PHIn(m,r1)=Zi(m)+Ri(m)*exp(2*(0,1)*r1*Pi/PHIn)(m,r1)=WMAP(Zrad(m,r1))m = 1 , N

! A(1:Nbig+Nj*Ns,1:Nbig+Nj*Ns)=MA(1:Nbig+Nj*Ns,1:Nbig+Nj*Ns,M)

! B(1:Nbig+Nj*Ns)=MB(1:Nbig+Nj*Ns,M)(1:Nbig+Nj*Ns)=MX(1:Nbig+Nj*Ns,M)i1 = 1 , Nr1 = 0 , PHIn=Wrad(i1,r1)-Wi(m)=dlog(CDABS((Wrad(i1,r1)-DCONJG(Wi(m)))/argl))I2 = 1 , NL2 = 1 , L=Wrad(i1,r1)-Wi(I2)=DCONJG(Wrad(i1,r1))-Wi(I2)=argl1**(-L2) - argl2**(-L2)=dummyr + dreal(dummyc)*X(IND2(I2,L2,1))-dimag(dummyc)*X(IND2(I2,L2,2))(m,i1,r1,1)=dummyr(m,i1,r1,0)=r1*PHIREs=0s = 1 , Ns=getparamA(slots(s,4))=getparamB(slots(s,4))(slots(s,4)==wh) thenii = 1 , Njk = Nk,1,-1=X(IND3(ii,s))*dsin(k*aa*dimag(Zrad(i1,r1)))*dexp(k*aa*(dreal(Zrad(i1,r1))-bb))*((1-dexp(-2*k*aa*(dreal(Zrad(i1,r1)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,ii,k)=UU+addd(slots(s,4)==oh) thenii = 1 , Njk = Nk,1,-1=X(IND3(ii,s))*dsin(k*aa*dimag(Zrad(i1,r1)))*dexp(-k*aa*dreal(Zrad(i1,r1)))*((1-dexp(-2*k*aa*(bb-dreal(Zrad(i1,r1)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,ii,k)=UU+addd(slots(s,4)==hw) thenii = 1 , Njk = Nk,1,-1=X(IND3(ii,s))*dsin(k*aa*dreal(Zrad(i1,r1)))*dexp(k*aa*(dImag(Zrad(i1,r1))-bb))*((1-dexp(-2*k*aa*(dImag(Zrad(i1,r1)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,ii,k)=UU+addd(slots(s,4)==ow) thenii = 1 , Njk = Nk,1,-1=X(IND3(ii,s))*dsin(k*aa*dreal(Zrad(i1,r1)))*dexp(-k*aa*dImag(Zrad(i1,r1)))*((1-dexp(-2*k*aa*(bb-dImag(Zrad(i1,r1)))))/(1-dexp(-2*k*aa*bb)))*BJCS(s,ii,k)=UU+addd((metodnk==2)) then=Nk+1=aa*dReal(Zrad(i1,r1))=-aa*dImag(Zrad(i1,r1))=-2*aa*(bb-dImag(Zrad(i1,r1)))=0=0=0=aa*slots(s,1)/2=-2*aa*bb=-pI*Ii/2-pI/4dQdAGI(f1 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)=-ADDD*Ii*2.8284271247461900976033774484194D0/sqrt(PI)=uu+addd(m,i1,r1,1)=UU(m,i1,r1,0)=r1*PHIREs(1:N,1:N,0:PHIN,1)=potrad(1:N,1:N,0:PHIN,1)+psirad(1:N,1:N,0:PHIN,1)(1:N,1:N,0:PHIN,0)=potrad(1:N,1:N,0:PHIN,0)=0r1=0,PHIn(1:N,1:N)=Aim(1:N,1:N)+sumrad(1:N,1:N,r1,1)/(1+Phin)subroutine getaimFUNCTION IND1(I1,L1),INTENT(IN):: I1,L1 ! I1 - номер проводника; L1 - порядок мультиполя

IND1=(I1-1)*2*L+L1FUNCTION IND1FUNCTION IND2(I2,L2,M2),INTENT(IN):: I2,L2,M2 ! I2 - номер проводника; L2 - порядок мультиполя; M2=(1,2) - поляризация мультиполя (ReФ,ImФ)

IND2=(I2-1)*2*L+2*(L2-1)+M2FUNCTION IND2FUNCTION IND3(j3,s3),INTENT(IN):: j3,s3      =2*L*N+Nj*(s3-1)+j3FUNCTION IND3(8) FUNCTION getparamA(s)(8),INTENT(IN)::s((s==wh) .or. (s==oh)) then=Pi/2/Hb=Pi/WbifFUNCTION getparamA(8) FUNCTION getparamB(s)(8),INTENT(IN)::s( (s==wh) .or. (s==oh) ) then=Wb=2*HbFUNCTION getparamB(8) FUNCTION WMAP(Z) ! Отображение на верхнюю полуплоскость (w=u+iv), v>=0(8),INTENT(IN):: Z ! z=x+iy - комплексная координата в прямоугольнике (0<x<Wb); (0<y<2Hb)

COMPLEX(8) ZEJDNZEJDN ! IMSL: dn(z,m)=-ZEJDN(Z*K1/Hb,AM)/DSQRT(1-AM)FUNCTION(8) FUNCTION AMELL(K1) ! Корень K(m)/K(1-m)=Wb/Hb, где m=AMELL, K(m) - полный эллиптический интеграл Лежандра I рода(8),INTENT(OUT):: K1 ! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода, m=k^2 - параметр интеграла

REAL(8) DELK,R,X1,X2,X,Y1,Y2,YDELK ! IMSL: K(m)=Wb/Hb=0; X2=1; X=0.5=DELK(X)/DELK(1-X)(Y==R) THEN(Y>R) THEN=X; Y2=Y(Y<R) THEN=X; Y1=Y(X1==0 .OR. X2==1) THEN=(X1+X2)/2=X1 + (X2-X1)/(Y2-Y1)*(R-Y1)(X<X1+(X2-X1)/10) X=X1+(X2-X1)/10(X>X2-(X2-X1)/10) X=X2-(X2-X1)/10((X2-X1)/X <= 4*EPSILON(X)) EXIT=DELK(1-X)=XFUNCTION AMELL(8) function ggggg(x)(8) X ,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10(4) a,b,ARG1ARG1ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10=DSIN(arg9+(arg8)*X)=arg1+1=DBESJN(arg1+1,X*ARG7)=arg1-1=DBESJN(arg1-1,X*ARG7)=dexp(arg4*x)= arg2*(arg3+arg5)*arg6

!WRITE (3, *) x,gggggfunction ggggg(8) function f1(x)(8) X,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10(4) ARG1ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10= (DSIN(X*ARG10)*DEXP(X*ARG2)*(1-DEXP(X*ARG3))-DSIN(X*ARG4)*DEXP(X*ARG5)*(1-DEXP(X*ARG6)))*DCOS(ARG7*X-ARG9)/((1-DEXP(ARG8*X))*X*DSQRT(X*ARG7))

! WRITE (3, *) f1function f1(8) function f2(x)(8) X,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9(4) ARG1ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10= DSIN(X*ARG10)*DEXP(X*ARG2)*(1-DEXP(X*ARG3))*DCOS(ARG7*X-ARG9)/((1-DEXP(ARG4*X))*X*DSQRT(X*ARG7))

! WRITE (3, *) f2function f2

Листинг программы на языке fortran, вычисляющей коэффициент связи резонатора на основе линии с отверстием в диэлектрике

program ANROD5PORTLIBdfIMSL(4),parameter:: N=2,L=10,Nbig=N*2*L,Ns=1(8):: h2,Ri(N),Eps,Eps2(Ns),Rs(Ns) ,ERRABS,ERRREL,ERREST(1):: draws(8):: Zi(N),Zs(Ns),Ws(Ns),Wi(N),Wil(N,0:2*L),Zil(N,0:2*L),dummyc,Iim(n,n),Zm(n),zm2(n),Iim2(n,n)(8):: K1,AM,dx,A(Nbig+L*Ns*4,Nbig+L*Ns*4),B(Nbig+L*Ns*4),X(Nbig+L*Ns*4),MX(Nbig+L*Ns*4,N),aa,bb,dummyr,PsiInt,psiext,psislot,addd,ksv(8)::dz,dw,Zsi(Ns,2*L), Wsi(Ns,2*L),Zse(Ns,2*L),Wse(Ns,2*L),Zsb(Ns,2*L),Wsb(Ns,2*L)(4) k,i1,i2,i3,s,s1,s2,k3,ii,r1,mm(8),PARAMETER::Pi=3.14159265358979323846264338327950288419716939937510d0, Zc=376.73035d0         ,vellight=2.99792E+8(8),PARAMETER::Wb=20,Hb=5,eps0=1,parameter::PHIn=100 ,WALLN=500(8),parameter::PHIREs=2*Pi/PHIn ,WALLresW=Wb/WALLn,WALLresH=2*hb/WALLn(8)::Zrad(N,0:PHIn),Wrad(N,0:PHIn),argl1,argl2,Zsss(Ns,0:phin),Wsss(Ns,0:phin)(8)::potrad(N,N,0:PHIn,0:1),psirad(N,N,0:PHIn,0:1),UU ,sumrad(N,N,0:PHIn,0:1),pots(N,Ns,0:phin,0:1),psis(N,Ns,0:phin,0:1),sums(N,Ns,0:phin,0:1),pot(N,Walln,walln,0:2)(8)::SINNN1,SINNN2,EXP1,EXP2,EXP3,EXP4,EXP5,Aim(n,n),Aim1(n,n),Aim2(n,n),em(n),K0,argl(8) ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10(3, FILE = 'output.TXT')UMACH (-2, 3)PGOPT (-1, 255)ERSET (0, 0, 0)(3, *) CPSEC(),' SECONDS '=2*Hb(1)=(5,5)(2)=(15,5)(1)=(10,6)=100(1)=1(1)=1(2)=1(1)=1=0.0001=AMELL(K1)I1 =1, N(I1)=WMAP( Zi(I1) )L1=0,2*L(I1,L1)=Zi(I1)+Ri(I1)*exp((0,2)*L1*Pi/(2*L+1))(I1,L1)=WMAP(Zil(I1,L1))I1 =1, Ns(I1)=WMAP( Zs(I1) )s =1, NsK3 =1, 2*L(s,K3)=zs(s)+Rs(s)*exp((0,2)*k3*Pi/(2*L+1))(s,K3)=WMAP(Zsb(s,K3))(s,K3)=zs(s)+(Rs(s)-dx)*exp((0,2)*k3*Pi/(2*L+1))(s,K3)=WMAP(Zsi(s,K3))(s,K3)=zs(s)+(Rs(s)+dx)*exp((0,2)*k3*Pi/(2*L+1))(s,K3)=WMAP(zse(s,K3))getaim(eps)DEVCRG(n,Aim,n,Zm,Iim2,n)DLINRG(n,aim,n,aim1,n)getaim(eps0)=matmul(aim1,aim)DEVCRG(n,Aim2,n,Zm2,Iim,n)=real(zm2)*eps=(em(2)-em(1))/(em(2)+em(1))(3, *) dabs(ksv)(3, *) CPSEC(),' SECONDS '(3)getaim(er)(8) erM =1, N     =0=0=0I1 =1, NL1 =1, 2*L(IND1(I1,L1))=dlog(cdabs((Wil(I1,L1)-Wi(M))/(Wil(I1,L1)-conjg(Wi(M)))*(Wil(I1,0)-conjg(Wi(M)))/(Wil(I1,0)-Wi(M))))I2 =1, N L2 =1, L=(Wil(I1,L1)-Wi(I2))**(-L2)-(conjg(Wil(I1,L1))-Wi(I2))**(-L2)-(Wil(I1,0)-Wi(I2))**(-L2)+(conjg(Wil(I1,0))-Wi(I2))**(-L2)(IND1(I1,L1),IND2(I2,L2,1))= dReal(dummyc)(IND1(I1,L1),IND2(I2,L2,2))=-dImag(dummyc)    I2 =1, NsL2 =1, L=(Wil(I1,L1)-Ws(I2))**(-L2)-(conjg(Wil(I1,L1))-Ws(I2))**(-L2)-(Wil(I1,0)-Ws(I2))**(-L2)+(conjg(Wil(I1,0))-Ws(I2))**(-L2)(IND1(I1,L1),nbig+IND2(I2,L2,1))= dReal(dummyc)(IND1(I1,L1),nbig+IND2(I2,L2,2))=-dImag(dummyc)(IND1(I1,L1),nbig+ns*L*2+IND2(I2,L2,1))= 0(IND1(I1,L1),nbig+ns*L*2+IND2(I2,L2,2))=0                    I1 =1, Ns    L1 =1, 2*L(nbig+IND1(I1,L1))=-(dlog(cdabs((Wsb(I1,L1)-Wi(M))/(Wsb(I1,L1)-conjg(Wi(M))) )) -dlog(cdabs((Wsi(I1,L1)-Wi(M))/(Wsi(I1,L1)-conjg(Wi(M))) ))-(er/eps2(i1))*dlog(cdabs((Wse(I1,L1)-Wi(M))/(Wse(I1,L1)-conjg(Wi(M))) ))+(er/eps2(i1))*dlog(cdabs((Wsb(I1,L1)-Wi(M))/(Wsb(I1,L1)-conjg(Wi(M))) )) )I2 =1, N L2 =1, L=((Wsb(I1,L1)-Wi(I2))**(-L2)-(conjg(Wsb(I1,L1))-Wi(I2))**(-L2))=dummyc-((Wsi(I1,L1)-Wi(I2))**(-L2)-(conjg(Wsi(I1,L1))-Wi(I2))**(-L2))=dummyc-(er/eps2(i1))*((Wse(I1,L1)-Wi(I2))**(-L2)-(conjg(Wse(I1,L1))-Wi(I2))**(-L2))=dummyc+(er/eps2(i1))*((Wsb(I1,L1)-Wi(I2))**(-L2)-(conjg(Wsb(I1,L1))-Wi(I2))**(-L2))(nbig+IND1(I1,L1),IND2(I2,L2,1))= dReal(dummyc)(nbig+IND1(I1,L1),IND2(I2,L2,2))=-dImag(dummyc)    I2 =1, NsL2 =1, L=((Wsb(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsb(I1,L1))-Ws(I2))**(-L2))=dummyc-(er/eps2(i1))*((Wse(I1,L1)-Ws(I2))**(-L2)-(conjg(Wse(I1,L1))-Ws(I2))**(-L2))=dummyc+(er/eps2(i1))*((Wsb(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsb(I1,L1))-Ws(I2))**(-L2))(i1/=i2) then=dummyc-((Wsi(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsi(I1,L1))-Ws(I2))**(-L2))(nbig+IND1(I1,L1),nbig+IND2(I2,L2,1))= dReal(dummyc)(nbig+IND1(I1,L1),nbig+IND2(I2,L2,2))=-dImag(dummyc)(i1==i2) then=-((Zsi(I1,L1)-Zs(I2)) )**(L2)(nbig+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,1))= dReal(dummyc)(nbig+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,2))=-dImag(dummyc)         I2 =1, NsL2 =1, L(i1==i2) then=((Wsb(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsb(I1,L1))-Ws(I2))**(-L2))(nbig+2*L*Ns+IND1(I1,L1),nbig+IND2(I2,L2,1))= dReal(dummyc)(nbig+2*L*Ns+IND1(I1,L1),nbig+IND2(I2,L2,2))=-dImag(dummyc)=-(zsb(I1,L1)-zs(I2))**(L2)         (nbig+2*L*Ns+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,1))= dReal(dummyc)(nbig+2*L*Ns+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,2))=-dImag(dummyc)                   DLSARG(Nbig+4*L*Ns,A,Nbig+4*L*Ns,B,1,X)(1:Nbig+4*L*Ns,M)=X(1:Nbig+4*L*Ns)m=1,Nr1=0,PHIn(m,r1)=Zi(m)+Ri(m)*exp(2*(0,1)*r1*Pi/PHIn)(m,r1)=WMAP(Zrad(m,r1))m = 1 , N(1:Nbig+4*L*Ns)=MX(1:Nbig+4*L*Ns,M)i1 = 1 , nr1 = 0 , phin=0=Zrad(i1,r1)=Wrad(i1,r1)=dw-Wi(m)=CDABS((dw-DCONJG(Wi(m)))/arg1)=dlog(argl)I2 = 1 , NL2 = 1 , L=dw-Wi(I2)=DCONJG(dw)-Wi(I2)=arg1**(-L2) - arg2**(-L2)=dummyr + dreal(dummyc)*X(IND2(I2,L2,1))-dimag(dummyc)*X(IND2(I2,L2,2))I2 = 1 , NsL2 = 1 , L=dw-Ws(I2)=DCONJG(dw)-Ws(I2)=arg1**(-L2) - arg2**(-L2)=dummyr + dreal(dummyc)*X(nbig+IND2(I2,L2,1))-dimag(dummyc)*X(nbig+IND2(I2,L2,2)) (m,i1,r1,1)=dummyr(m,i1,r1,0)=r1*PHIREs=0r1=0,PHIn(1:N,1:N)=Aim(1:N,1:N)+potrad(1:N,1:N,r1,1)/(1+Phin)subroutine getaimFUNCTION IND1(I1,L1),INTENT(IN):: I1,L1 ! I1 - номер проводника; L1 - порядок мультиполя

IND1=(I1-1)*2*L+L1FUNCTION IND1FUNCTION IND2(I2,L2,M2),INTENT(IN):: I2,L2,M2 ! I2 - номер проводника; L2 - порядок мультиполя; M2=(1,2) - поляризация мультиполя (ReФ,ImФ)

IND2=(I2-1)*2*L+2*(L2-1)+M2FUNCTION IND2

COMPLEX(8) FUNCTION WMAP(Z) ! Отображение на верхнюю полуплоскость (w=u+iv), v>=0(8),INTENT(IN):: Z ! z=x+iy - комплексная координата в прямоугольнике (0<x<Wb); (0<y<2Hb)

COMPLEX(8) ZEJDNZEJDN ! IMSL: dn(z,m)=-ZEJDN(Z*K1/Hb,AM)/DSQRT(1-AM)FUNCTION(8) FUNCTION AMELL(K1) ! Корень K(m)/K(1-m)=Wb/Hb, где m=AMELL, K(m) - полный эллиптический интеграл Лежандра I рода(8),INTENT(OUT):: K1 ! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода, m=k^2 - параметр интеграла

REAL(8) DELK,R,X1,X2,X,Y1,Y2,YDELK ! IMSL: K(m)=Wb/Hb=0; X2=1; X=0.5=DELK(X)/DELK(1-X)(Y==R) THEN(Y>R) THEN=X; Y2=Y(Y<R) THEN=X; Y1=Y(X1==0 .OR. X2==1) THEN=(X1+X2)/2=X1 + (X2-X1)/(Y2-Y1)*(R-Y1)(X<X1+(X2-X1)/10) X=X1+(X2-X1)/10(X>X2-(X2-X1)/10)         X=X2-(X2-X1)/10((X2-X1)/X <= 4*EPSILON(X)) EXIT=DELK(1-X)=XFUNCTION AMELL