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

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

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

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



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

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

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

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

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



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

12.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 с.

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

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

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

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

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

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

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

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

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

Приложение А

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

SUBROUTINE anrod(er,hb,wb,l,n,error,outdata)

USE DFIMSL

IMPLICIT NONE

INTEGER, 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-м проводнике

REAL(8) Ri(N)! Ri - радиусы проводников

REAL(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 -го проводников

REAL(8)::AM ! AM=m=k^2 параметр эллиптического интеграла

REAL(8)::K1 ! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода

REAL,PARAMETER::Pi=3.1415926536, Zc=376.73035,vellight=2.99792E+8 ! константы в единицах СИ

INTEGER ERROR,I,J,I1,L1,I2,L2,Nbig,M ! вспомогательные переменные

REAL:: Zr

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

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

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

REAL(8)::dummyr

COMPLEX(8) dummyc

!-----------------------------------------

do I=1,n ! извлечение входных данных

zi(I)=(1,0)*outdata(I)+(0,1)*outdata(I+n)

ri(i)=outdata(I+2*n)

enddo

IF(Er<=0) THEN; ERROR=1; RETURN; ENDIF ! ERROR=1

IF(Hb<=0) THEN; ERROR=2; RETURN; ENDIF ! ERROR=2

IF(Wb<=0) THEN; ERROR=3; RETURN; ENDIF ! ERROR=3

IF(ANY(Ri<=0)) THEN; ERROR=4; RETURN; ENDIF ! ERROR=4

IF( N<=0 ) THEN; ERROR=5; RETURN; ENDIF ! ERROR=5

IF( L<=0 ) THEN; ERROR=6; RETURN; ENDIF ! ERROR=6

Zr=Zc/(2*Pi*SQRT(Er)) !

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

IF( ALLOCATED(Wi) .OR. ALLOCATED(Wil)) DEALLOCATE(Wi,Wil)

ALLOCATE(Wi(N),Wil(N,0:2*L),STAT=ERROR)

IF( ERROR/=0 ) THEN; ERROR=7; RETURN; ENDIF ! ERROR=7

DO I1=1,N ! Заготовки w-координат:

Wi(I1) =WMAP( Zi(I1) ) ! w-координаты отображений центра I1-окружности

DO L1=0,2*L

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

ENDDO

ENDDO

Nbig=N*2*L

IF( ALLOCATED(A) .OR. ALLOCATED(B) .OR. ALLOCATED(X)) DEALLOCATE(A,B,X)

ALLOCATE(A(Nbig,Nbig),B(Nbig),X(Nbig),STAT=ERROR)

IF( ERROR/=0 ) THEN; ERROR=8; RETURN; ENDIF ! ERROR=8

! ВСЕ МОДЫ ----------------------------------------------------------------------

DO M=1,N ! M - все моды (т.е. заряды поочередно расположены по всем проводникам)

A=0; B=0; X=0

DO I1=1,N ! I1 - проводник №1 (потенциал)

DO L1=1,2*L ! L1 - точки на проводнике №1

B(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)) ))

DO I2=1,N ! I2 - проводник №2 (заряды)

DO 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)

A(IND1(I1,L1),IND2(I2,L2,1))= DREAL(dummyc)

A(IND1(I1,L1),IND2(I2,L2,2))=-DIMAG(dummyc)

ENDDO ! L2 - порядок мультиполя на проводнике №2

ENDDO ! I2 - проводник №2 (заряды)

ENDDO ! L1 - точки на проводнике №1

ENDDO ! I1 - проводник №1 (потенциал)

CALL ERSET(3,0,0) ! Для ошибок уровня=3 в процедурах IMSL отключить (0) PRINT и STOP

