Оптимизация двухдипольной излучающей системы

Министерство образования и науки Российской Федерации НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ Кафедра РПиРПУ Пояснительная записка к курсовой работе по дисциплине «Устройства СВЧ и антенны»

Новосибирск, 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