Оптимизация двухдипольной излучающей системы
Министерство образования и науки Российской Федерации НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ Кафедра РПиРПУ Пояснительная записка к курсовой работе по дисциплине «Устройства СВЧ и антенны». Построение и анализ диаграмм направленности двухдипольной излучающей системы Плоскость магнитного вектора: Рис. 6. Диаграмм направленности для всей плоскости XOY в полярной системе координат… Читать ещё >
Оптимизация двухдипольной излучающей системы (реферат, курсовая, диплом, контрольная)
Министерство образования и науки Российской Федерации НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ Кафедра РПиРПУ Пояснительная записка к курсовой работе по дисциплине «Устройства СВЧ и антенны»
Новосибирск, 2014
1.Оптимизация двухдипольной излучающей системы Задание. Найти оптимальные размеры двухдипольной излучающей системы с нумерацией диполей согласно рисунку 1 при условии, что 2-й диполь является пассивным токонесущим (т.е. короткозамкнутым), частота сигнала равна 870 МГц, волновое сопротивление питающего коаксиального кабеля, радиус проводников диполей равен 1 мм.
Рис. 1. Двухдипольная излучающая система Выполнение задания. По модифицированной Фортран-программе optimizac2dipoln. ex определяется оптимальные значения размеров. В исходных данных все диполи вначале берутся полуволновыми. Для начало вычислим длину волны:
Итак с клавиатуры вводятся следующие значения:
После начала вычисления на экран выводятся текущие значения целевой функции, массива трех неизвестных переменных и другая связанная с процессом оптимизации информация. Затем на экран выводятся следующие значения:
Полученные результаты означают что:
оптимальная длина возбудителя 2l1 = 2*78.69 мм;
оптимальная длина короткозамкнутого рефлектора 2l1 = 2*99,45 мм;
оптимальное расстояние между ними d = 60,56 мм:
Величина R01 показывает значение входного сопротивления возбудителя системы. Видим, что оно практически равно. Величина mod и fm (радианы) показывает значение модуля и фазы. Они необходимы для построения и анализа диаграмм направленности проектируемой излучающей системы в плоскости как электрического, так и магнитного вектора.
Построение и анализ диаграмм направленности двухдипольной излучающей системы Плоскость магнитного вектора:
Построим диаграммы направленности в пакете mathcad 14:
оптимальный излучающая система Рис. 2. Нормированная диаграмм направленности в плоскости H.
Рис. 3. Диаграмм направленности в полярной системе координат,
в плоскости H.
Плоскость электрического вектора:
Угол и изменяется от 0 до 180.
mp[цcut = 90] = 60,07
mp[цcut = 270] = 21,90
FE (и = 90) = mp[цcut = 90]/mp[цcut = 270] = 60,07/21,90= 0.332
Рис. 4 Нормированная диаграмма направленности в плоскости Е 270
Рис. 5 Нормированная диаграмма направленности Е 90
Рис. 6. Диаграмм направленности для всей плоскости XOY в полярной системе координат, в плоскости E.
Приложение
program optimizac2dipoln
c calculation lengths & distance between TWO DIPOLES
c when the cable impedance «Rcab» is prescribed
dimension x1(9), h (9,10), x (9), x0(9)
real la, mod
write (*,*) '******* Modified Powell_s method *******'
write (*,*) '- - - - - The steepest-descent method — - - - - -'
write (*,*) '- - - - The conjugate-gradient method — - - - -'
write (*,*) '. .. .. optimization of TWO DIPOLES.. .. .. .'
write (*,*) 'sapros wire radius & lambda, both in «mm» '
read (*,*) wr, la
write (*,*) ' wr=', wr, ' la=', la
write (*,*) 'sapros cable impedance [Ohms] '
read (*,*) Rcab
write (*,*) ' Rcab=', Rcab
write (*,*) 'sapros number of variable N'
read (*,*) n
write (*,*) 'N=', n
write (*,*) 'sapros tochnostei: E1, E2'
read (*,*) E1, E2
write (*,*) ' E1=', E1, ' E2=', E2
write (*,*) 'sapros extremum of celewaja-function FH'
read (*,*) fh
write (*,*) ' FH=', fh
7 format (e12.4)
pi = 4.*atan (1.)
write (*, 8)
8 format (5x,'Sapros initial point/array x1, a11 in «mm» ')
read (*, 9) (x1(i), i=1,n)
9 format (9e12.4)
t=1.618
j=1
it=0
il=0
ib=1
iw=0
59 do 12 ig=ib, n
do 12 i=1,n
if (iw.eq.0) goto 10
h (i, ig)=h (i, ig+1)
goto 12
10 if (i.eq.ig) goto 11
h (i, ig)=0
goto 12
11 h (i, ig)=1
12 continue
if (i1.ne.0) goto 50
iq=0
do 112 i=1,n
112 x (i)=x1(i)
call cel2di1(wr, la, Rcab, n, x, cf)
cf0 = cf
write (*, 80)
80 format (3x, 'Goal function in initial point/array CF0= ')
write (*, 7) cf0
read (*,*)
write (*, 81)
81 format (/, 1x,'Table 1',/, 5x,'IT', 5x,'J', 5x,'IQ', 5x,'FO', 5x,'X1')
fo=cf0
i1=1
50 f11=fo
wl=0
is2=0
it=it+1
do 70 j=1,n+1
do 13 i=1,n
x1(i)=x (i)
if (j.ne.1) goto 13
x0(i)=x (i)
ww=x (i)
13 continue
xx1=0
do 14 ii=1,n
14 xx1=xx1+x1(ii)**2
u=e1*sqrt (xx1)
ir=-1
no1=0
no2=0
no3=0
q=u
b1=fo
b=0
52 iq=iq+1
do 17 i=1,n
17 x (i)=x1(i)+q*h (i, j)
call cel2di1(wr, la, Rcab, n, x, cf)
cf0 = cf
18 if (no2.eq.1) goto 19
a1=cf0
a=q
goto 53
19 if (no3.eq.1) goto 20
f2=cf0
goto 54
20 f1=cf0
goto 54
53 if (a1.le.b1) goto 21
if (no1.ne.0) goto 22
a=0
b=q
c1=a1
a1=b1
b1=c1
ir=-1
goto 53
21 c=b
b1=a1
b=a
q=t*(t*b-c)
ir=ir+1
no1=1
goto 52
22 f1=b1
q1=b
no2=1
55 no3=0
q2=c+(a-c)/t
q=q2
goto 52
57 no3=1
q1=c+(a-c)/t**2
q=q1
goto 52
54 if (abs (f2-f1).le.e2*abs (fo)) goto 62
if (f2.gt.f1) goto 23
c=q1
q1=q2
f1=f2
goto 55
23 a=q2
q2=q1
f2=f1
goto 57
62 f=cf0
write (*, 82) it, j, iq, f,(x (i), i=1,n)
82 format (5x, I3,5x, I2,5x, I4,5x, e12.4,5x, 9e12.4,/)
56 is2=is2+ir
if (j.ge.n+1) goto 24
dt=fo-f
fo=f
if (dt.le.wl) goto 24
wl=dt
ib=j
24 if (j.ne.n) goto 70
yy=0
do 25 i=1,n
25 yy=yy+(x (i)-x0(i))**2
y=sqrt (yy)
if (cf0*(fh-cf0)*abs (y).ge.0) goto 58
do 63 i=1,n
63 h (i, n+1)=(x (i)-x0(i))/y
70 continue
if (it.gt.3) goto 91
write (*, 72)
72 format (2x,'Number of iteration IT= ', i3)
read (*,*)
write (*, 73)
73 format (2x,'Matrix of search directions H= ')
do 74 ii1=1,n
74 write (*, 75) (h (ii1,j1), j1=1,n+1)
75 format (1x, 9e12.4)
read (*,*)
91 if (f.lt.f11) goto 28
58 write (*, 26)
26 format (2x,'The values of variables X1= ')
write (*, 9)(x (i), i=1,n)
read (*,*)
write (*, 27)
27 format (2x,'The value of goal function CF0= ')
write (*, 9) cf0
read (*,*)
c ———- Calculation of related values ——————;
call rinxin (x (1), wr, la, r11, x11)
call rinxin (x (2), wr, la, r22, x22)
call r12×12(x (1), x (2), x (3), la, r12, x12)
bbb1 = r12**2*r22/(r22**2+x22**2)
bbb2 = r22*x12**2/(r22**2+x22**2)
bbb3 = 2.*r12*x12*x22/(r22**2+x22**2)
r01 = r11-bbb1+bbb2-bbb3
mod = sqrt ((r12**2+x12**2)/(r22**2+x22**2))
fm = pi+atan (x12/r12)-atan (x22/r22)
write (*,*) 'Related values'
write (*,*) ' R01=', r01,' mod=', mod,' fm=', fm
read (*,*)
stop
28 if (it.eq.1) goto 30
if (is2.gt.1) goto 29
e1=e1/10
goto 30
29 if (0.5*is2/(n+1).le.1) goto 30
e1=e1*0.5*is2/(n+1)
30 if (q.lt.0) goto 31
if (4*wl*(fo-f).lt.(f11-fo-wl)**2) goto 31
iw=1
fo=f
goto 59
31 fo=f
goto 50
stop
end
subroutine cel2di1(wr, la, Rcab, n, x, cf)
c Goal function abs (X01) & abs (R01-Rcab)
c for TWO dipoles when the second dipole is passiv
c Both the dipoles have various lengths: L1 and L2
dimension x (n)
real l1, l2,la
l1 = x (1)
l2 = x (2)
d = x (3)
call rinxin (l1,wr, la, r11, x11)
call rinxin (l2,wr, la, r22, x22)
call r12×12(l1,l2,d, la, r12, x12)
den = r22**2+x22**2
a1 = 2.*r12*x12*r22/den
a2 = r12**2*x22/den
a3 = x12**2*x22/den
b1 = r12**2*r22/den
b2 = r22*x12**2/den
b3 = 2.*r12*x12*x22/den
R01 = r11-b1+b2-b3
cf = abs (x11-a1+a2-a3)+abs (R01-Rcab)
return
end
subroutine rinxin (l, wr, la, Rin, Xin)
cimpedances of classical dipole by method of inducted
celectomotive force when the wire radius is take into account
cThe half-lenght «l» is introduced by input parameters
cThis is by copy of 'dipoself'
real la, l
pi = 4.*atan (1.)
ak = 2.*pi/la
zh = l/1000.
aim1 = -sin (ak*l)**2/l
aim2 = -sin (ak*l)**2/l
a3 = 2.*cos (ak*l)
aim3 = -a3*ak*sin (ak*l)
do 100 i=1,999,2
aim1 = aim1−4.*sin (ak*(l-i*zh))*sin (ak*(l-i*zh))/(l-i*zh)
aim2 = aim2−4.*sin (ak*(l+i*zh))*sin (ak*(l-i*zh))/(l+i*zh)
100 aim3 = aim3−4.*a3*sin (ak*i*zh)*sin (ak*(l-i*zh))/(i*zh)
do 101 i=2,999,2
aim1 = aim1−2.*sin (ak*(l-i*zh))*sin (ak*(l-i*zh))/(l-i*zh)
aim2 = aim2−2.*sin (ak*(l+i*zh))*sin (ak*(l-i*zh))/(l+i*zh)
101 aim3 = aim3−2.*a3*sin (ak*i*zh)*sin (ak*(l-i*zh))/(i*zh)
Ra = -20.*zh*(aim1+aim2-aim3)
r10 = sqrt (wr**2+l**2)
aim11 = -sin (ak*r10)*sin (ak*l)/r10
r20 = sqrt (wr**2+l**2)
aim22 = -sin (ak*r20)*sin (ak*l)/r20
r30 = wr
aim33 = -sin (ak*r30)*a3*sin (ak*l)/r30
do 110 i=1,999,2
r1 = sqtr (wr**2+(i*zh-l)**2)
r2 = sqtr (wr**2+(i*zh+l)**2)
r3 = sqtr (wr**2+(i*zh)**2)
aim11 = aim11−4.*sin (ak*r1)*sin (ak*(l-i*zh))/r1
aim22 = aim22−4.*sin (ak*r2)*sin (ak*(l-i*zh))/r2
110 aim33 = aim33−4.*sin (ak*r3)*a3*sin (ak*(l-i*zh))/r3
do 111 i=2,999,2
r1 = sqtr (wr**2+(i*zh-l)**2)
r2 = sqtr (wr**2+(i*zh+l)**2)
r3 = sqtr (wr**2+(i*zh)**2)
aim11 = aim11−2.*sin (ak*r1)*sin (ak*(l-i*zh))/r1
aim22 = aim22−2.*sin (ak*r2)*sin (ak*(l-i*zh))/r2
111 aim33 = aim33−2.*sin (ak*r3)*a3*sin (ak*(l-i*zh))/r3
Raa = -20.*zh*(aim11+aim22-aim33)
re11 = cos (ak*r10)*sin (ak*l)/r10
re22 = cos (ak*r20)*sin (ak*l)/r20
re33 = cos (ak*r30)*a3*sin (ak*l)/r30
c ————- Calculation of image part of Za —————-;
c initial r10, r20, r30 are all the same as previous
do 120 i=1,999,2
r1 = sqtr (wr**2+(i*zh-l)**2)
r2 = sqtr (wr**2+(i*zh+l)**2)
r3 = sqtr (wr**2+(i*zh)**2)
re11 = re11+4.*cos (ak*r1)*sin (ak*(l-i*zh))/r1
re22 = re22+4.*cos (ak*r2)*sin (ak*(l-i*zh))/r2
120 re33 = re33+4.*cos (ak*r3)*a3*sin (ak*(l-i*zh))/r3
do 121 i=2,999,2
r1 = sqtr (wr**2+(i*zh-l)**2)
r2 = sqtr (wr**2+(i*zh+l)**2)
r3 = sqtr (wr**2+(i*zh)**2)
re11 = re11+2.*cos (ak*r1)*sin (ak*(l-i*zh))/r1
re22 = re22+2.*cos (ak*r2)*sin (ak*(l-i*zh))/r2
121 re33 = re33+2.*cos (ak*r3)*a3*sin (ak*(l-i*zh))/r3
Xaa = 20.*zh*(re11+re22-re33)
c ————- end calculation Xaa —————————-;
c ———— Input resistance & reactance ———;
c (referred to at the current at the input terminals)
Rin = Raa/(sin (ak*l)**2)
Xin = Xaa/(sin (ak*l)**2)
return
end
subroutine r12×12(l1,l2,d, la, R12, X12)
c Mutual impedance of TWO clasical dipoles by method of
c induced electromotive force when the wire radius is NULL
c Both the dipoles have various lengths: L1 and L2
c——— Z12 Z12 Z12 Z12 Z12 Z12 Z12——-;
real la, l1, l2
pi = 4.*atan (1.)
ak = 2.*pi/la
zh = l1/1000.
a3 = 2.*cos (ak*l2)
r10 = sqrt (d**2+l2**2)
aim11 = -sin (ak*r10)*sin (ak*l1)/r10
r20 = sqrt (d**2+l2**2)
aim22 = -sin (ak*r20)*sin (ak*l1)/r20
r30 = d
aim33 = -sin (ak*r30)*a3*sin (ak*l1)/r30
do 110 i=1,999,2
r1 = sqrt (d**2+(i*zh-l2)**2)
r2 = sqrt (d**2+(i*zh+l2)**2)
r3 = sqrt (d**2+(i*zh)**2)
aim11 = aim11−4.*sin (ak*r1)*sin (ak*(l1-i*zh))/r1
aim22 = aim22−4.*sin (ak*r2)*sin (ak*(l1-i*zh))/r2
110 aim33 = aim33−4.*sin (ak*r3)*a3*sin (ak*(l1-i*zh))/r3
do 111 i=2,999,2
r1 = sqrt (d**2+(i*zh-l2)**2)
r2 = sqrt (d**2+(i*zh+l2)**2)
r3 = sqrt (d**2+(i*zh)**2)
aim11 = aim11−2.*sin (ak*r1)*sin (ak*(l1-i*zh))/r1
aim22 = aim22−2.*sin (ak*r2)*sin (ak*(l1-i*zh))/r2
111 aim33 = aim33−2.*sin (ak*r3)*a3*sin (ak*(l1-i*zh))/r3
Rm = -20.*zh*(aim11+aim22-aim33)
re11 = cos (ak*r10)*sin (ak*l1)/r10
re22 = cos (ak*r20)*sin (ak*l1)/r20
re33 = cos (ak*r30)*a3*sin (ak*l1)/r30
c ————- Calculation of image part of Zm —————-;
c initial r10, r20, r30 are all the same as previous
do 120 i=1,999,2
r1 = sqrt (d**2+(i*zh-l2)**2)
r2 = sqrt (d**2+(i*zh+l2)**2)
r3 = sqrt (d**2+(i*zh)**2)
re11 = re11+4.*cos (ak*r1)*sin (ak*(l1-i*zh))/r1
re22 = re22+4.*cos (ak*r2)*sin (ak*(l1-i*zh))/r2
120 re33 = re33+4.*cos (ak*r3)*a3*sin (ak*(l1-i*zh))/r3
do 121 i=2,999,2
r1 = sqrt (d**2+(i*zh-l2)**2)
r2 = sqrt (d**2+(i*zh+l2)**2)
r3 = sqrt (d**2+(i*zh)**2)
re11 = re11+2.*cos (ak*r1)*sin (ak*(l1-i*zh))/r1
re22 = re22+2.*cos (ak*r2)*sin (ak*(l1-i*zh))/r2
121 re33 = re33+2.*cos (ak*r3)*a3*sin (ak*(l1-i*zh))/r3
Xm = 20.*zh*(re11+re22-re33)
c ————- end calculation Xm —————————-;
R12 = Rm/(sin (ak*l1)*sin (ak*l2))
X12 = Xm/(sin (ak*l1)*sin (ak*l2))
return
end
program rpdip12h
c Radiation patterns of two dipoles (H-plane cut)
c when the wire radius is take into account in the induced EMF'
include 'fgraph.fi'
integer kk
character tx (6)*2,fff (10)*6
dimension asd (6), gr (10,1000)
real la, l1, l2,m, mp
write (*,*)' Strength Simpson integration'
write (*,*)'Sapros lengths, both in mm.'
read (*,*) l1, l2
write (*,*)' l1=', l1,' l2=', l2
write (*,*)'Sapros wire radius (wr) & Lambda (la), both in mm.'
read (*,*) wr, la
write (*,*)' wr=', wr,' la=', la
write (*,*)'Sapros coupling: m<1, ef (radian)'
read (*,*) m, ef
write (*,*)' m=', m,' ef=', ef
write (*,*)'Sapros full distance (d) in mm'
read (*,*) d
write (*,*)' d=', d
write (*,*)'Sapros max. pattern «mp» (first=1, then-value)'
read (*,*) mp
write (*,*) ' mp=', mp
write (*,*)'Sapros of H-plane angles (degrees)'
read (*,*) fin, fih, fib
write (*,*) ' fin=', fin,' fih=', fih,' fib=', fib
write (*,*)'Sapros type driver: 1 — print; 2 — propusk'
read (*,*) isnak
write (*,*) ' isnak=', isnak
pi=4.atan (1.)
ak=2.*pi/la
c ———————— H-plane cut ————————————————;
zh1=l1/1000.
zh2=l2/1000.
fmax=0.05
kk=1
fi=fin
deni11=sin (ak*l1)
deni21=sin (ak*l2)
begini11=1.
endi11=0.0
begini21=1.
endi21=0.0
1 ai11=0.0
ai21=0.0
do 100 i=1,999,2
ai11=ai11+4*sin (ak*(l1-i*zh1))/deni11
100ai21=ai21+4*sin (ak*(l2-i*zh2))/deni21
do 101 i=2,999,2
ai11=ai11+4*sin (ak*(l1-i*zh1))/deni11
101ai21=ai21+4*sin (ak*(l2-i*zh2))/deni21
ai11=zh1*(ai11+begini11+endi11)/3
ai21=zh2*(ai21+begini21+endi21)/3
c ————————— memento: ef is used in radians ————————————;
ah=ai11*cos (ak*d/2*sin (fi*pi/180))+
* m*ai21*cos (ak*d/2*sin (fi*pi/180)+ef)
bh=ai11*cos (ak*d/2*sin (fi*pi/180))+
* m*ai21*cos (ak*d/2*sin (fi*pi/180)+ef)
f=sqrt (ah**2+bh**2)
if (f.gt.fmax) goto 301
goto 302
301fmax=f
302continue
pattern=f/mp
if (isnak .eq. 1) goto 201
goto 202
201write (*,*)' A N G L E (degrees)=', fi
write (*,*)' P A T T E R N =', pattern
write (*,*)'Sapros continue — press enter'
read (*,*)
202continue
gr (1,kk)=pattern
gr (2,kk)=0.1
gr (3,kk)=0.2
gr (4,kk)=0.3
gr (5,kk)=0.4
gr (6,kk)=0.5
gr (7,kk)=0.7071
gr (8,kk)=0.6
gr (9,kk)=0.8
gr (10,kk)=0.999
kk=kk+1
fi=fi+fih
if (fi .le. fib) goto 1
tx (1)= 'l1'
tx (2)= 'l2'
tx (3)= 'd'
tx (4)= 'la'
tx (5)= 'mp'
tx (6)= 'wr'
asd (1)=l1
asd (2)=l2
asd (3)=d
asd (4)=la
asd (5)=mp
asd (6)=wr
fff (1)='H-cut '
fff (2)=' '
fff (3)=' '
fff (4)=' '
fff (5)=' rp '
fff (6)='dip12h'
fff (7)=' '
fff (8)=' '
fff (9)=' '
fff (10)=' '
write (*,*)'maxpattern «mp» =', fmax
write (*,*)'Sapros continue — press enter'
read (*,*)
stop
end
program rpdip12e
cRadiation patterns of two dipoles (E-plane cut)
cwhen the wire radius is take into account in the induced EMF'
include 'fgraph.fi'
integer kk
character tx (6)*2,fff (10)*6
dimension asd (6), gr (10,1000)
real la, l1, l2,m, mp
write (*,*)' Strength simpson integration'
write (*,*)'Sapros lenghs, both in mm'
read (*,*) l1, l2
write (*,*)' l1=', l1,' l2=', l2
write (*,*)'Sapros wire radius & lambda, both in mm'
read (*,*) wr, la
write (*,*) ' wr=', wr,' la=', la
write (*,*)'Sapros coupling: m<1, ef (radian)'
read (*,*) m, ef
write (*,*) ' m=', m,' ef=', ef
write (*,*)'Sap. full dist.(d) mm & «ficut» (degrees)'
read (*,*) d, ficut
write (*,*)' d=', d,' ficut=', ficut
write (*,*)'Sapros max. pattern «mp» (first=1, then-value)'
read (*,*) mp
write (*,*)'mp=', mp
write (*,*)'Sapros of E-plane angles (degrees)'
read (*,*) tetan, tetah, tetab
write (*,*)' tetan=', tetan,' tetah=', tetah,' tetab=', tetab
write (*,*)'Sapros type driver: 1- print; 2 — propust'
read (*,*) isnak
write (*,*)' isnak=', isnak
pi= 4.*atan (1.)
ak = 2.*pi/la
c——————- E-plan cut ————————-;
zh1 = l1/1000.
zh2 = l2/1000.
fmax = 0.05
kk = 1
teta=tetan
deni11 = sin (ak*l1)
deni12 = deni11
deni21 = sin (ak*l2)
deni22=deni21
begini11=1.
begini12=0.
begini21=1.
begini22=0.
endi11=0.0
endi12=0.
endi21=0.
endi22=0.
1t=teta*pi/180.
ai11=0.0
ai12=0.0
ai21=0.0
ai22=0.0
do 100 i=1,999,2
ai11=ai11+4*sin (ak*(l1-i*zh1))*cos (ak*i*zh1*cos (t))/deni11
ai12=ai12+4*sin (ak*(l1-i*zh1))*sin (ak*i*zh1*cos (t))/deni12
ai21=ai21+4*sin (ak*(l2-i*zh2))*cos (ak*i*zh2*cos (t))/deni21
ai22=ai22+4*sin (ak*(l2-i*zh2))*sin (ak*i*zh2*cos (t))/deni22
100continue
do 101 i=2,999,2
ai11=ai11+2*sin (ak*(l1-i*zh1))*cos (ak*i*zh1*cos (t))/deni11
ai12=ai12+2*sin (ak*(l1-i*zh1))*sin (ak*i*zh1*cos (t))/deni11
ai21=ai21+2*sin (ak*(l2-i*zh2))*cos (ak*i*zh2*cos (t))/deni21
ai22=ai22+2*sin (ak*(l2-i*zh2))*sin (ak*i*zh2*cos (t))/deni22
101continue
ai11=zh1*(ai11+begini11+endi11)/3.
ai12=zh1*(ai12+begini12+endi12)/3.
ai21=zh2*(ai21+begini21+endi21)/3.
ai22=zh2*(ai22+begini22+endi22)/3.
c—————- memento: ef is used as radians —————-;
psi=ak*d*sin (t)*sin (ficut*pi-180.)/2.
ae=ai11*cos (psi)+ai12*sin (psi)+m*ai21*cos (psi+ef)
ae=ae-m*ai22*sin (psi+ef)
be=ai12*cos (psi)-ai11*sin (psi)+m*ai21*sin (psi+ef)
be=be+m*ai22*cos (psi+ef)
f=sqrt (ae**2+be**2)*sin (t)
if (f. gt. fmax) goto 301
goto 302
301fmax=f
302 continue
pattern = f/mp
if (isnak .eq. 1) goto 201
goto 202
201 write (*,*)' A N G L E (degrees) =', teta
write (*,*)' P A T T E R N =', pattern
write (*,*)'Sapros continue — press enter'
read (*,*)
202continue
gr (1,kk) = pattern
gr (2,kk)=0.1
gr (3,kk)=0.2
gr (4,kk)=0.3
gr (5,kk)=0.4
gr (6,kk)=0.5
gr (7,kk)=0.7071
gr (8,kk)=0.6
gr (9,kk)=0.8
gr (10,kk)=0.999
kk = kk+1
teta=teta+tetah
if (teta .le. tetab) goto 1
tx (1) = '11'
tx (2) = '12'
tx (3) = 'd'
tx (4) = 'la'
tx (5) = 'mp'
tx (6) = 'fi'
asd (1) = l1
asd (2) = l2
asd (3) = d
asd (4) = la
asd (5) = mp
asd (6) = ficut
fff (1) = 'E-cut'
fff (2) = ' '
fff (3) = ' '
fff (4) = ' '
fff (5) = ' rp '
fff (6) = 'dip12e'
fff (7) = ' '
fff (8) = ' '
fff (9) = ' '
fff (10) = ' '
write (*,*)'maxpattern «mp» =', fmax
write (*,*)'Sapros continue — press enter'
read (*,*)
ccall calc (tetan, tetab, gr, asd, 1000, tx, fff, 10)
stop
end