CALL ERSET(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

IF(IERCD()==1) THEN; ERROR= 9; RETURN; ENDIF ! ERROR= 9 Poor matrix A in DLSARG

IF(IERCD()==2) THEN; ERROR=10; RETURN; ENDIF ! ERROR=10 Singular matrix A in DLSARG

CALL ERSET(3,2,2) ! Для ошибок уровня=3 в процедурах IMSL установить по умолчанию (2) PRINT и STOP

CALL ERSET(4,2,2) ! Для ошибок уровня=4 в процедурах IMSL установить по умолчанию (2) PRINT и STOP

DO I1=1,N ! I1 - проводник №1

DO L1=1,2*L ! L1 - точки на проводнике №1

dummyr = DLOG(CDABS( (Wil(I1,L1)-DCONJG(Wi(M))) / (Wil(I1,L1)-Wi(M)) ))

DO I2=1,N ! I2 - проводник №2 (заряды)

DO L2=1,L ! L2 - порядок мультиполя на проводнике №2

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

dummyr = dummyr + DREAL(dummyc)*X(IND2(I2,L2,1)) - DIMAG(dummyc)*X(IND2(I2,L2,2))

ENDDO ! L2 - порядок мультиполя на проводнике №2

ENDDO ! I2 - проводник №2 (заряды)

ENDDO ! L1 - точки на проводнике №1

Uim(I1,M)=dummyr ! Напряжение на I1-ом проводнике, когда единичный ток течёт по M-му проводнику

ENDDO ! I1 - проводник №1

ENDDO ! M - моды

Uim=Zr*Uim

DEALLOCATE(A,B,X,STAT=ERROR)

IF( ERROR/=0 ) THEN; ERROR=11; RETURN; ENDIF ! ERROR=11

DEALLOCATE(Wi,Wil,STAT=ERROR)

IF( ERROR/=0 ) THEN; ERROR=12; RETURN; ENDIF ! ERROR=12

lim=(1/vellight*sqrt(Er))*UIM !вычисление погонной индуктивности

call DLINRG(n,uim,n,cim,n) !вычисление погонной ёмкости (обратная матрица)

Cim=(sqrt(Er)/vellight)*cim !вычисление погонной ёмкости

call DEVCRG(n,Uim,n,Zm,Iim,n) ! вычисление амплитуд токов и волновых сопротивлений (Zm-[U])*[I]=0

do j=1,n ! вывод полученных данных в массив

do I=1,n

outdata(I+n*(j-1)+3*n)=uim(i,j)

outdata(I+n*(j-1)+3*n+n*n)=Cim(i,j)

outdata(I+n*(j-1)+3*n+2*n*n)=Lim(i,j)

outdata(I+n*(j-1)+4*n+3*n*n)=real(Iim(i,j))

enddo

enddo

do j=1,n

outdata(3*n*n+3*n+j)=real(zm(j))

enddo

! конец___________________________________________

CONTAINS !вложенные процедуры

INTEGER FUNCTION IND1(I1,L1)

INTEGER,INTENT(IN):: I1,L1 ! I1 - номер проводника; L1 - порядок мультиполя

IND1=(I1-1)*2*L+L1

END FUNCTION IND1

INTEGER FUNCTION IND2(I2,L2,M2)

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

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

END FUNCTION IND2

COMPLEX(8) FUNCTION WMAP(Z) ! Отображение на верхнюю полуплоскость (w=u+iv), v>=0

COMPLEX(8),INTENT(IN):: Z ! z=x+iy - комплексная координата в прямоугольнике (0<x<Wb); (0<y<2Hb)

COMPLEX(8) ZEJDN

EXTERNAL ZEJDN ! IMSL: dn(z,m)

WMAP=-ZEJDN(Z*K1/Hb,AM)/DSQRT(1-AM)

END FUNCTION

REAL(8) FUNCTION AMELL(K1) ! Корень K(m)/K(1-m)=Wb/Hb, где m=AMELL, K(m) - полный эллиптический интеграл Лежандра I рода

REAL(8),INTENT(OUT):: K1! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода, m=k^2 - параметр интеграла

REAL(8) DELK,R,X1,X2,X,Y1,Y2,Y

EXTERNAL DELK ! IMSL: K(m)

R=Wb/Hb

X1=0; X2=1; X=0.5

DO

Y=DELK(X)/DELK(1-X)

IF(Y==R) THEN

EXIT

ELSEIF(Y>R) THEN

X2=X; Y2=Y

ELSEIF(Y<R) THEN

X1=X; Y1=Y

ENDIF

IF(X1==0 .OR. X2==1) THEN

X=(X1+X2)/2

ELSE

X=X1 + (X2-X1)/(Y2-Y1)*(R-Y1)

IF(X<X1+(X2-X1)/10) X=X1+(X2-X1)/10

IF(X>X2-(X2-X1)/10) X=X2-(X2-X1)/10

ENDIF

IF((X2-X1)/X <= 4*EPSILON(X)) EXIT

ENDDO

K1=DELK(1-X)

AMELL=X

END FUNCTION AMELL

END SUBROUTINE ANrod

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

…// визуализация, объявление переменных, вычисление параметров эксперимента…

for j:=0 to numpoints do // цикл по всем точкам эксперимента

begin

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

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

form9.StaticText2.Caption:=inttostr(j);

for i:=0 to numrods-1 do begin // зануление массивов ввода-вывода

outputd[i,j]:=0;

newoutputd[i,j]:=0;

end;

repeat //поиск значения L, при котором не происходит ошибки начиная с максимально заданного

application.ProcessMessages;

form9.StaticText4.caption:=inttostr(l);

if form9.tag=1 then break;

try

for i:=0 to numrods-1 do begin //цикл заполнения массивов ввода-вывода

outdata[i]:=rodsgeom[0,i,j]; // rodsgeom - хранилище данных о геометрии проводников

outdata[i+numrods]:=rodsgeom[1,i,j];

outdata[i+2*numrods]:=rodsgeom[2,i,j];

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

end;

// первый вызов процедуры 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

outputd[i,j]:=(outdata[i+3*numrods+3*numrods*numrods]);

end;

for i:=0 to numrods-1 do begin // изменяем геометрию проводников

outdataprev[i]:= outdata[I+3*numrods+3*numrods*numrods]; //outdataprev - старые волновые сопротивления,outdata[...] -новые

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;

dec(l); // уменьшение L на случай, если произошла ошибка

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

for i:=0 to numrods-1 do begin

dz:=(outdata[i+3*numrods+3*numrods*numrods]-outdataprev[i]); // изменение волнового сопротивления для каждой моды

if dz<>0 then

newoutputd[i,j]:=outdataprev[i]/dz ; // в newoutputd помещаются добротности

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

end; // конец эксперимента

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

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

program ANROD41

USE PORTLIB

USE dfIMSL

INTEGER(4),parameter:: N=2,L=5,Ns=1,Nk=500,Nj=5,Nbig=N*2*L

REAL(8):: h2,Ri(N),Eps,Sb ,ERRABS,ERRREL,ERREST

logical(1):: draws

COMPLEX(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)

INTEGER,PARAMETER:: ow=101, hw=102,wh=103, oh=104 , metod=1, metodnk=1

REAL(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

COMPLEX(8):: Zslot(Ns,Nj),Wslot(Ns,Nj),zex(Ns,Nj),Wext(Ns,Nj),Zint(Ns,Nj),Wint(Ns,Nj),argl

integer(4) k,i1,i2,i3,s,s1,s2,k3,ii,r1,mm,ARG1

REAL(8),PARAMETER::Pi=3.14159265358979323846264338327950288419716939937510d0, Zc=376.73035d0 ,vellight=2.99792E+8

REAL(8),PARAMETER::Wb=20,Hb=5,eps0=1

integer,parameter::PHIn=100,WALLN=500

COMPLEX(8),parameter::PHIREs=2*Pi/PHIn,WALLresW=Wb/WALLn,WALLresH=Wb/WALLn

COMPLEX(8)::Zrad(N,0:PHIn),Wrad(N,0:PHIn),argl1,argl2,Zwall(Ns,0:WALLN),Wwall(Ns,0:WALLN)

REAL(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)

REAL(8)::SINNN1,SINNN2,EXP1,EXP2,EXP3,EXP4,EXP5,Aim(n,n),Aim1(n,n),Aim2(n,n),em(n),K0

REAL(8) ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

external ggggg,f1,f2

COMMON ARG1

COMMON ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

OPEN (3, FILE = 'output.TXT')

CALL UMACH (-2, 3)

CALL PGOPT (-1, 255)

CALL ERSET (0, 0, 0)

WRITE (3, *) CPSEC(),' SECONDS '

!WRITE (3, '(A50)') 'START***OK'

draws=.false.

h2=2*Hb

Zi(1)=(5,5)

Zi(2)=(15,5)

Eps=100

Ri(1)=1

Ri(2)=1

Sb=2

ERRABS=0

ERRREL=1e-10

errest=0

INTERVAL=1

dx=0.001

do mm=20,180

slots(1,1)=Sb

slots(1,2)=wb/2

slots(1,3)=0

slots(1,4)=ow

! WRITE (3, '(A50)') 'LINE PARAMETERS *** OK'

AM=AMELL(K1)

do I1 =1, N

Wi(I1)=WMAP( Zi(I1) )

do L1=0,2*L

Zil(I1,L1)=Zi(I1)+Ri(I1)*exp((0,2)*L1*Pi/(2*L+1))

Wil(I1,L1)=WMAP(Zil(I1,L1))

enddo

enddo

!WRITE (3, '(A50)') 'CONFORM IMAGES OF RODS POINT *** OK'

do s =1, Ns

if ( slots(s,4)==wh ) then

do K3 =1, Nj

Zslot(s,K3)=slots(s,2)+(0,1)*(slots(s,3)-slots(s,1)/2+K3*slots(s,1)/(Nj+1))

Wslot(s,K3)=WMAP(Zslot(s,K3))

Zint(s,K3)=Zslot(s,K3)-dx

Wint(s,K3)=WMAP(Zint(s,K3))

zex(s,K3)=Zslot(s,K3)+dx

Wext(s,K3)=WMAP(zex(s,K3))

enddo

endif

if ( slots(s,4)==oh ) then

do K3 =1, Nj

Zslot(s,K3)=slots(s,2)+(0,1)*(slots(s,3)-slots(s,1)/2+K3*slots(s,1)/(Nj+1))

Wslot(s,K3)=WMAP(Zslot(s,K3))

Zint(s,K3)=Zslot(s,K3)+dx

Wint(s,K3)=WMAP(Zint(s,K3))

zex(s,K3)=Zslot(s,K3)-dx

Wext(s,K3)=WMAP(zex(s,K3))

enddo

endif

if ( slots(s,4)==ow ) then

do K3 =1, Nj

Zslot(s,K3)=slots(s,2)-slots(s,1)/2+K3*slots(s,1)/(Nj+1)+(0,1)*(slots(s,3))

Wslot(s,K3)=WMAP(Zslot(s,K3))

Zint(s,K3)=Zslot(s,K3)+(0,1)*dx

Wint(s,K3)=WMAP(Zint(s,K3))

zex(s,K3)=Zslot(s,K3)-(0,1)*dx

Wext(s,K3)=WMAP(zex(s,K3))

enddo

endif

if ( slots(s,4)==hw ) then

do K3 =1, Nj

Zslot(s,K3)=slots(s,2)-slots(s,1)/2+K3*slots(s,1)/(Nj+1)+(0,1)*(slots(s,3))

Wslot(s,K3)=WMAP(Zslot(s,K3))

Zint(s,K3)=Zslot(s,K3)-(0,1)*dx

Wint(s,K3)=WMAP(Zint(s,K3))

zex(s,K3)=Zslot(s,K3)+(0,1)*dx

Wext(s,K3)=WMAP(zex(s,K3))

enddo

endif

enddo

! WRITE (3, '(A50)') 'SLOT POINTS (BORDER) *** OK'

do s=1,Ns

aa=getparamA(slots(s,4))

bb=getparamB(slots(s,4))

if ((slots(s,4)==oh) .or. (slots(s,4)==wh) ) then

do k =1,Nk

do i3 =1, Nj

BFF=DBESJN(i3,aa*k*slots(s,1)/2)

CSSS=dcos(aa*k*slots(s,3)+Pi*i3/2)

BJCS(s,i3,k)=-i3*2*BFF*CSSS/k

enddo

enddo

else

do k =1,Nk

do i3 =1, Nj

BFF=DBESJN(i3,aa*k*slots(s,1)/2)

CSSS=dcos(aa*k*slots(s,2)+Pi*i3/2)

BJCS(s,i3,k)=-i3*2*BFF*CSSS/k

enddo

enddo

endif

enddo

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

! WRITE (3, *) 'START FLOODING MATRICES A,B [AX=B] FOR RODS POINTS'

call getaim(eps)

call DEVCRG(n,Aim,n,Zm,Iim2,n)

call DLINRG(n,aim,n,aim1,n)

call getaim(eps0)

aim2=matmul(aim1,aim)

call DEVCRG(n,Aim2,n,Zm2,Iim,n)

em=real(zm2)*eps

ksv=(em(2)-em(1))/(em(2)+em(1))

WRITE (3, *) dabs(ksv)

enddo !mm

if (draws==.true.) then

do i1=1,N

do m=1,N

WRITE (3, *)

WRITE (3, *)

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

WRITE (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)

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

enddo

enddo

! read(3,*)

do s = 1 , Ns

do r1 = 0 , WALLn

Zwall(s,r1)=WALLresW*r1

Wwall(s,r1)=WMAP(Zwall(s,r1))

enddo

enddo

do m = 1 , N

X(1:Nbig+Nj*Ns)=MX(1:Nbig+Nj*Ns,M)

do r1 = 0 , WALLn

do s1 = 1 , Ns

argl=(Wwall(s1,r1)-Wi(m))

dummyr=dlog(cdabs((Wwall(s1,r1)-dconjg(Wi(m)))/argl))

do I2 = 1 , N

do L2 = 1 , L

argl1=(Wwall(s1,r1)-Wi(I2))

argl2=dconjg(Wwall(s1,r1))-Wi(I2)

dummyc=(argl1)**(-L2) - (argl2)**(-L2)

dummyr=dummyr + dreal(dummyc)*X(IND2(I2,L2,1)) -dimag(dummyc)*X(IND2(I2,L2,2))

enddo

enddo

potwall(m,s1,r1,1)=dummyr

potwall(m,s1,r1,0)=r1*wallREsW

enddo

UU=0

do s = 1 , Ns

aa=getparamA(slots(s,4))

bb=getparamB(slots(s,4))

if (slots(s,4)==wh) then

do ii = 1 , Nj

do k = Nk,1,-1

UU=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)

enddo

enddo

endif

if (slots(s,4)==oh) then

do ii = 1 , Nj

do k = Nk,1,-1

UU=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)

