プログラム

1 parameter (nmax=15)

2 implicit complex (c)

3 complex c(nmax,nmax+1)

4 common akL,akr0,akd

5

6 pi = acos(-1.0d0)

7 c = 3.141592....

8 eta = 120*pi

9 c = sqrt(mu0/epsilon0)

10

11 c step 1

12 aL = 0.25

13 r0 = 0.015625

14 N0 = 3

15 fs = 0.75

16 fe = 1.25

17 df = 0.0625

18 read(*,*) aL,r0,N0,fs,fe,df

19 N0 = 2*(N0/2)+1

20 ncenter = N0/2 + 1

21 write(*,1000) aL,r0,N0,fs,fe,df

22 1000 format('# L, r0 = ',2f10.5/

23 + '# N = ',i5/

24 + '# frequency (starte, end, delta) : ',3f10.5)

25

26 c step 3

27 do 100 f = fs, fe, df

28 c step 3a

29 akL = 2*pi*aL*f

30 akr0 = 2*pi*r0*f

31 akd = 2*akL/(n0+1)*f

32 c step 3b

33 do 10 n=1, n0

34 do 12 m=1,n-1

35 12 c(m,n) = c(n,m)

36 do 10 m=n, n0

37 c(m,n) = (cf(m-n-2)+cf(n-m)+cf(n-m-2)+cf(m-n)

38 + -2*cos(akd)*(cf(m-n-1)+cf(n-m-1)))/(sin(akd))**2

39 10 continue

40 c step 3c

41 do 20 m=1, n0

42 if(m .ne. ncenter) then

43 c(m,n0+1) = 0.0

44 else

45 c(m,n0+1) = cmplx(0.0,4*pi/eta)

46 endif

47 20 continue

48 c step 3d

49 call cswepm(c,nmax,n0,1,1.0e-10, ier)

50 c step 3e 51 cz = 1.0/c(ncenter,n0+1)

52 az = abs(cz)

53 phase = atan2(aimag(cz),real(cz))*180.0/pi

54 write(*,1100) f, az, phase

55 1100 format(f10.5,f20.5,f10.2)

56 100 continue

57 stop

58 end

59

60 complex function cgs(u)

61 implicit complex (c)

62 common akL,akr0,akd

63 common /block1/ n

64

65 z = n*akd+u

66 r = sqrt(z*z+akr0*akr0)

67 cg = cmplx(cos(r)/r,-sin(r)/r)

68 cgs = cg * sin(u)

69 return

70 end

71

72 complex function cf(L)

73 implicit complex (c)

74 common akL,akr0,akd

75 common /block1/ n

76 external cgs

77 n = L

78 call cromb(0.0,akd,cgs,1.0e-6, cs, ier)

79 cf = cs

80 return

81 end

82

83 C -- SIMALTANIOUS LINEAR EQUATIONS

84 C ( SWEEP-OUT METHOD )

85 C

86 subroutine cswepm(c,nmax,nsize,nright,eps,ier)

87 implicit complex (c)

88 complex c(nmax,nmax)

89 INTEGER PIVOT

90 if(eps.le.0.0) eps=1.0e-6

91 ier = 0

92 C

93 C STEP 2

94 DO 100 I=1,nsize

95 C STEP 3.1

96 PIVOT=I

97 cMAX=c(I,I)

98 amax=abs(cmax)

99 DO 30 J=I+1,nsize

100 WORK=ABS(c(J,I)) 101 IF( WORK .GT. AMAX ) THEN

102 PIVOT=J

103 cmax=c(j,i)

104 AMAX=WORK

105 ENDIF

106 30 CONTINUE

107 C STEP 3.2

108 IF( PIVOT .NE. I ) THEN

109 C EXCHANGE: c(I,*), c(PIVOT,*)

110 DO 32 J=I,nsize+nright

111 cWORK =c(I,J)

112 c(I,J) =c(PIVOT,J)

113 c(PIVOT,J)=cWORK

114 32 CONTINUE

115 ENDIF

116 C STEP 4

117 c if((amax .lt. eps) .and. (i .gt. 1)) then

118 if( amax .lt. eps ) then

119 ier = -1000

120 write(*,*) 'CSWEPM: ERROR! ier & PIVOT=', ier, AMAX

121 return

122 end if

123 cWORK=1.0/c(I,I)

124 DO 40 J=I+1,nsize+nright