enddo

enddo

endif

if (slots(s,4)==hw) then

do ii = 1 , Nj

do k = Nk,1,-1

UU=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)

enddo

enddo

endif

if (slots(s,4)==ow) then

do ii = 1 , Nj

if ((metod==1) ) then

do k = Nk,1,-1

addd=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=UU+addd

enddo

if ((metodnk==2)) then

k0=Nk+1

ARG10=aa*dreal(Zwall(s,r1))

ARG2=-aa*dImag(Zwall(s,r1))

ARG3=-2*aa*(bb-dImag(Zwall(s,r1)))

ARG4=-2*aa*bb

ARG7=aa*slots(s,1)/2

ARG9=-pI*Ii/2-pI/4

CALL dQdAGI(f2 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)

ADDD=-ADDD*Ii*2.8284271247461900976033774484194D0/sqrt(PI)

uu=uu+addd

endif

endif

if (metod==2) then

k0=0

arg1=ii

!ARG3=

ARG4=dImag(Zwall(s,r1))

ARG7=slots(s,1)/2

ARG9=pI*Ii/2

arg8=slots(s,2)-dreal(Zwall(s,r1))

CALL DQDAGI(ggggg,K0,INTERVAL,ERRABS,ERRREL,addd,ERREST)

uu=uu+X(IND3(ii,s))*addd*arg7/2

ENDIF

enddo

endif

psiwall(m,s,r1,1)=UU

psiwall(m,s,r1,0)=r1*wallREsW

enddo

enddo

enddo

sumwall(1:N,1:Ns,0:wallN,1)=potwall(1:N,1:Ns,0:wallN,1)+psiwall(1:N,1:Ns,0:wallN,1)

sumwall(1:N,1:Ns,0:wallN,0)=potwall(1:N,1:Ns,0:wallN,0)

do s1=1,Ns

do m=1,N

WRITE (3, *)

WRITE (3, *)

WRITE (3, *)'S= ',s1

WRITE (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)

enddo

enddo

endif

WRITE (3, *) CPSEC(),' SECONDS '

CLOSE (3)

CONTAINS

subroutine getaim(er)

real(8) er

do M =1, N

do I1 =1, N

do L1 =1, 2*L

B(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),')'

do I2 =1, N

do L2 =1, L

dummyc=(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)

A(IND1(I1,L1),IND2(I2,L2,1))= dReal(dummyc)

! WRITE (3, *) 'A(',IND1(I1,L1),IND2(I2,L2,1),')'

A(IND1(I1,L1),IND2(I2,L2,2))=-dImag(dummyc)

! WRITE (3, *) 'A(',IND1(I1,L1),IND2(I2,L2,2),')'

enddo

enddo

do s =1, Ns

aa=getparamA(slots(s,4))

bb=getparamB(slots(s,4))

if ( slots(s,4)==wh ) then

do i3 =1, Nj

PsiInt=0

do k =Nk,1,-1