125 40 c(I,J)=c(I,J)*cWORK

126 C STEP 5

127 DO 50 J=I+1,nsize

128 cWORK=-c(J,I)

129 DO 50 K=I+1,nsize+nright

130 c(J,K)=c(J,K)+c(I,K)*cWORK

131 50 CONTINUE

132 C

133 100 CONTINUE

134 C

135 C STEP 7

136 DO 70 I=nsize,2,-1

137 DO 70 J=1,I-1

138 do 70 k=1,nright

139 c(J,nsize+k)=c(J,nsize+k)-c(I,nsize+k)*c(J,I)

140 70 CONTINUE

141 return

142 END

143

144 C -- NUMERICAL INTEGRATION BY ROMBERG'S RULE

145 C

146 C ROMBERG'S INTEGRAL ( SUB PROGRAM )

147 C

148 SUBROUTINE cROMB(LOWER,UPPER,cFUNC,EPS,cROM,IER)

149 C INTEGRATION BY THE ROMBERG'S RULE.

150 C [LOWER,UPPER] RANGE OF INTEGRAL 151 C FUNC INTEGRAND

152 C EPS ERROR CRITERION

153 C ROMBRG RESULT OF THE INTEGRAL

154 C IER = 0 LOWER=UPPER

155 C K NO. OF DIVISION: 2**K (K<=20)

156 C -K NOT ACCURATE (WITHIN EPS)

157 implicit complex (c)

158 PARAMETER (MAXMUM=20)

159 REAL LOWER

160 complex cwork(20)

161 DELTA=UPPER-LOWER

162 IF(DELTA.EQ.0.0) THEN

163 IER=0

164 crom=0.0

165 RETURN

166 END IF

167 IF(EPS.LE.0.0) EPS=1.0E-6

168 C STEP 2

169 csum=0.5*(cfunc(LOWER)+cfunc(UPPER))

170 cs=csum

171 ND=1

172 C STEP 3

173 DO 300 K=1,MAXMUM

174 ND=ND*2

175 C STEP 3.1

176 DO 310 J=1,ND,2

177 csum=csum+cfunc(LOWER+(DELTA/ND)*J)

178 310 CONTINUE

179 cwork(K)=cs

180 cs=csum/ND

181 C STEP 3.2

182 L=1

183 DO 320 N=1,K

184 L=L*4

185 cw=cs

186 cs=(L*cs-cwork(N))/(L-1)

187 cwork(N)=cw

188 320 CONTINUE

189 C STEP 3.3

190 IF(K.GT.3) THEN

191 IF(ABS((cwork(K)-cs)).LE.ABS(EPS*cs)) THEN

192 C STEP 4

193 IER=K

194 crom=cs*DELTA

195 RETURN

196 END IF

197 IF(abs(cwork(K)).EQ.abs(cwork(K-1))) GO TO 900

198 END IF

199 300 CONTINUE

200 900 IER=-K 201 crom=cs*DELTA

202 RETURN

203 END

• 複素数型の変数はすべてcで始まる名前を用いることとし、 各プログラム単位の先頭で
  implicit complex (c)
  を宣言している。
• 分割数、したがって、展開関数は最大１５を想定して、 第１行でnmax=15としている。
• ４行目のcommon文は$F_\ell$の計算においてパラメータを 関数に渡すために使用する。
• 41行からの部分で給電部の電圧は電流の展開関数の内、 中心の部分を除いて０であるので、それをif文で判定している。
• 60行からのfunction cgs(u)は $G(\ell\Delta+u)\sin u$を計算する関数である。 u以外のパラメータはcommon文によりメインルーチンからと 関数cfより受け取っている。

  $G(z)$の計算には
  

  $$G(z) = \frac{e^{-jr}}{r} = \frac{\cos r}{r} - j \frac{\sin r}{r}$$
  
  
  を利用している。( $r=sqrt{z^2+r_0^2}$)

  72行からのfunction cf(k)は
  

  $$F_k = \int_0^\Delta G(k\Delta+u)\sin u du$$
  
  
  を計算している。 ロンベルグ積分(crombにより数値積分を行なっている。 k以外のパラメータはcommon文を用いてメインルーチンより受け取っている。

  83行からはガウスの消去法により連立方程式を解くサブルーチンである。 また、 144行以降はロンベルグ積分による数値積分のためのサブルーチンである。