addd=(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=PsiInt+addd

enddo

A(IND1(I1,L1),IND3(i3,s))=PsiInt

! WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'

enddo

endif

if ( slots(s,4)==oh ) then

do i3 =1, Nj

PsiInt=0

do k =Nk,1,-1

addd=(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=PsiInt+addd

enddo

A(IND1(I1,L1),IND3(i3,s))=PsiInt

! WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'

enddo

endif

if ( slots(s,4)==hw ) then

do i3 =1, Nj

PsiInt=0

! if ((metod==1) .or. (metod==2)) then

do k =Nk,1,-1

addd=(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=PsiInt+addd

enddo

! endif

! if (metod==2) then

!ENDIF

A(IND1(I1,L1),IND3(i3,s))=PsiInt

!WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'

enddo

endif

if ( slots(s,4)==ow ) then

do i3 =1, Nj

PsiInt=0

!if (metod==1) then

do k =Nk,1,-1

SINNN1=dsin(k*aa*dReal(Zil(I1,L1)))

EXP1=dexp(-k*aa*dImag(Zil(I1,L1)))

EXP2=1-dexp(-2*k*aa*(bb-dImag(Zil(I1,L1))))

EXP3=1-dexp(-2*k*aa*bb)

SINNN2=dsin(k*aa*dReal(Zil(I1,0)))

EXP4=dexp(-k*aa*dImag(Zil(I1,0)))

EXP5=1-dexp(-2*k*aa*(bb-dImag(Zil(I1,0))))

!EXP6=1-dexp(-2*k*aa*bb)

addd=(SINNN1*EXP1*EXP2-SINNN2*EXP4*EXP5)*BJCS(s,i3,k)/EXP3

PsiInt=PsiInt+addd

enddo

if ((metodnk==2)) then

k0=Nk+1

ARG10=aa*dReal(Zil(I1,L1))

ARG2=-aa*dImag(Zil(I1,L1))

ARG3=-2*aa*(bb-dImag(Zil(I1,L1)))

ARG4=aa*dReal(Zil(I1,0))

arg5=-aa*dImag(Zil(I1,0))

arg6=-2*aa*(bb-dImag(Zil(I1,0)))

ARG7=aa*slots(s,1)/2

arg8=-2*aa*bb

ARG9=-pI*I3/2-pI/4

CALL dQdAGI(f1 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)

ADDD=-ADDD*I3*2.8284271247461900976033774484194D0/sqrt(PI)

psiint=psiint+addd

endif

! endif

A(IND1(I1,L1),IND3(i3,s))=PsiInt

! WRITE (3, *) 'A(',IND1(I1,L1),IND3(i3,s),')'

enddo

endif

enddo

enddo

enddo

do s1 =1, Ns

do s2 =1, Ns

do i1 =1, Nj

do i2 =1, Nj

A(IND3(i1,s1),IND3(i2,s2))=0

enddo

enddo

enddo

enddo

! WRITE (3, *) 'START FLOODING MATRICES A,B [AX=B] FOR SLOTS POINTS'

do s =1, Ns

aa=getparamA(slots(s,4))

bb=getparamB(slots(s,4))

do i =1, Nj

argl=Wint(s,i)-Wi(M)

B(IND3(i,s))=dlog(cdabs((Wint(s,i)-dconjg(Wi(M)))/argl))

! WRITE (3, *) dummyr

do I2 =1, N

do L2 =1, L

dummyc=-((Wint(s,i)-Wi(I2))**(-L2)-(dconjg(Wint(s,i))-Wi(I2))**(-L2))

! WRITE (3, *) dummyc

A(IND3(i,s),IND2(I2,L2,1))= dreal(dummyc)

! WRITE (3, *) 'A(',IND3(i,s),IND2(I2,L2,1),')'

A(IND3(i,s),IND2(I2,L2,2))=-dimag(dummyc)

! WRITE (3, *) 'A(',IND3(i,s),IND2(I2,L2,2),')'

enddo

enddo

if ( slots(s,4)==wh ) then

do i3 =1, Nj

PsiInt=0

PsiExt=0

PsiSlot=0

PsiSlot=dsin(i3*dacos( 2*(dImag(Zslot(s,i))-slots(s,3))/slots(s,1)) )

do k =Nk,1,-1

addd=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=PsiInt+addd

addd=dsin(k*aa*dimag(zex(s,i)))*dexp(k*aa*(bb-dreal(zex(s,i))))*BJCS(s,i3,k)

PsiExt=PsiExt+addd

enddo

A(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

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

enddo

endif

if ( slots(s,4)==oh ) then

do i3 =1, Nj

PsiInt=0

PsiExt=0

PsiSlot=0

PsiSlot=dsin(i3*dacos( 2*(dimag(Zslot(s,i))-slots(s,3))/slots(s,1)) )

do k =Nk,1,-1

addd=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=PsiInt+addd

addd=dsin(k*aa*dimag(zex(s,i)))*dexp(k*aa*dreal(zex(s,i)))*BJCS(s,i3,k)

PsiExt=PsiExt+addd

enddo

A(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

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

enddo

endif

if ( slots(s,4)==hw ) then

do i3 =1, Nj

PsiInt=0

PsiExt=0

PsiSlot=0

PsiSlot=dsin(i3*dacos( 2*(dreal(Zslot(s,i))-slots(s,2))/slots(s,1)) )

do k =Nk,1,-1

addd=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=PsiInt+addd

addd=dsin(k*aa*dreal(zex(s,i)))*dexp(k*aa*(bb-dimag(zex(s,i))))*BJCS(s,i3,k)

PsiExt=PsiExt+addd

enddo

A(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

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

enddo

endif

if ( slots(s,4)==ow ) then

do i3 =1, Nj

PsiInt=0

PsiExt=0

PsiSlot=0

PsiSlot=dsin(i3*dacos( 2*(dreal(Zslot(s,i))-slots(s,2))/slots(s,1)) )

do k =Nk,1,-1

SINNN1=dsin(k*aa*dreal(Zint(s,i)))

EXP1=dexp(-k*aa*dImag(Zint(s,i)))

EXP2=1-dexp(-2*k*aa*(bb-dImag(Zint(s,i))))

EXP3=1-dexp(-2*k*aa*bb)

addd=SINNN1*EXP1*EXP2*BJCS(s,i3,k)/EXP3

PsiInt=PsiInt+addd

enddo

if ((metodnk==2)) then

k0=Nk+1

ARG10=aa*dReal(Zint(s,i))

ARG2=-aa*dImag(Zint(s,i))

ARG3=-2*aa*(bb-dImag(Zint(s,i)))

ARG4=0

arg5=0

arg6=0

ARG7=aa*slots(s,1)/2

arg8=-2*aa*bb

ARG9=-pI*I3/2-pI/4

CALL dQdAGI(f1 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)

ADDD=-ADDD*I3*2.8284271247461900976033774484194D0/sqrt(PI)

psiint=psiint+addd

endif

if ((metod==1)) then

do k =Nk,1,-1

SINNN2=dsin(k*aa*dreal(zex(s,i)))

EXP4=dexp(k*aa*dimag(zex(s,i)))

addd=SINNN2*EXP4*BJCS(s,i3,k)

PsiExt=PsiExt+addd

enddo

!WRITE (*, *) i3,psiext

!psiext=0

if ((metodnk==2)) then

k0=Nk+1

ARG10=aa*dreal(zex(s,i))

ARG2=-aa*dImag(zex(s,i))

ARG3=-2*aa*(bb-dImag(zex(s,i)))

ARG4=-2*aa*bb

ARG7=aa*slots(s,1)/2

ARG9=-pI*I3/2-pI/4

CALL dQdAGI(f2 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)

ADDD=-ADDD*I3*2.8284271247461900976033774484194D0/sqrt(PI)

psiext=psiext+addd

endif

endif

if (metod==2) then

k0=0

arg1=i3

!ARG3=dreal(zex(s,i))

ARG4=dimag(zex(s,i))

ARG7=slots(s,1)/2

ARG9=pI*I3/2

arg8=slots(s,2)-dreal(zex(s,i))

CALL DQDAGI(ggggg,K0,INTERVAL,ERRABS,ERRREL,addd,ERREST)

PsiExt=addd*arg7/2

WRITE (*, *) i3,psiext

ENDIF

A(IND3(i,s),IND3(i3,s))=(1+1/Er)*PsiSlot-PsiInt-(1/Er)*PsiExt

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

enddo

endif

enddo

enddo

! WRITE (3, *) 'solving linear equation******'

CALL DLSARG(Nbig+Nj*Ns,A,Nbig+Nj*Ns,B,1,X)

! WRITE (3, *) 'Ok'

CALL 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)

MX(1:Nbig+Nj*Ns,M)=X(1:Nbig+Nj*Ns)

enddo

do m=1,N

do r1=0,PHIn

Zrad(m,r1)=Zi(m)+Ri(m)*exp(2*(0,1)*r1*Pi/PHIn)

Wrad(m,r1)=WMAP(Zrad(m,r1))

enddo

enddo

do 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)

X(1:Nbig+Nj*Ns)=MX(1:Nbig+Nj*Ns,M)

do i1 = 1 , N

do r1 = 0 , PHIn

argl=Wrad(i1,r1)-Wi(m)

dummyr=dlog(CDABS((Wrad(i1,r1)-DCONJG(Wi(m)))/argl))

do I2 = 1 , N

do L2 = 1 , L

argl1=Wrad(i1,r1)-Wi(I2)

argl2=DCONJG(Wrad(i1,r1))-Wi(I2)

dummyc=argl1**(-L2) - argl2**(-L2)

dummyr=dummyr + dreal(dummyc)*X(IND2(I2,L2,1))-dimag(dummyc)*X(IND2(I2,L2,2))

enddo

enddo

potrad(m,i1,r1,1)=dummyr

potrad(m,i1,r1,0)=r1*PHIREs

UU=0

do s = 1 , Ns

aa=getparamA(slots(s,4))

bb=getparamB(slots(s,4))

if (slots(s,4)==wh) then

do ii = 1 , Nj

do k = Nk,1,-1

addd=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=UU+addd

enddo

enddo

endif

if (slots(s,4)==oh) then

do ii = 1 , Nj

do k = Nk,1,-1

addd=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=UU+addd

enddo

enddo

endif

if (slots(s,4)==hw) then

do ii = 1 , Nj

do k = Nk,1,-1

addd=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=UU+addd

enddo

enddo

endif

if (slots(s,4)==ow) then

do ii = 1 , Nj

do k = Nk,1,-1

addd=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=UU+addd

enddo

if ((metodnk==2)) then

k0=Nk+1

ARG10=aa*dReal(Zrad(i1,r1))

ARG2=-aa*dImag(Zrad(i1,r1))

ARG3=-2*aa*(bb-dImag(Zrad(i1,r1)))

ARG4=0

arg5=0

arg6=0

ARG7=aa*slots(s,1)/2

arg8=-2*aa*bb

ARG9=-pI*Ii/2-pI/4

CALL dQdAGI(f1 , K0, INTERV, ERRABS, ERRREL, addd, ERREST)

ADDD=-ADDD*Ii*2.8284271247461900976033774484194D0/sqrt(PI)

uu=uu+addd

endif

enddo

endif

enddo

psirad(m,i1,r1,1)=UU

psirad(m,i1,r1,0)=r1*PHIREs

enddo

enddo

enddo

sumrad(1:N,1:N,0:PHIN,1)=potrad(1:N,1:N,0:PHIN,1)+psirad(1:N,1:N,0:PHIN,1)

sumrad(1:N,1:N,0:PHIN,0)=potrad(1:N,1:N,0:PHIN,0)

Aim=0

do r1=0,PHIn

Aim(1:N,1:N)=Aim(1:N,1:N)+sumrad(1:N,1:N,r1,1)/(1+Phin)

enddo

end subroutine getaim

INTEGER FUNCTION IND1(I1,L1)

INTEGER,INTENT(IN):: I1,L1 ! I1 - номер проводника; L1 - порядок мультиполя

IND1=(I1-1)*2*L+L1

END FUNCTION IND1

INTEGER FUNCTION IND2(I2,L2,M2)

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

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

END FUNCTION IND2

INTEGER FUNCTION IND3(j3,s3)

INTEGER,INTENT(IN):: j3,s3

IND3=2*L*N+Nj*(s3-1)+j3

END FUNCTION IND3

real(8) FUNCTION getparamA(s)

real(8),INTENT(IN)::s

if ((s==wh) .or. (s==oh)) then

getparamA=Pi/2/Hb

else

getparamA=Pi/Wb

end if

END FUNCTION getparamA

real(8) FUNCTION getparamB(s)

real(8),INTENT(IN)::s

if ( (s==wh) .or. (s==oh) ) then

getParamB=Wb

else

getParamB=2*Hb

endif

END FUNCTION getparamB

COMPLEX(8) FUNCTION WMAP(Z) ! Отображение на верхнюю полуплоскость (w=u+iv), v>=0

COMPLEX(8),INTENT(IN):: Z ! z=x+iy - комплексная координата в прямоугольнике (0<x<Wb); (0<y<2Hb)

COMPLEX(8) ZEJDN

EXTERNAL ZEJDN ! IMSL: dn(z,m)

WMAP=-ZEJDN(Z*K1/Hb,AM)/DSQRT(1-AM)

END FUNCTION

REAL(8) FUNCTION AMELL(K1) ! Корень K(m)/K(1-m)=Wb/Hb, где m=AMELL, K(m) - полный эллиптический интеграл Лежандра I рода

REAL(8),INTENT(OUT):: K1 ! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода, m=k^2 - параметр интеграла

REAL(8) DELK,R,X1,X2,X,Y1,Y2,Y

EXTERNAL DELK ! IMSL: K(m)

R=Wb/Hb

X1=0; X2=1; X=0.5

DO

Y=DELK(X)/DELK(1-X)

IF(Y==R) THEN

EXIT

ELSEIF(Y>R) THEN

X2=X; Y2=Y

ELSEIF(Y<R) THEN

X1=X; Y1=Y

ENDIF

IF(X1==0 .OR. X2==1) THEN

X=(X1+X2)/2

ELSE

X=X1 + (X2-X1)/(Y2-Y1)*(R-Y1)

IF(X<X1+(X2-X1)/10) X=X1+(X2-X1)/10

IF(X>X2-(X2-X1)/10) X=X2-(X2-X1)/10

ENDIF

IF((X2-X1)/X <= 4*EPSILON(X)) EXIT

ENDDO

K1=DELK(1-X)

AMELL=X

END FUNCTION AMELL

END

real(8) function ggggg(x)

REAL(8) X ,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

integer(4) a,b,ARG1

COMMON ARG1

COMMON ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

arg2=DSIN(arg9+(arg8)*X)

a=arg1+1

arg3=DBESJN(arg1+1,X*ARG7)

b=arg1-1

arg5=DBESJN(arg1-1,X*ARG7)

arg6=dexp(arg4*x)

ggggg = arg2*(arg3+arg5)*arg6

!WRITE (3, *) x,ggggg

end function ggggg

real(8) function f1(x)

REAL(8) X,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

integer(4) ARG1

COMMON ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

F1 = (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, *) f1

END function f1

real(8) function f2(x)

REAL(8) X,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9

integer(4) ARG1

COMMON ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

f2 = DSIN(X*ARG10)*DEXP(X*ARG2)*(1-DEXP(X*ARG3))*DCOS(ARG7*X-ARG9)/((1-DEXP(ARG4*X))*X*DSQRT(X*ARG7))

! WRITE (3, *) f2

end function f2

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

program ANROD5

USE PORTLIB

USE dfIMSL

INTEGER(4),parameter:: N=2,L=10,Nbig=N*2*L,Ns=1

REAL(8):: h2,Ri(N),Eps,Eps2(Ns),Rs(Ns) ,ERRABS,ERRREL,ERREST

logical(1):: draws

COMPLEX(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)

REAL(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

COMPLEX(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)

integer(4) k,i1,i2,i3,s,s1,s2,k3,ii,r1,mm

REAL(8),PARAMETER::Pi=3.14159265358979323846264338327950288419716939937510d0, Zc=376.73035d0 ,vellight=2.99792E+8

REAL(8),PARAMETER::Wb=20,Hb=5,eps0=1

integer,parameter::PHIn=100 ,WALLN=500

COMPLEX(8),parameter::PHIREs=2*Pi/PHIn ,WALLresW=Wb/WALLn,WALLresH=2*hb/WALLn

COMPLEX(8)::Zrad(N,0:PHIn),Wrad(N,0:PHIn),argl1,argl2,Zsss(Ns,0:phin),Wsss(Ns,0:phin)

REAL(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)

REAL(8)::SINNN1,SINNN2,EXP1,EXP2,EXP3,EXP4,EXP5,Aim(n,n),Aim1(n,n),Aim2(n,n),em(n),K0,argl

COMPLEX(8) ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,arg10

OPEN (3, FILE = 'output.TXT')

CALL UMACH (-2, 3)

CALL PGOPT (-1, 255)

CALL ERSET (0, 0, 0)

WRITE (3, *) CPSEC(),' SECONDS '

h2=2*Hb

Zi(1)=(5,5)

Zi(2)=(15,5)

Zs(1)=(10,6)

Eps=100

eps2(1)=1

Ri(1)=1

Ri(2)=1

Rs(1)=1

dx=0.0001

AM=AMELL(K1)

do I1 =1, N

Wi(I1)=WMAP( Zi(I1) )

do L1=0,2*L

Zil(I1,L1)=Zi(I1)+Ri(I1)*exp((0,2)*L1*Pi/(2*L+1))

Wil(I1,L1)=WMAP(Zil(I1,L1))

enddo

enddo

do I1 =1, Ns

Ws(I1)=WMAP( Zs(I1) )

enddo

do s =1, Ns

do K3 =1, 2*L

Zsb(s,K3)=zs(s)+Rs(s)*exp((0,2)*k3*Pi/(2*L+1))

Wsb(s,K3)=WMAP(Zsb(s,K3))

Zsi(s,K3)=zs(s)+(Rs(s)-dx)*exp((0,2)*k3*Pi/(2*L+1))

Wsi(s,K3)=WMAP(Zsi(s,K3))

zse(s,K3)=zs(s)+(Rs(s)+dx)*exp((0,2)*k3*Pi/(2*L+1))

Wse(s,K3)=WMAP(zse(s,K3))

enddo

enddo

call getaim(eps)

call DEVCRG(n,Aim,n,Zm,Iim2,n)

call DLINRG(n,aim,n,aim1,n)

call getaim(eps0)

aim2=matmul(aim1,aim)

call DEVCRG(n,Aim2,n,Zm2,Iim,n)

em=real(zm2)*eps

ksv=(em(2)-em(1))/(em(2)+em(1))

WRITE (3, *) dabs(ksv)

WRITE (3, *) CPSEC(),' SECONDS '

CLOSE (3)

CONTAINS

subroutine getaim(er)

real(8) er

do M =1, N

a=0

b=0

x=0

do I1 =1, N

do L1 =1, 2*L

B(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))))

do I2 =1, N

do L2 =1, L

dummyc=(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)

A(IND1(I1,L1),IND2(I2,L2,1))= dReal(dummyc)

A(IND1(I1,L1),IND2(I2,L2,2))=-dImag(dummyc)

enddo

enddo

do I2 =1, Ns

do L2 =1, L

dummyc=(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)

A(IND1(I1,L1),nbig+IND2(I2,L2,1))= dReal(dummyc)

A(IND1(I1,L1),nbig+IND2(I2,L2,2))=-dImag(dummyc)

A(IND1(I1,L1),nbig+ns*L*2+IND2(I2,L2,1))= 0

A(IND1(I1,L1),nbig+ns*L*2+IND2(I2,L2,2))=0

enddo

enddo

enddo

enddo

do I1 =1, Ns

do L1 =1, 2*L

B(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))) )) )

do I2 =1, N

do L2 =1, L

dummyc=((Wsb(I1,L1)-Wi(I2))**(-L2)-(conjg(Wsb(I1,L1))-Wi(I2))**(-L2))

dummyc=dummyc-((Wsi(I1,L1)-Wi(I2))**(-L2)-(conjg(Wsi(I1,L1))-Wi(I2))**(-L2))

dummyc=dummyc-(er/eps2(i1))*((Wse(I1,L1)-Wi(I2))**(-L2)-(conjg(Wse(I1,L1))-Wi(I2))**(-L2))

dummyc=dummyc+(er/eps2(i1))*((Wsb(I1,L1)-Wi(I2))**(-L2)-(conjg(Wsb(I1,L1))-Wi(I2))**(-L2))

A(nbig+IND1(I1,L1),IND2(I2,L2,1))= dReal(dummyc)

A(nbig+IND1(I1,L1),IND2(I2,L2,2))=-dImag(dummyc)

enddo

enddo

do I2 =1, Ns

do L2 =1, L

dummyc=((Wsb(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsb(I1,L1))-Ws(I2))**(-L2))

dummyc=dummyc-(er/eps2(i1))*((Wse(I1,L1)-Ws(I2))**(-L2)-(conjg(Wse(I1,L1))-Ws(I2))**(-L2))

dummyc=dummyc+(er/eps2(i1))*((Wsb(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsb(I1,L1))-Ws(I2))**(-L2))

if (i1/=i2) then

dummyc=dummyc-((Wsi(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsi(I1,L1))-Ws(I2))**(-L2))

endif

A(nbig+IND1(I1,L1),nbig+IND2(I2,L2,1))= dReal(dummyc)

A(nbig+IND1(I1,L1),nbig+IND2(I2,L2,2))=-dImag(dummyc)

if (i1==i2) then

dummyc=-((Zsi(I1,L1)-Zs(I2)) )**(L2)

A(nbig+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,1))= dReal(dummyc)

A(nbig+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,2))=-dImag(dummyc)

endif

enddo

enddo

do I2 =1, Ns

do L2 =1, L

if (i1==i2) then

dummyc=((Wsb(I1,L1)-Ws(I2))**(-L2)-(conjg(Wsb(I1,L1))-Ws(I2))**(-L2))

A(nbig+2*L*Ns+IND1(I1,L1),nbig+IND2(I2,L2,1))= dReal(dummyc)

A(nbig+2*L*Ns+IND1(I1,L1),nbig+IND2(I2,L2,2))=-dImag(dummyc)

dummyc=-(zsb(I1,L1)-zs(I2))**(L2)

A(nbig+2*L*Ns+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,1))= dReal(dummyc)

A(nbig+2*L*Ns+IND1(I1,L1),nbig+2*L*Ns+IND2(I2,L2,2))=-dImag(dummyc)

endif

enddo

enddo

enddo

enddo

CALL DLSARG(Nbig+4*L*Ns,A,Nbig+4*L*Ns,B,1,X)

MX(1:Nbig+4*L*Ns,M)=X(1:Nbig+4*L*Ns)

enddo

do m=1,N

do r1=0,PHIn

Zrad(m,r1)=Zi(m)+Ri(m)*exp(2*(0,1)*r1*Pi/PHIn)

Wrad(m,r1)=WMAP(Zrad(m,r1))

enddo

enddo

do m = 1 , N

X(1:Nbig+4*L*Ns)=MX(1:Nbig+4*L*Ns,M)

do i1 = 1 , n

do r1 = 0 , phin

dummyr=0

dz=Zrad(i1,r1)

dw=Wrad(i1,r1)

arg1=dw-Wi(m)

argl=CDABS((dw-DCONJG(Wi(m)))/arg1)

dummyr=dlog(argl)

do I2 = 1 , N

do L2 = 1 , L

arg1=dw-Wi(I2)

arg2=DCONJG(dw)-Wi(I2)

dummyc=arg1**(-L2) - arg2**(-L2)

dummyr=dummyr + dreal(dummyc)*X(IND2(I2,L2,1))-dimag(dummyc)*X(IND2(I2,L2,2))

enddo

enddo

do I2 = 1 , Ns

do L2 = 1 , L

arg1=dw-Ws(I2)

arg2=DCONJG(dw)-Ws(I2)

dummyc=arg1**(-L2) - arg2**(-L2)

dummyr=dummyr + dreal(dummyc)*X(nbig+IND2(I2,L2,1))-dimag(dummyc)*X(nbig+IND2(I2,L2,2))

enddo

enddo

potrad(m,i1,r1,1)=dummyr

potrad(m,i1,r1,0)=r1*PHIREs

enddo

enddo

enddo

enddo

Aim=0

do r1=0,PHIn

Aim(1:N,1:N)=Aim(1:N,1:N)+potrad(1:N,1:N,r1,1)/(1+Phin)

enddo

end subroutine getaim

INTEGER FUNCTION IND1(I1,L1)

INTEGER,INTENT(IN):: I1,L1 ! I1 - номер проводника; L1 - порядок мультиполя

IND1=(I1-1)*2*L+L1

END FUNCTION IND1

INTEGER FUNCTION IND2(I2,L2,M2)

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

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

END FUNCTION IND2

COMPLEX(8) FUNCTION WMAP(Z) ! Отображение на верхнюю полуплоскость (w=u+iv), v>=0

COMPLEX(8),INTENT(IN):: Z ! z=x+iy - комплексная координата в прямоугольнике (0<x<Wb); (0<y<2Hb)

COMPLEX(8) ZEJDN

EXTERNAL ZEJDN ! IMSL: dn(z,m)

WMAP=-ZEJDN(Z*K1/Hb,AM)/DSQRT(1-AM)

END FUNCTION

REAL(8) FUNCTION AMELL(K1) ! Корень K(m)/K(1-m)=Wb/Hb, где m=AMELL, K(m) - полный эллиптический интеграл Лежандра I рода

REAL(8),INTENT(OUT):: K1 ! K1=K(1-m) - связанный полный эллиптический интеграл Лежандра I рода, m=k^2 - параметр интеграла

REAL(8) DELK,R,X1,X2,X,Y1,Y2,Y

EXTERNAL DELK ! IMSL: K(m)

R=Wb/Hb

X1=0; X2=1; X=0.5

DO

Y=DELK(X)/DELK(1-X)

IF(Y==R) THEN

EXIT

ELSEIF(Y>R) THEN

X2=X; Y2=Y

ELSEIF(Y<R) THEN

X1=X; Y1=Y

ENDIF

IF(X1==0 .OR. X2==1) THEN

X=(X1+X2)/2

ELSE

X=X1 + (X2-X1)/(Y2-Y1)*(R-Y1)

IF(X<X1+(X2-X1)/10) X=X1+(X2-X1)/10

IF(X>X2-(X2-X1)/10) X=X2-(X2-X1)/10

ENDIF

IF((X2-X1)/X <= 4*EPSILON(X)) EXIT

ENDDO

K1=DELK(1-X)

AMELL=X

END FUNCTION AMELL

END

Размещено на Allbest.ru