Ứng dụng phương pháp moment trong bài toán phân tích các kết cấu điện từ phẳng được kích thích bởi sóng chạy
Trích yếu luận án - Tên luận án: ứng dụng phương pháp moment trong bài toán phân tích các kết cấu điện từ phẳng được kích thích bởi sóng chạy. - Ngành khoa học của luận án: Thông tin vô tuyến, phát thanh và vô tuyến truyền hình. Mã số chuyên ngành: 2.07.02 - Tên cơ sở đào tạo: Trường Đại học Bách Khoa Hà Nội. a) Đối tượng nghiên cứu của luận án: Trong những thập kỷ 80 - 90 của thế kỷ XX, thế giới đã được chứng kiến những ứng dụng của vi mạch tích hợp trong các thiết bị điện tử, thông tin liên lạc phục vụ an ninh quốc phòng và đời sống hàng ngày. Việc sử dụng các vi mạch tích hợp (kết cấu mạch dải và khe dải là các thành phần cơ bản) có −u điểm dễ dàng và linh hoạt trong thiết kế mạch và nâng cao tính khai thác của kết cấu. Một ứng dụng quan trọng trong lĩnh vực siêu cao tần đó là các kết cấu truyền dẫn sóng chu kỳ (còn gọi là "kết cấu chu kỳ"). Sự quan tâm đến các kết cấu dẫn sóng loại này nhờ hai tính chất cơ bản của chúng là: (i) các đặc tính lọc thông băng và chặn băng tần; (ii) hỗ trợ các sóng có vận tốc pha nhỏ hơn vận tốc ánh sáng (sóng chậm). Luận án này đi sâu vào hướng nghiên cứu tổng hợp và phân tích tính chất thứ hai của kết cấu chu kỳ đó là tính chất hỗ trợ các sóng chậm sử dụng các kết cấu mạch dải phẳng và kết cấu sóng rò phẳng được kích thích bởi sóng chạy. b) Mục đích nghiên cứu: - Trên thực tế để tạo ra các đồ thị phương hướng (sóng thứ cấp) theo yêu cầu, bề mặt kết cấu thường có dạng hết sức phức tạp. Do vậy việc phân tích các kết cấu này gặp rất nhiều khó khăn đặc biệt phải tính toán đối với các phương trình đường cong hình học rất phức tạp. Nhiều nhà khoa học như Aizenberg G. Z.; Yampolski V. G.; Cheriosin O. N.; Tereshin O. N.; Sedov V. M. và Chaplin A. F. trong các nghiên cứu của mình cũng đã rất cố gắng để giải quyết bài toán tổng hợp để tìm ra mô hình đường cong của kết cấu có hình dạng bất kỳ được kích thích bởi sóng chạy. Tuy nhiên không phải là đối với bài toán nào cũng ra được nghiệm vì sử dụng phương pháp tính nghiệm là phương pháp bình phương nhỏ nhất chỉ cho phép tính toán đối với các phép toán giải tích và nhiều khi phương trình tích phân lại có dạng không khả tích. - Các phương pháp số như phương pháp phần tử hữu hạn (Finite element method), phương pháp sai phân hữu hạn (Finite difference method) chưa phát huy được hiệu quả. Luận án đã giải quyết bài toán tổng hợp, phân tích và mô phỏng các kết cấu có hình dạng bất kỳ được kích thích bởi sóng chạy thành các kết cấu phẳng (kết cấu mạch dải và sóng rò) sử dụng phương pháp số cho phép nhận được kết quả chính xác với thời gian ngắn. c) Các kết quả chính và kết luận: Luận án đã giải quyết được một số điểm đột phá như sau: - Thực hiện bài toán tổng hợp nhằm đưa một kết cấu có hình dạng bất kỳ có trở kháng bề mặt là đại lượng ảo chuyển thành một kết cấu phẳng có trở kháng bề mặt là đại lượng phức bảo đảm được mọi tính chất điện từ trường của kết cấu ban đầu. - Thực hiện bài toán phân tích kết cấu phẳng có trở kháng bề mặt là đại lượng phức để đánh giá kết quả khi chuyển kết cấu có hình dạng bất kỳ thành kết cấu phẳng. - Sử dụng phương pháp moment (MoM) với hàm cơ sở miền con để phân tích kết cấu. Đây là phương pháp tính toán sử dụng lý thuyết rời rạc để làm giảm nhẹ đáng kể bài toán về mối tương quan của các đại lượng vật lý trong môi trường tự do được biểu diễn qua các phương trình Maxwell và các điều kiện bờ, để biến đổi thành các phương trình tích phân có miền đ−ợc giới hạn và đủ nhỏ. Kích thước nhỏ của miền là vô cùng quan trọng vì phù hợp với kích cỡ RAM của máy tính luôn không phải là một nguồn tài nguyên dồi dào. Đây chính là −u điểm của MoM so với các phương pháp số khác. Đặc biệt MoM rất thuận tiện khi khảo sát các kết cấu phẳng. Những kết quả này cho phép mở rộng phạm vi ứng dụng của bài toán tới phạm vi rộng rãi hơn. - Nghiên cứu hai dạng bài toán đặc biệt chưa được nghiên cứu trong thực tế đó là: + Kết cấu có dạng như kết cấu sóng rò được kích thích bởi sóng chạy dưới góc tới θi bất kỳ trên bề mặt kết cấu. + Kết cấu có dạng như kết cấu sóng mặt (kiểu kết cấu mạch dải) được kích thích liên tục bởi sóng chạy dưới góc tới bất kỳ. - Các chương trình Matlab và Fortran được sử dụng để thực hiện bài toán mô phỏng bằng MoM. Thời gian mô phỏng trên máy tính nhanh hơn so với các kết quả nghiên cứu trước kia. - Với những kết quả đã đạt được, có thể nhận thấy rằng khả năng mô phỏng bằng phương pháp số đối với kết cấu mạch dải và sóng rò là khá chính xác. d) ứng dụng của hai dạng bài toán và kết cấu nghiên cứu - Giảm nhẹ kích thư−ớc các cấu tử nhờ áp dụng những kết cấu mới nh− kết cấu mạch dải và sóng rò một cách phù hợp. - Dễ dàng được sản xuất với chi phí thấp nhờ công nghệ cấy hàng ngàn các cấu tử siêu cao tần sóng được đưa vào cùng một quá trình. - Các kết cấu nghiên cứu rất mỏng và nhẹ. Việc gắn chúng lên thân các thiết bị không gây ảnh hưởng đến bề mặt của thiết bị. - Kết hợp các kết cấu sóng chậm này với các phần tử hay mạch tích cực để có anten tích cực. Mục lục Lời cam đoan .2 mục lục 3 danh mục các ký hiệu, các chữ viết tắt .5 danh mục các hình vẽ 7 mở đầu 9 chương 1: kết cấu điện từ được kích thích bởi sóng chạy 12 1.1. Giới thiệu về các kết cấu được kích thích bởi sóng chạy 12 1.1.1. Kết cấu sóng rò 12 1.1.2. Kết cấu sóng mặt 17 1.1.3. Các quan điểm phân tích kết cấu điện từ được kích thích bởi sóng chạy: .20 1.1.4. Những hạn chế trong bài toán phân tích các kết cấu được kích thích bởi sóng chạy và phương hướng giải quyết .24 1.2. Bài toán tổng hợp kết cấu sóng chạy (kết cấu impedance) 26 1.2.1. Xác định hàm số mặt cong của bề mặt kết cấu impedance và phân bố trở kháng bề mặt .26 1.2.2. Xây dựng mô hình mô phỏng kết cấu impedance có hình dạng bất kỳ .28 1.3. Bài toán phân tích kết cấu sóng chạy (kết cấu impedance) có hình dạng mặt cắt (Profile) bất kỳ 32 1.3.1. Phương trình tích phân đối với các bề mặt trở kháng có mặt cắt biến đổi ít. .32 1.3.2. Bài toán phân tích .34 1.3.3. Đánh giá sai số của phương pháp tổng hợp 37 1.4. Xây dựng kết cấu phẳng được kích thích bởi sóng chạy sử dụng kết cấu mạch dải và kết cấu khe trên hốc cộng hưởng .41 1.4.1. Đặt vấn đề 41 1.4.2. Tính chất điện từ của cấu trúc răng lược và cấu trúc gấp khúc 42 1.4.3. Các kết cấu được nghiên cứu .45 1.5. Kết luận 46 Chương 2: phân tích kết cấu sóng rò được kích thích bởi sóng chạy bằng phương pháp moment .48 2.1. Phương trình tích phân cho kết cấu khe có hình dạng bất kỳ trên hốc cộng h−ởng được kích thích bởi sóng chạy .48 2.1.1. Xác định phương trình điều kiện biên 48 2.1.2. Xác định trường bức xạ trong miền I .49 2.1.3. Xác định trường bức xạ trong miền II 51 2.2. Giải quyết bài toán bằng phương pháp moment 52 2.2.1. Nghiên cứu cấu trúc .52 2.2.2. Chọn hàm cơ sở và thiết lập phương trình ma trận .52 2.2.3. Xác định trường bức xạ 57 2.3. Kết quả mô phỏng .59 2.4. Kết luận 67 Chương 3: phân tích kết cấu sóng mặt (kết cấu mạch dải) kích thích bởi sóng chạy bằng phương pháp moment .68 3.1. Giới thiệu kết cấu mạch dải 68 3.2. Bài toán tổng quát phân tích kết cấu mạch dải có hình dạng bất kỳ sử dụng phương pháp moment 70 3.2.1. Xác định phương trình điều kiện biên và các thành phần của hàm Green 70 3.2.2. Xác định sự phân bố dòng trên bề mặt cấu trúc 71 3.2.3. Xác định phương trình ma trận và ma trận trở kháng 73 3.2.4. Xác định trường tán xạ và mặt cắt tán xạ ngược .74 3.2.5. Các kết quả mô phỏng .75 3.3. Phân tích kết cấu mạch dải hẹp hình dạng bất kỳ đ−ợc kích thích bởi sóng chạy bằng phương pháp moment .79 3.3.1. Những căn cứ xây dựng kết cấu mạch dải hẹp có hình dạng bất kỳ 79 3.3.2. Xác định phương trình điều kiện biên 79 3.3.3. Xác định sự phân bố dòng trên bề mặt kết cấu 80 3.3.4. Chọn hàm cơ sở và xác định phương trình ma trận 81 3.3.5. Xác định ma trận trở kháng 83 3.3.6. Xác định các tích phân Sommerfeld 87 3.3.7. Các kết quả mô phỏng 92 3.4. Kết luận 97 chương 4: kết luận .98 4.1. Nhận xét các kết quả đạt được .98 4.2. ứng dụng của kết cấu điện từ được kích thích bởi sóng chạy .99 4.3. Hướng nghiên cứu trong tương lai .101 danh mục công trình của tác giả .102 tài liệu tham khảo .103 phụ lục 1: giới thiệu phương pháp moment 105 Phụ lục 2: hàm số biểu diễn mặt cong z0(y) của kết cấu 115 Phụ lục 3: phân bố trở kháng trên bề mặt của kết cấu 116 Phụ lục 4: Dạng hình học của kết cấu được nghiên cứu .117 Phụ lục 5: Chương trình Matlab tính toán cấu trúc sóng rò kiểu khe hẹp có hình dạng bất kỳ trên hốc cộng hưởng được kích thích bởi sóng chạy .119 Phụ lục 6: phân tích hàm green, mặt cắt bức xạ và Hiệu ứng biên của kết cấu mạch dải .126 Phụ lục 7: xác định tích phân Sommerfeld đoạn cuối 134 Phụ lục 8: Chương trình fortran tính toán kết cấu mạch dải tổng quát được kích thích bởi sóng chạy .136 Phụ lục 9: Chương trình fortran tính toán kết cấu mạch dải hẹp hình dạng bất kỳ kích thích bởi sóng chạy .150
Bạn đang xem trước 20 trang tài liệu Ứng dụng phương pháp moment trong bài toán phân tích các kết cấu điện từ phẳng được kích thích bởi sóng chạy, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
j*K2*d)*cos(Theta)*Ptx
22
Ety= (377.0/(2.0*pi))*cexp(cj*K2*d)*cos(Theta)*Pty
Epx= (377.0/(2.0*pi))*cexp(cj*K2*d)*cos(Theta)*Ppx
Epy= (377.0/(2.0*pi))*cexp(cj*K2*d)*cos(Theta)*Ppy
ccccc
147
F1=dx*dy*(sin(Ky*dy2)/(Ky*dy2))*((sin(Kx*dx2)/(Kx*dx2))**2)
F2=dx*dy*(sin(Kx*dx2)/(Kx*dx2))*((sin(Ky*dy2)/(Ky*dy2))**2)
do 455 ip=1,mm
do 454 iq=1,nnp1
irow= ((ip-1)*nnp1) + iq
Ph= cexp(cj*((-Kx*x(ip))+(-Ky*y(iq))+(Ky*dy2)))
V(irow)= Di * ((Etx*Et)+(Epx*Ep)) *F1*Ph
if ((ip.eq.6).and.(iq.eq.7)) Vmid=V(irow)
ccccc write(6,470)ip,iq,V(irow)
454 continue
455 continue
do 465 ip=1,mmp1
do 460 iq=1,nn
irow= ((ip-1)*nn) + iq + nmax1
Ph= cexp(cj*((-Kx*x(ip))+(-Ky*y(iq))+(Kx*dx2)))
V(irow)= Di * ((Ety*Et)+(Epy*Ep)) *F2*Ph
ccccc if ((ip.eq.6).and.(iq.eq.6)) Vmid=V(irow)
ccccc write(6,470)ip,iq,V(irow)
460 continue
465 continue
470 format(i3,i3,f12.8,f12.8)
CCCCC
CCCCC Now define shape function
CCCCC
500 Sx=1
Sy=1
do 550 m=1,mm
do 540 n=1,nnp1
dist=sqrt((x(m)**2)+((y(n)-dy2)**2))
if (dist.gt.0.023) Sx(m,n)=0
540 continue
550 continue
do 580 m=1,mmp1
do 570 n=1,nn
dist=sqrt(((x(m)-dx2)**2)+(y(n)**2))
if (dist.gt.0.023) Sy(m,n)=0
570 continue
580 continue
if (IF.eq.2000) write(6,605)((Sx(m,n),m=1,mm),n=nnp1,1,-1)
605 format(12I1)
if (IF.eq.2000) write(6,606)((Sy(m,n),m=1,mmp1),n=nn,1,-1)
606 format(13I1)
irownew=0
do 650 ip=1,mm
do 640 iq=1,nnp1
icolnew=0
if (Sx(ip,iq).eq.0) goto 640
23
irow= ((ip-1)*nnp1) + iq
irownew=irownew+1
VN(irownew)=V(irow)
do 620 m=1,mm
do 610 n=1,nnp1
if (Sx(m,n).eq.0) goto 610
icol= ((m-1)*nnp1) + n
icolnew=icolnew+1
ZZN(irownew,icolnew)= ZZ(irow,icol)
148
610 continue
620 continue
do 630 m=1,mmp1
do 625 n=1,nn
if (Sy(m,n).eq.0) goto 625
icol= ((m-1)*nn) + n + nmax1
icolnew=icolnew+1
ZZN(irownew,icolnew)= ZZ(irow,icol)
625 continue
630 continue
640 continue
650 continue
do 700 ip=1,mmp1
do 690 iq=1,nn
icolnew=0
if (Sy(ip,iq).eq.0) goto 690
irow= ((ip-1)*nn) + iq + nmax1
irownew=irownew+1
VN(irownew)= V(irow)
do 660 m=1,mm
do 655 n=1,nnp1
if (Sx(m,n).eq.0) goto 655
icol= ((m-1)*nnp1) + n
icolnew=icolnew+1
ZZN(irownew,icolnew)= ZZ(irow,icol)
655 continue
660 continue
do 680 m=1,mmp1
do 670 n=1,nn
if (Sy(m,n).eq.0) goto 670
icol= ((m-1)*nn) + n + nmax1
icolnew= icolnew+1
ZZN(irownew,icolnew)= ZZ(irow,icol)
670 continue
680 continue
690 continue
700 continue
J=(0.0,0.0)
750 call l2acg(irownew,ZZN,nmax,VN,1,J,fac,ipvt,wk)
CCCCC
CCCCC Now find RCS at angle (Theta,Phi)
CCCCC
24
irow=0
Et=(0.0,0.0)
Ep=(0.0,0.0)
do 1100 ip=1,mm
do 1050 iq=1,nnp1
if (Sx(ip,iq).eq.0) goto 1050
irow= irow+1
ccccc write(6,1040)float(ip),float(iq),cabs(J(irow))
1040 format(f8.2,f8.2,f14.5)
if ((ip.eq.6).and.(iq.eq.7)) Jmid=J(irow)
Ph= cexp(cj*((-Kx*x(ip))+(-Ky*y(iq))+(Ky*dy2)))
Et= Et+(J(irow)*Etx*F1*Ph)
Ep= Ep+(J(irow)*Epx*F1*Ph)
1050 continue
149
1100 continue
do 1200 ip=1,mmp1
do 1150 iq=1,nn
if (Sy(ip,iq).eq.0) goto 1150
irow= irow+1
ccccc write(6,1040)float(ip),float(iq),cabs(J(irow))
ccccc if ((ip.eq.6).and.(iq.eq.6)) Jmid=J(irow)
Ph= cexp(cj*((-Kx*x(ip))+(-Ky*y(iq))+(Kx*dx2)))
Et= Et+(J(irow)*Ety*F2*Ph)
Ep= Ep+(J(irow)*Epy*F2*Ph)
1150 continue
1200 continue
1300 format(i3,i3,f12.8,f12.8)
RCS= 4.0*pi*cabs(Et)*cabs(Et)
ccccc RCS= 4.0*pi*cabs(Ep)*cabs(Ep)
RCS= 10.0*alog10(RCS)
JN=Jmid/Vmid
write(6,1500)(fr/(1.0E+09)),RCS,JN
1400 format('Fr=',f6.2,' RCS = ',f12.6,' dBsm')
1500 format(f8.4,f12.6,2x,e12.6,2x,e12.6)
ccccc write(6,1600)Vmid,Jmid
1600 format(4e14.6)
2000 continue
stop
end
ccccc
150
Phô lôc 9: Ch−¬ng tr×nh fortran tÝnh to¸n kÕt cÊu
m¹ch d¶i hÑp h×nh d¹ng bÊt kú kÝch thÝch bëi sãng ch¹y
PROGRAM CALCULATION_Zmn
C
EXTERNAL FUNSUB
INTRINSIC CDEXP,CMPLX,MOD
INTEGER NUMFUN,NW,MINPTS,MAXPTS,KEY,RESTAR,FUNCT
INTEGER WRKSUB,EMAX,NUMNUL,NUM
PARAMETER (NUM=21,NUMNUL=16)
PARAMETER (NUMFUN=1,WRKSUB=100,EMAX=25)
C
PARAMETER (NW=1+5*NUMFUN+4*WRKSUB*NUMFUN+3*WRKSUB+(EMAX+1)**2+
+(NUMNUL+1+NUM)*NUMFUN)
C
PARAMETER (KEY=2)
C
INTEGER(2) tmphour, tmpminute, tmpsecond, tmphund
INTEGER(2) tmphour1, tmpminute1, tmpsecond1, tmphund1
C paramenter for
PARAMETER (NM=271)
INTEGER PP,QQ,JJ
DOUBLE PRECISION Xn(NM),Yn(NM),Xm(NM),Ym(NM),l0,RA,AA,ZZ,DELT
DOUBLE PRECISION Znm(NM,NM),InvZnm(NM,NM),PIS(NM,NM),PIO(NM,NM)
DOUBLE PRECISION Iunit(NM,NM),SS,Vc,Vmr(NM),Vmi(NM),Sm(NM),theta,d
DOUBLE PRECISION EPXI0,EPXIR,OMEGA,Fi,BETAY,K0,KY,F11,F12,F21,F22
DOUBLE PRECISION Imr(NM),Imi(NM),Muy0
C
DOUBLE PRECISION FTT,PTT,CTH,STH,THETAT,RAD,ETA,ETMM(181)
DOUBLE COMPLEX CRT,Ji,Im(NM)
C
DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,GAMMA,PERIOD,
+RESULT(NUMFUN),ABSERR(NUMFUN),WORK(NW),PI
INTEGER NEVAL,IFAIL,IWORK(3*WRKSUB)
LOGICAL DONE(NUMFUN)
C Caculate the initial time of program
CALL GETTIM(tmphour, tmpminute, tmpsecond, tmphund)
tmphour1 = tmphour
tmpminute1 = tmpminute
tmpsecond1=tmpsecond
tmphund1=tmphund
PI=2*ASIN(1.D0)
C Parameter for calculate integration!
PERIOD=2*PI
GAMMA=1.0D0
A(1)=0.
B(1)=3
RESTAR=0
MINPTS=42
MAXPTS=550
EPSABS=0
EPSREL=1.D-9
C Parameter for calculate main program!
PP = 1
QQ = 2
151
l0 = 0.003
AA = 0.25E-3
d = 0.01
DELT = 1.0E-5
F11 = 0.
F12 = 0.
F21 = 0.
F22 = 0.
EPXI0 = 8.85E-12
EPXIR = 2.2
Vc = 3.0E+8
Fi = 2.5E+9
OMEGA = 2*PI*Fi
Muy0 =1.26E-6
K0 = 2*PI*Fi/Vc
C KY=K0/BETAY
BETAY = 1.*INT((NM-1)/(PP+QQ))/(NM-1)
KY = K0/BETAY
C
WRITE (*,*) 'A=',A(1),' B=',B(1)
WRITE (*,*) 'GAMMA=',GAMMA,' PERIOD=',PERIOD
WRITE (*,*) 'MAXPTS=',MAXPTS,' WRKSUB=',WRKSUB,' NW=',NW
WRITE (*,*) 'KEY=',KEY,' EPSREL=',EPSREL
C SET Znm =0.
DO I=1,NM
DO J=1,NM
Znm(I,J)=0.0
ENDDO
ENDDO
C
C CALCULATE Znm
DO I=1,NM
DO J=1,NM
CALL COORDINATE(NM,Xn,Yn,Xm,Ym,PP,QQ,l0)
RA = RS(NM,Xn,Yn,Xm,Ym,AA,I,J)
C
C CALCULATE FUNCT=1
C
FUNCT=1
ZZ=0.01
CALL DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS,
+MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,DONE,
+WORK,IWORK,FUNCT,RA,ZZ)
F11=RESULT(1)
ZZ=0.01-DELT
CALL DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS,
+MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,DONE,
+WORK,IWORK,FUNCT,RA,ZZ)
F12=RESULT(1)
PIS(I,J)=(F11-F12)/DELT/2/PI/EPXI0/OMEGA
C CALCULATE FUNCT=2
FUNCT=2
ZZ=0.01
CALL DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS,
+MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,DONE,
+WORK,IWORK,FUNCT,RA,ZZ)
F21=RESULT(1)
152
ZZ=0.01-DELT
CALL DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS,
+MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,DONE,
+WORK,IWORK,FUNCT,RA,ZZ)
F22=RESULT(1)
PIO(I,J)=(1-EPXIR)*(F21-F22)/DELT/2/PI/EPXI0/OMEGA
C WRITE(*,*)'RA(',I,J,'=',PIS(I,J),PIO(I,J)
ENDDO
ENDDO
C
C CALCULATE Znm
C
DO I=1,NM
DO J=1,1
Znm(I,J)=PIS(I,J+1)-2*COS(KY*l0)*PIS(I,J)
+ +(KY+1)*(2*COS(KY*l0)*PIO(I,J)-PIO(I,J+1))
Znm(I,J)=Znm(I,J)*KY/SIN(KY*l0)/(KY**2+1)
ENDDO
ENDDO
DO I=1,NM
DO J=2,NM-1
Znm(I,J)=PIS(I,J+1)+PIS(I,J-1)-2*COS(KY*l0)*PIS(I,J)
+ +(KY+1)*(2*COS(KY*l0)*PIO(I,J)-PIO(I,J+1)-PIO(I,J-1))
Znm(I,J)=Znm(I,J)*KY/SIN(KY*l0)/(KY**2+1)
C WRITE(*,*)'Znm(',I,J,'=',Znm(I,J)
ENDDO
ENDDO
DO I=1,NM
DO J=NM,NM
Znm(I,J)=PIS(I,J-1)-2*COS(KY*l0)*PIS(I,J)
+ +(KY+1)*(2*COS(KY*l0)*PIO(I,J)-PIO(I,J-1))
Znm(I,J)=Znm(I,J)*KY/SIN(KY*l0)/(KY**2+1)
ENDDO
ENDDO
C
C INVERSE MATRIX Znm -> InvZnm
C
NRA = NM
NCA = NM
LDA = NM
LDGINV = NM
TOL = 10.*AMACH(4)
CALL DLSGRR (NRA, NCA, Znm, LDA, TOL, IRANK, InvZnm, LDGINV)
DO I=1,NM
DO J=1,NM
WRITE(*,*)InvZnm(I,J)
ENDDO
ENDDO
C CHECK Znm.InvZnm=Iunitnm
DO I=1,NM
DO J=1,NM
SS=0.0
DO JJ=1,NM
SS=SS+ Znm(JJ,I)*InvZnm(J,JJ)
ENDDO
IF (ABS(SS) .LT. 3.E-4) THEN
SS=0.0
153
ENDIF
Iunit(I,J)=SS
ENDDO
ENDDO
C Calculate Vm
theta = PI/2.
DO I=1,NM
Sm(I)=DSQRT(Xm(I)**2+Ym(I)**2)
ENDDO
DO I=1,NM
Vmi(I)= -(4*PI/OMEGA/Muy0)*(COS(KY*d)*COS(KY*SIN(theta)*Sm(I))
& - SIN(KY*d)*SIN(KY*SIN(theta)*Sm(I)))
Vmr(I)= -(4*PI/OMEGA/Muy0)*(SIN(KY*d)*COS(KY*SIN(theta)*Sm(I))
& +COS(KY*d)*SIN(KY*SIN(theta)*Sm(I)))
ENDDO
C Calculate Inm
DO I=1,NM
SS=0.0
DO J=1,NM
SS=SS+InvZnm(I,J)*Vmi(J)
ENDDO
Imi(I)=SS
Imr(I)=Vmr(I)
ENDDO
C SAVE Znm ON FILE Znm.TXT
C
OPEN(10, FILE='Inm.TXT')
REWIND(10)
DO I=1,NM
WRITE(10,*)I,' ',Imr(I),' ',Imi(I)
ENDDO
ENDFILE(10)
C MAIN CALCULATE
Ji=CMPLX(0.,1.)
OPEN(11, FILE='Im.TXT')
REWIND(11)
DO I=1,NM
Im(I)=CMPLX(Imr(I),Imi(I))
WRITE(11,*)I,' ',Im(I)
ENDDO
ENDFILE(11)
RAD = PI/180
ETA = 120*PI
DO I=1,181
THETAT = (I-1)*RAD
CTH = COS(THETAT)
STH = SIN(THETAT)
IF (DABS(CTH) .LE. 1.E-3) THEN
FTT=1.
ELSE
FTT=SIN(KY*l0*CTH*0.5)/(KY*l0*CTH*0.5)
ENDIF
CRT=0.
DO J=1,NM
CRT=CRT+CDEXP(Ji*KY*Sm(J)*CTH)*FTT*Im(J)*l0
ENDDO
PTT=CDABS(CRT)*STH*STH*ETA*0.5
154
ETMM(I)=PTT
ENDDO
DO I=1,181
PTT=ETMM(I)**2*KY*2
IF (PTT .LE. 1.E-10) PTT=1.E-10
ETMM(I)=10.*DLOG10(PTT)
ENDDO
C SAVE ETMM ON FILE ETMM.TXT!
OPEN(20, FILE='ETMM.TXT')
REWIND(20)
DO I=1,181
WRITE(20,*)I,' ',ETMM(I)
ENDDO
ENDFILE(20)
C TEST!
FUNCT=3
CALL DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS,
+MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,IFAIL,DONE,
+WORK,IWORK,FUNCT,RA,ZZ)
WRITE (*,*)
WRITE (*,*) 'Result:',FUNCT
WRITE (*,*) RESULT
WRITE (*,*) ABSERR
C WRITE (*,*) 'NEVAL=',NEVAL,' IFAIL=',IFAIL
CALL GETTIM(tmphour, tmpminute, tmpsecond, tmphund)
WRITE(*,*)'THE TIME LASTED:',(tmphour-tmphour1),'h',INT((tmpminute
& *60+tmpsecond-tmpminute1*60-tmpsecond1)/60),'ph',MOD(
& 1.*(tmpminute*60+tmpsecond-tmpminute1*60-tmpsecond1),60.),'s'
END
C ***************************************
SUBROUTINE FUNSUB (X,NUMFUN,NP,FUNVLS,FUNCT,RA,ZZ)
INTEGER NUMFUN,I,NP,J,K,FUNCT
DOUBLE PRECISION X(NP),GAMA0(NP),GAMA1(NP),FUNVLS(NP,1),Y,RA,ZZ
DOUBLE PRECISION DEXP,BSJ0,k0,epsilonR,d,delta1,PI,Vc,Fi
+ r,a
INTRINSIC DEXP,DCOSH,DSINH,DTANH,DSQRT,DABS
EXTERNAL DBSJ0,CONST
PI = 3.1415
Vc=3.0E+8
Fi=2.5E+9
k0=2*PI*Fi/Vc
epsilonR=2.2
d=0.01
C delta1=1.E-5
delta1=ZZ-d
C a=0.0025
a=RA
DO I=1,NP
Y=0.0
DO J=1,NUMFUN
FUNVLS(I,J)=Y
ENDDO
ENDDO
C CHOOSE FUNCTION PI_S (FUNCT=1)
155
IF (FUNCT .EQ. 1) THEN
DO I=1,NP
GAMA0(I)=DSQRT(DABS(X(I)*X(I)-k0*k0))
GAMA1(I)=DSQRT(DABS(X(I)*X(I)-epsilonR*k0*k0))
Y=DBSJ0(a*X(I))*DEXP(-GAMA0(I)*delta1)*X(I)/
+ (GAMA0(I)+GAMA1(I)*COSH(d*GAMA1(I))/SINH(d*GAMA1(I)))
DO J=1,NUMFUN
FUNVLS(I,J)=Y
ENDDO
ENDDO
C CHOOSE FUNCTION PIO (FUNCT=2)
ELSE IF (FUNCT .EQ. 2) THEN
DO I=1,NP
GAMA0(I)=DSQRT(DABS(X(I)*X(I)-k0*k0))
GAMA1(I)=DSQRT(DABS(X(I)*X(I)-epsilonR*k0*k0))
Y=DBSJ0(a*X(I))*GAMA0(I)*DEXP(-GAMA0(I)*delta1)*X(I)/
+ (GAMA0(I)+GAMA1(I)*COSH(d*GAMA1(I))/SINH(d*GAMA1(I)))/
+ (epsilonR*GAMA0(I)+GAMA1(I)*TANH(d*GAMA1(I)))
DO J=1,NUMFUN
FUNVLS(I,J)=Y
ENDDO
ENDDO
END IF
IF (FUNCT .EQ. 3) THEN
DO I=1,NP
Y=(2*SIN(X(I))+X(I)*COS(X(I)))
DO J=1,NUMFUN
FUNVLS(I,J)=FUNCT*Y/(J+X(I)**2)
ENDDO
ENDDO
ENDIF
RETURN
END
C
C ************************************************
C Subroutine to determine the coordinate of the points.
SUBROUTINE COORDINATE(NM,Xn,Yn,Xm,Ym,P,Q,l0)
INTEGER I,J,P,Q,N,M,NODE
DOUBLE PRECISION Xnm(NM),Ynm(NM),Xn(NM),Yn(NM),Xm(NM),Ym(NM),l0
NODE=0
DO I=1,NM
NODE = INT(I/(P+Q))+1
IF (MOD(NODE,2) .EQ. 1) THEN
IF (I .LE. ((NODE-1)*(P+Q)+P)) THEN
Xnm(I) = 0
Ynm(I) = I-(NODE-1)*(P+Q)+(NODE-1)*P
ELSE
Xnm(I)= -(I-P-(NODE-1)*(P+Q))
Ynm(I)= P+(NODE-1)*P
END IF
ELSE
IF (I .LE. ((NODE-1)*(P+Q)+P)) THEN
Xnm(I) = -Q
Ynm(I) = I-(NODE-1)*(P+Q)+(NODE-1)*P
ELSE
Xnm(I) = (I-P-(NODE-1)*(P+Q))-Q
156
Ynm(I) = P+(NODE-1)*P
END IF
ENDIF
END DO
Xn(1)=0.
Yn(1)=0.
DO I=1,NM-1
Xn(I+1) = l0*Xnm(I)
Yn(I+1) = l0*Ynm(I)
ENDDO
DO I=1,NM
IF ((MOD((I-1),(P+Q))) .LT. P) THEN
Xm(I) = Xn(I)
Ym(I) =l0/2 + Yn(I)
ELSE
IF (MOD(INT((I-1)/(P+Q)),2) .EQ. 0) THEN
Xm(I) = -l0/2 + Xn(I)
Ym(I) = Yn(I)
ELSE
Xm(I) = l0/2 + Xn(I)
Ym(I) = Yn(I)
ENDIF
ENDIF
ENDDO
END
C *************************************
C FUNCTION TO CALCULATE RADIUS VECTOR!!!
C
DOUBLE PRECISION FUNCTION RS(NM,Xn,Yn,Xm,Ym,AA,J,K)
INTEGER J,K
DOUBLE PRECISION Xn(NM),Yn(NM),Xm(NM),Ym(NM),AA
RS = DSQRT(AA**2+(Xm(J)-Xn(K))**2+(Ym(J)-Yn(K))**2)
RETURN
END
C *************************************
SUBROUTINE DQAINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,
+MINPTS,MAXPTS,EPSABS,EPSREL,NW,KEY,RESTAR,RESULT,ABSERR,NEVAL,
+IFAIL,DONE,WORK,IWORK,FUNCT,RA,ZZ)
C Global variables.
EXTERNAL FUNSUB
INTEGER WRKSUB
INTEGER NUMFUN,MINPTS,MAXPTS,NW,RESTAR,NEVAL,IFAIL,IWORK(3*WRKSUB)
+,KEY,EMAX,FUNCT
DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,RESULT(NUMFUN),
+ABSERR(NUMFUN),WORK(NW),PERIOD,GAMMA,RA,ZZ
LOGICAL DONE(NUMFUN)
C
C Local variables.
C
INTEGER NSUB,LENW,NUM,NUMNUL
INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14
PARAMETER (NUM=21,NUMNUL=16)
C
CALL DCHINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,MINPTS,MAXPTS,EPSABS,
+EPSREL,WRKSUB,NW,EMAX,KEY,RESTAR,NEVAL,IFAIL,FUNCT)
IF (IFAIL.NE.0) THEN
RETURN
157
END IF
C
C Split up the work space.
C
I1=1
I2=I1+WRKSUB*NUMFUN
I3=I2+WRKSUB*NUMFUN
I4=I3+WRKSUB
I5=I4+WRKSUB
I6=I5+WRKSUB
I7=I6+WRKSUB*NUMFUN
I8=I7+WRKSUB*NUMFUN
I9=I8+NUMFUN
I10=I9+NUMFUN
I11=I10+NUMFUN
I12=I11+(EMAX+1)*(EMAX+1)
I13=I12+NUMFUN
I14=I13+NUMFUN
C
IF (RESTAR.EQ.1) THEN
NSUB=WORK(NW)
END IF
C
LENW=(NUMNUL+1+NUM)*NUMFUN
CALL DADINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,MINPTS,
+MAXPTS,EPSABS,EPSREL,RESTAR,LENW,KEY,RESULT,ABSERR,NEVAL,NSUB,
+IFAIL,DONE,WORK(I1),WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6),
+WORK(I7),WORK(I8),WORK(I9),WORK(I10),WORK(I11),WORK(I12),WORK(I13)
+,WORK(I14),IWORK(1),IWORK(1+WRKSUB),IWORK(1+2*WRKSUB),FUNCT,RA,ZZ)
WORK(NW)=NSUB
RETURN
END
C *****************************************
SUBROUTINE DADINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,WRKSUB,EMAX,
+MINPTS,MAXPTS,EPSABS,EPSREL,RESTAR,LENW,KEY,RESULT,ABSERR,NEVAL,
+NSUB,IFAIL,DONE,VALUES,ERRORS,GREATE,C,D,U,E,T,CT,EXTERR,BETA,
+UNEW,UERR,WORK2,WORST,LIST,UPOINT,FUNCT,RA,ZZ)
C
C Global variables.
C
EXTERNAL FUNSUB,D1MACH
INTEGER NUMFUN,LENW,RESTAR,NEVAL,WRKSUB,NSUB,IFAIL,
+LIST(WRKSUB),KEY,WORST(WRKSUB),UPOINT(WRKSUB),EMAX,
+MAXPTS,MINPTS,FUNCT
DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,RESULT(NUMFUN),
+ABSERR(NUMFUN),VALUES(NUMFUN,WRKSUB),PERIOD,ERRORS(NUMFUN,WRKSUB),
+GREATE(WRKSUB),WORK2(LENW),C(WRKSUB),D(WRKSUB),U(NUMFUN,WRKSUB),
+E(NUMFUN,WRKSUB),GAMMA,EXTERR(NUMFUN),T(NUMFUN),CT(NUMFUN),
+BETA(0:EMAX,0:EMAX),UNEW(NUMFUN),UERR(NUMFUN),D1MACH,RA,ZZ
LOGICAL DONE(NUMFUN)
C
C Local variables.
C
INTEGER I,J,K,JJ,NUMINT,ONE,NUMU,VACANT,UPDATE
INTEGER SBRGNS,I1,I2,I3,L1,L2,L3
INTEGER NDIV,POINTR,INDEX,TOP,FLAG,NUM,NUMNUL
158
DOUBLE PRECISION AOLD,BOLD,GREAT,P,LEFT,MAXEER,EPMACH,CC
LOGICAL MORE
SAVE MAXEER,NUMU,LEFT,P,CC
PARAMETER (NUMINT=2,NUM=21,NUMNUL=16)
C
C Define machine epsilon.
C
EPMACH=D1MACH(4)
C
C Set some pointers for WORK2.
C
L1=1
L2=L1+NUMNUL*NUMFUN
L3=L2+NUM*NUMFUN
C
IF (RESTAR.EQ.1) THEN
SBRGNS=NSUB
GO TO 90
ELSE
FLAG=0
CC=B(1)/PERIOD
END IF
C
C Initiate SBRGNS, DONE, MORE, P, A and B.
C
MORE=.TRUE.
SBRGNS=0
DO 10 J=1,NUMFUN
DONE(J)=.FALSE.
10 CONTINUE
P=PERIOD/2
A(2)=B(1)
B(2)=A(2)+P
LEFT=B(2)
C
C Initialize E and U
C
DO 20 J=1,NUMFUN
DO 20 I=1,WRKSUB
E(J,I)=0
20 U(J,I)=0
C
C Apply the rule procedure over the two first intervals.
C
IF (A(1).EQ.B(1)) THEN
DO 30 I=1,NUMFUN
VALUES(I,1)=0
30 ERRORS(I,1)=0
GREATE(1)=0
WORST(1)=1
ONE=2
SBRGNS=1
UPOINT(1)=1
ELSE
ONE=1
END IF
DO 40 I=ONE,2
159
CALL DGAINT (A(I),B(I),NUMFUN,FUNSUB,DONE,MORE,EPMACH,
+ WORK2(L1),WORK2(L2),WORK2(L3),FLAG,VALUES(1,I),
+ ERRORS(1,I),GREATE(I),WORST(I),C(I),D(I),FUNCT,RA,ZZ)
NEVAL=NEVAL+NUM
SBRGNS=SBRGNS+1
UPOINT(I)=I
40 CONTINUE
DO 60 I=1,2
DO 50 J=1,NUMFUN
U(J,I)=VALUES(J,I)
E(J,I)=ERRORS(J,I)
50 CONTINUE
60 CONTINUE
NUMU=2
C
DO 80 I=1,NUMINT
GREAT=0.D0
DO 70 J=1,NUMFUN
IF (ERRORS(J,I).GT.GREAT) THEN
GREAT=ERRORS(J,I)
WORST(I)=J
END IF
70 CONTINUE
GREATE(I)=GREAT
C(I)=A(I)+(B(I)-A(I))/3
D(I)=B(I)-(B(I)-A(I))/3
INDEX=I
CALL DTRINT (2,INDEX,GREATE,LIST,I)
80 CONTINUE
C
UPDATE=-1
CALL DEXINF (NUMFUN,GAMMA,CC,KEY,NUMU-1,T,CT,UPDATE,EMAX,U,E,
+RESULT,ABSERR,EXTERR,BETA,MAXEER)
C
IF (MORE.OR.(NEVAL.LT.MINPTS)) THEN
GO TO 90
END IF
IFAIL=0
GO TO 170
C
90 TOP=LIST(1)
C
IF (GREATE(TOP).LT.(MAXEER)) THEN
C
IF (NUMU-1.GT.EMAX) THEN
IFAIL=-1
GO TO 170
ELSE
VACANT=0
NDIV=1
END IF
ELSE
VACANT=TOP
NDIV=3
END IF
160
C
IF (NEVAL+NDIV*NUM.LE.MAXPTS) THEN
C
IF (VACANT.EQ.0) THEN
POINTR=SBRGNS+1
NUMU=NUMU+1
UPDATE=NUMU
INDEX=POINTR
A(INDEX)=LEFT
LEFT=LEFT+P
B(INDEX)=LEFT
TOP=INDEX
DO 100 J=1,NUMFUN
UERR(J)=0
100 UNEW(J)=0
ELSE
UPDATE=UPOINT(VACANT)
C
AOLD=A(TOP)
BOLD=B(TOP)
C
POINTR=SBRGNS+2
CALL DTRINT (1,SBRGNS,GREATE,LIST,K)
C
I1=TOP
I2=POINTR-1
I3=POINTR
A(I1)=AOLD
B(I1)=C(TOP)
A(I2)=B(I1)
B(I2)=D(TOP)
A(I3)=B(I2)
B(I3)=BOLD
INDEX=VACANT
DO 110 J=1,NUMFUN
UERR(J)=-ERRORS(J,INDEX)
110 UNEW(J)=-VALUES(J,INDEX)
END IF
C
DO 130 I=1,NDIV
CALL DGAINT (A(INDEX),B(INDEX),NUMFUN,FUNSUB,DONE,MORE,
+ EPMACH,WORK2(L1),WORK2(L2),WORK2(L3),FLAG,
+ VALUES(1,INDEX),ERRORS(1,INDEX),GREATE(INDEX),
+ WORST(INDEX),C(INDEX),D(INDEX),FUNCT,RA,ZZ)
UPOINT(INDEX)=UPDATE
C
DO 120 J=1,NUMFUN
UNEW(J)=UNEW(J)+VALUES(J,INDEX)
UERR(J)=UERR(J)+ERRORS(J,INDEX)
120 CONTINUE
INDEX=SBRGNS+I+1
130 CONTINUE
NEVAL=NEVAL+NDIV*NUM
C
DO 140 J=1,NUMFUN
U(J,UPDATE)=U(J,UPDATE)+UNEW(J)
E(J,UPDATE)=E(J,UPDATE)+UERR(J)
161
140 CONTINUE
C
IF (VACANT.EQ.0) THEN
UPDATE=-1
END IF
CALL DEXINF (NUMFUN,GAMMA,CC,KEY,NUMU-1,T,CT,UPDATE,EMAX,U,E,
+ RESULT,ABSERR,EXTERR,BETA,MAXEER)
C
MORE=.FALSE.
DO 150 J=1,NUMFUN
IF (ABSERR(J).LE.EPSREL*ABS(RESULT(J)).OR.ABSERR(J).LE.
+ EPSABS) THEN
DONE(J)=.TRUE.
ELSE
MORE=.TRUE.
DONE(J)=.FALSE.
END IF
150 CONTINUE
C
INDEX=TOP
DO 160 I=1,NDIV
JJ=SBRGNS+I
CALL DTRINT (2,JJ,GREATE,LIST,INDEX)
INDEX=SBRGNS+I+1
160 CONTINUE
SBRGNS=SBRGNS+NDIV
C
IF (MORE.OR.(NEVAL.LT.MINPTS)) THEN
GO TO 90
END IF
C
ELSE
IFAIL=+1
END IF
C
170 NSUB=SBRGNS
RETURN
END
C ***************************************************
SUBROUTINE DCHINF (NUMFUN,FUNSUB,PERIOD,GAMMA,A,B,MINPTS,MAXPTS,
+EPSABS,EPSREL,WRKSUB,NW,EMAX,KEY,RESTAR,NEVAL,IFAIL,FUNCT)
C
INTEGER NUMFUN,MINPTS,MAXPTS,NW,RESTAR,EMAX,NEVAL,WRKSUB,KEY,
+IFAIL,FUNCT
DOUBLE PRECISION A(WRKSUB),B(WRKSUB),EPSABS,EPSREL,PERIOD,GAMMA
EXTERNAL DGAINF,FUNSUB
C
C Local variables.
C
INTEGER LIMIT,NUMNUL,NUM
DOUBLE PRECISION WIDTH
PARAMETER (NUMNUL=16,NUM=21)
IFAIL=0
C
IF (NUMFUN.LT.1) THEN
IFAIL=2
162
RETURN
END IF
C
WIDTH=B(1)-A(1)
IF (WIDTH.LT.0) THEN
IFAIL=3
RETURN
END IF
C
IF (PERIOD.LE.0) THEN
IFAIL=4
RETURN
END IF
C
IF ((MAXPTS.LT.MINPTS).OR.(MINPTS.LT.NUM)) THEN
IFAIL=5
RETURN
END IF
C
C Check on legal accuracy requests.
C
IF (EPSABS.LE.0.AND.EPSREL.LE.0) THEN
IFAIL=6
RETURN
END IF
C
C Check on legal WRKSUB.
C
IF (NUM*WRKSUB.LT.MAXPTS) THEN
IFAIL=7
RETURN
END IF
C
C Check on legal value of EMAX.
C
IF (EMAX.LT.1) THEN
IFAIL=8
RETURN
END IF
C
C Check on legal RESTAR.
C
IF (RESTAR.NE.0.AND.RESTAR.NE.1) THEN
IFAIL=9
RETURN
ELSE IF (RESTAR.EQ.0) THEN
NEVAL=0
END IF
C
C Check on legal KEY.
C
IF (KEY.NE.0.AND.KEY.NE.1.AND.KEY.NE.2) THEN
IFAIL=10
RETURN
END IF
C
C Check on big enough NW.
163
C
LIMIT=1+5*NUMFUN+4*WRKSUB*NUMFUN+3*WRKSUB+(EMAX+1)**2+(NUMNUL+1+
+NUM)*NUMFUN
IF (NW.LT.LIMIT) THEN
IFAIL=11
RETURN
END IF
IF ((GAMMA.LE.0).AND.(KEY.EQ.2)) THEN
CALL DGAINF (NUMFUN,FUNSUB,B(1),PERIOD,GAMMA,NEVAL,IFAIL)
RETURN
END IF
C
RETURN
END
C *********************************
SUBROUTINE DTRINT (DVFLAG,SBRGNS,GREATE,LIST,NEW)
C
C Global variables.
C
INTEGER DVFLAG,NEW,SBRGNS,LIST(*)
DOUBLE PRECISION GREATE(*)
C
C Local variables.
C
INTEGER SUBRGN,SUBTMP,PNTR
DOUBLE PRECISION GREAT
C
C
IF (DVFLAG.NE.2) THEN
IF (DVFLAG.EQ.1) THEN
PNTR=LIST(SBRGNS)
GREAT=GREATE(PNTR)
SBRGNS=SBRGNS-1
ELSE IF (DVFLAG.EQ.3) THEN
PNTR=LIST(1)
GREAT=GREATE(PNTR)
END IF
SUBRGN=1
10 SUBTMP=2*SUBRGN
IF (SUBTMP.LE.SBRGNS) THEN
IF (SUBTMP.NE.SBRGNS) THEN
C
C Find max. of left and right child.
C
IF (GREATE(LIST(SUBTMP)).LT.GREATE(LIST(SUBTMP+1))) THEN
SUBTMP=SUBTMP+1
END IF
END IF
C
C Compare max.child with parent.
C If parent is max., then done.
C
IF (GREAT.LT.GREATE(LIST(SUBTMP))) THEN
C
C Move the pointer at position subtmp up the heap.
164
C
LIST(SUBRGN)=LIST(SUBTMP)
SUBRGN=SUBTMP
GO TO 10
END IF
END IF
C
C Update the pointer.
C
IF (SBRGNS.GT.0) THEN
LIST(SUBRGN)=PNTR
END IF
ELSE IF (DVFLAG.EQ.2) THEN
C
C If DVFLAG = 2, find the position for the NEW region in the heap.
C
GREAT=GREATE(NEW)
SUBRGN=SBRGNS
20 SUBTMP=SUBRGN/2
IF (SUBTMP.GE.1) THEN
C
C Compare max.child with parent.
C If parent is max, then done.
C
IF (GREAT.GT.GREATE(LIST(SUBTMP))) THEN
C
C Move the pointer at position subtmp down the heap.
C
LIST(SUBRGN)=LIST(SUBTMP)
SUBRGN=SUBTMP
GO TO 20
END IF
END IF
C
C Set the pointer to the new region in the heap.
C
LIST(SUBRGN)=NEW
END IF
RETURN
END
C ***********************************
SUBROUTINE DGAINT (A,B,NUMFUN,FUNSUB,DONE,MORE,EPMACH,NULL,
+FUNVAL,BASABS,FLAG,BASVAL,RGNERR,GREATE,WORST,C,D,FUNCT,RA,ZZ)
C
C Global variables.
C
EXTERNAL FUNSUB
LOGICAL DONE(NUMFUN),MORE
INTEGER NUMFUN,FLAG,WORST,FUNCT
DOUBLE PRECISION A,B,BASVAL(NUMFUN),RGNERR(NUMFUN),NULL(16,
+NUMFUN),GREATE,C,D,BASABS(NUMFUN),FUNVAL(21,NUMFUN),EPMACH,RA,ZZ
C
C Local variables.
C
INTEGER I,I1,I2,J,K,WTLENG,DGE1,DEGREE,NUMNUL,DIMNUL
PARAMETER (NUMNUL=8,DIMNUL=16)
DOUBLE PRECISION X(21),RR(NUMNUL/2-1),R,NOISE,E(NUMNUL/2),EMAX,
165
+ALPHA,CONST,W(11),ABSCIS,HLEN,SAFETY,RCRIT,ABSSUM,SUM,DIFF,CENTR,
+G(11),FACTOR,NULLW(16,11)
PARAMETER (WTLENG=11,FACTOR=0.02D+0,SAFETY=10.D0,RCRIT=0.5D0,DGE1=
+18,DEGREE=41,ALPHA=(DEGREE-DGE1)/2.D0)
C
DATA (G(I),I=1,11)/0.9937521706203895D0,0.9672268385663062D0,0.
+9200993341504008D0,0.8533633645833172D0,0.7684399634756779D0,0.
+6671388041974123D0,0.5516188358872198D0,0.4243421202074387D0,0.
+2880213168024010D0,0.1455618541608950D0,0.000000000000000D0/
C
C Weights of the 21 point Gauss quadrature rule. Same remark
C about symmetry.
C
DATA (W(I),I=1,11)/0.01601722825777433D0,0.03695378977085249D0,0.
+05713442542685720D0,0.07610011362837930D0,0.09344442345603386D0,0.
+1087972991671483D0,0.1218314160537285D0,0.1322689386333374D0,0.
+1398873947910731D0,0.1445244039899700D0,0.1460811336496904D0/
C
DATA NULLW(1,1)/0.5859224694730026D-02/
DATA NULLW(1,2)/-0.2024707000741622D-01/
DATA NULLW(1,3)/0.3883580247121445D-01/
DATA NULLW(1,4)/-0.5965412595242497D-01/
DATA NULLW(1,5)/0.8114279498343020D-01/
DATA NULLW(1,6)/-0.1019231472030305D+00/
DATA NULLW(1,7)/0.1207652963052454D+00/
DATA NULLW(1,8)/-0.1366043796711610D+00/
DATA NULLW(1,9)/0.1485698262567817D+00/
DATA NULLW(1,10)/-0.1560150459859118D+00/
DATA NULLW(1,11)/0.1585416482170856D+00/
DATA NULLW(2,1)/0.1301976799706014D-01/
DATA NULLW(2,2)/-0.4379005851020851D-01/
DATA NULLW(2,3)/0.7990096087086482D-01/
DATA NULLW(2,4)/-0.1138307201442027D+00/
DATA NULLW(2,5)/0.1394263659262734D+00/
DATA NULLW(2,6)/-0.1520456605731098D+00/
DATA NULLW(2,7)/0.1489588260146731D+00/
DATA NULLW(2,8)/-0.1296181347851887D+00/
DATA NULLW(2,9)/0.9568420420614478D-01/
DATA NULLW(2,10)/-0.5078074459100106D-01/
DATA NULLW(2,11)/0.0000000000000000D+00/
DATA NULLW(3,1)/0.2158184987561927D-01/
DATA NULLW(3,2)/-0.6965227115767195D-01/
DATA NULLW(3,3)/0.1174438965053943D+00/
DATA NULLW(3,4)/-0.1473794001233916D+00/
DATA NULLW(3,5)/0.1481989091733945D+00/
DATA NULLW(3,6)/-0.1168273220215079D+00/
DATA NULLW(3,7)/0.5890214560028095D-01/
DATA NULLW(3,8)/0.1273585156460484D-01/
DATA NULLW(3,9)/-0.8133037350927629D-01/
DATA NULLW(3,10)/0.1304777802379009D+00/
DATA NULLW(3,11)/-0.1483021322906934D+00/
DATA NULLW(4,1)/0.3119657617222001D-01/
DATA NULLW(4,2)/-0.9516116459787523D-01/
DATA NULLW(4,3)/0.1431705707841666D+00/
DATA NULLW(4,4)/-0.1462171986707822D+00/
DATA NULLW(4,5)/0.9677919508585969D-01/
166
DATA NULLW(4,6)/-0.1075583592960879D-01/
DATA NULLW(4,7)/-0.7936141880173113D-01/
DATA NULLW(4,8)/0.1380749552472930D+00/
DATA NULLW(4,9)/-0.1417577117227090D+00/
DATA NULLW(4,10)/0.8867798725104829D-01/
DATA NULLW(4,11)/0.0000000000000000D+00/
DATA NULLW(5,1)/0.4157639307601386D-01/
DATA NULLW(5,2)/-0.1179114894020921D+00/
DATA NULLW(5,3)/0.1511566572815612D+00/
DATA NULLW(5,4)/-0.1073644825716617D+00/
DATA NULLW(5,5)/0.4174411212397235D-02/
DATA NULLW(5,6)/0.1012057375471417D+00/
DATA NULLW(5,7)/-0.1472858866418607D+00/
DATA NULLW(5,8)/0.1063772962797608D+00/
DATA NULLW(5,9)/-0.2323645744823691D-02/
DATA NULLW(5,10)/-0.1030910117645103D+00/
DATA NULLW(5,11)/0.1469720414561505D+00/
DATA NULLW(6,1)/0.5246625962340516D-01/
DATA NULLW(6,2)/-0.1358182960749584D+00/
DATA NULLW(6,3)/0.1386355746642898D+00/
DATA NULLW(6,4)/-0.3967547044517777D-01/
DATA NULLW(6,5)/-0.8983329656061084D-01/
DATA NULLW(6,6)/0.1471801958801420D+00/
DATA NULLW(6,7)/-0.8524048604745531D-01/
DATA NULLW(6,8)/-0.4617298114857220D-01/
DATA NULLW(6,9)/0.1397282757969823D+00/
DATA NULLW(6,10)/-0.1185867937861861D+00/
DATA NULLW(6,11)/0.0000000000000000D+00/
DATA NULLW(7,1)/0.6363031447247006D-01/
DATA NULLW(7,2)/-0.1472028208379138D+00/
DATA NULLW(7,3)/0.1063757528761779D+00/
DATA NULLW(7,4)/0.3881687506910375D-01/
DATA NULLW(7,5)/-0.1432999618142209D+00/
DATA NULLW(7,6)/0.9698969424297650D-01/
DATA NULLW(7,7)/0.5209450556829039D-01/
DATA NULLW(7,8)/-0.1455658951943161D+00/
DATA NULLW(7,9)/0.8343286549711612D-01/
DATA NULLW(7,10)/0.6800562635441401D-01/
DATA NULLW(7,11)/-0.1465539124681842D+00/
DATA NULLW(8,1)/0.7484599067063009D-01/
DATA NULLW(8,2)/-0.1508776892345244D+00/
DATA NULLW(8,3)/0.5853283458897407D-01/
DATA NULLW(8,4)/0.1062457136342151D+00/
DATA NULLW(8,5)/-0.1318732622123368D+00/
DATA NULLW(8,6)/-0.1673118490576078D-01/
DATA NULLW(8,7)/0.1428979325476485D+00/
DATA NULLW(8,8)/-0.7818432195969258D-01/
DATA NULLW(8,9)/-0.9112601052788798D-01/
DATA NULLW(8,10)/0.1382849496064090D+00/
DATA NULLW(8,11)/0.0000000000000000D+00/
DATA NULLW(9,1)/0.8590167508061712D-01/
DATA NULLW(9,2)/-0.1462121283895959D+00/
DATA NULLW(9,3)/0.1973066649848703D-02/
DATA NULLW(9,4)/0.1434120884358229D+00/
DATA NULLW(9,5)/-0.6050045998747565D-01/
DATA NULLW(9,6)/-0.1192968264851738D+00/
DATA NULLW(9,7)/0.1063582979588903D+00/
167
DATA NULLW(9,8)/0.7871971420591506D-01/
DATA NULLW(9,9)/-0.1360664606736734D+00/
DATA NULLW(9,10)/-0.2747426951466367D-01/
DATA NULLW(9,11)/0.1463706054390223D+00/
DATA NULLW(10,1)/0.9659618304508728D-01/
DATA NULLW(10,2)/-0.1331667766592828D+00/
DATA NULLW(10,3)/-0.5483608118503819D-01/
DATA NULLW(10,4)/0.1395396581193282D+00/
DATA NULLW(10,5)/0.3842219337177220D-01/
DATA NULLW(10,6)/-0.1430606163131484D+00/
DATA NULLW(10,7)/-0.2498840938693774D-01/
DATA NULLW(10,8)/0.1451753725543706D+00/
DATA NULLW(10,9)/0.1236834040097046D-01/
DATA NULLW(10,10)/-0.1461902970879641D+00/
DATA NULLW(10,11)/0.0000000000000000D+00/
DATA NULLW(11,1)/0.1067392105869384D+00/
DATA NULLW(11,2)/-0.1122928178247447D+00/
DATA NULLW(11,3)/-0.1031959020477783D+00/
DATA NULLW(11,4)/0.9558143541939021D-01/
DATA NULLW(11,5)/0.1196951864405313D+00/
DATA NULLW(11,6)/-0.7225984000378730D-01/
DATA NULLW(11,7)/-0.1339424732473705D+00/
DATA NULLW(11,8)/0.4492456833975673D-01/
DATA NULLW(11,9)/0.1431238351576778D+00/
DATA NULLW(11,10)/-0.1523606417516131D-01/
DATA NULLW(11,11)/-0.1462742772906911D+00/
DATA NULLW(12,1)/0.1161523191050730D+00/
DATA NULLW(12,2)/-0.8469391287377601D-01/
DATA NULLW(12,3)/-0.1355896186086413D+00/
DATA NULLW(12,4)/0.2408868272651161D-01/
DATA NULLW(12,5)/0.1460359069105494D+00/
DATA NULLW(12,6)/0.4632194803727984D-01/
DATA NULLW(12,7)/-0.1231814607655799D+00/
DATA NULLW(12,8)/-0.1068762140630544D+00/
DATA NULLW(12,9)/0.7029919038187181D-01/
DATA NULLW(12,10)/0.1416700606593806D+00/
DATA NULLW(12,11)/0.0000000000000000D+00/
DATA NULLW(13,1)/0.1246701955330830D+00/
DATA NULLW(13,2)/-0.5195253287985397D-01/
DATA NULLW(13,3)/-0.1469123277046623D+00/
DATA NULLW(13,4)/-0.5433985749387003D-01/
DATA NULLW(13,5)/0.1052913655334742D+00/
DATA NULLW(13,6)/0.1341759283651172D+00/
DATA NULLW(13,7)/-0.2310968036052471D-02/
DATA NULLW(13,8)/-0.1358135230198954D+00/
DATA NULLW(13,9)/-0.1024826424015055D+00/
DATA NULLW(13,10)/0.5656562071840163D-01/
DATA NULLW(13,11)/0.1462174827723311D+00/
DATA NULLW(14,1)/0.1321420280297885D+00/
DATA NULLW(14,2)/-0.1602500237425218D-01/
DATA NULLW(14,3)/-0.1353193782985104D+00/
DATA NULLW(14,4)/-0.1170027382391999D+00/
DATA NULLW(14,5)/0.1614011431557236D-01/
DATA NULLW(14,6)/0.1330641979525424D+00/
DATA NULLW(14,7)/0.1205891076794731D+00/
DATA NULLW(14,8)/-0.8640919108146020D-02/
DATA NULLW(14,9)/-0.1294253464460428D+00/
168
DATA NULLW(14,10)/-0.1251272318395094D+00/
DATA NULLW(14,11)/0.0000000000000000D+00/
DATA NULLW(15,1)/0.1384328909795413D+00/
DATA NULLW(15,2)/0.2088816813445404D-01/
DATA NULLW(15,3)/-0.1025551817987029D+00/
DATA NULLW(15,4)/-0.1456993480020755D+00/
DATA NULLW(15,5)/-0.8041833842953963D-01/
DATA NULLW(15,6)/0.4369891359834745D-01/
DATA NULLW(15,7)/0.1355713757017371D+00/
DATA NULLW(15,8)/0.1284341564046552D+00/
DATA NULLW(15,9)/0.2777799655524739D-01/
DATA NULLW(15,10)/-0.9304002613430708D-01/
DATA NULLW(15,11)/-0.1461812140165950D+00/
DATA NULLW(16,1)/0.1434250647895144D+00/
DATA NULLW(16,2)/0.5648832390790171D-01/
DATA NULLW(16,3)/-0.5370731005946019D-01/
DATA NULLW(16,4)/-0.1320553898020212D+00/
DATA NULLW(16,5)/-0.1399117916675364D+00/
DATA NULLW(16,6)/-0.7464504070837483D-01/
DATA NULLW(16,7)/0.2922259459390760D-01/
DATA NULLW(16,8)/0.1177993871727999D+00/
DATA NULLW(16,9)/0.1454239155303997D+00/
DATA NULLW(16,10)/0.9797588771824906D-01/
DATA NULLW(16,11)/0.0000000000000000D+00/
C
CONST=SAFETY*RCRIT**(1-ALPHA)
C Compute half-length and center of interval.
C
HLEN=(B-A)/2
CENTR=A+HLEN
X(11)=CENTR
C
C Compute all abscissas
C
DO 10 I=1,10
ABSCIS=HLEN*G(I)
X(I)=CENTR-ABSCIS
X(22-I)=CENTR+ABSCIS
10 CONTINUE
IF ((X(21).EQ.X(20)).OR.(X(1).EQ.X(2))) THEN
FLAG=1
END IF
C
C Evaluate all NUMFUN functions in the 21 evaluation points.
C
DO I=1,21
DO J=1,NUMFUN
FUNVAL(I,J)=0.0
ENDDO
ENDDO
CALL FUNSUB (X,NUMFUN,21,FUNVAL,FUNCT,RA,ZZ)
C
C Apply the basic rule: first center point
C
DO J=1,NUMFUN
BASVAL(J)=0.0
169
ENDDO
DO J=1,NUMFUN
BASVAL(J)=W(11)*FUNVAL(11,J)
BASABS(J)=ABS(BASVAL(J))
ENDDO
C
C Apply all nullrules: first center point
C
DO J=1,NUMFUN
DO I=1,NUMNUL
NULL(I,J)=NULLW(I,11)*FUNVAL(11,J)
ENDDO
ENDDO
C
C Compute the basic rule contributions from the other points.
C
SUM = 0.0
DO 60 J=1,NUMFUN
DO 50 I=1,10
SUM=FUNVAL(I,J)+FUNVAL(22-I,J)
ABSSUM=ABS(FUNVAL(I,J))+ABS(FUNVAL(22-I,J))
BASVAL(J)=W(I)*SUM+BASVAL(J)
BASABS(J)=BASABS(J)+W(I)*ABSSUM
50 CONTINUE
60 CONTINUE
C
C
DO 90 J=1,NUMFUN
DO 80 K=1,10
SUM=FUNVAL(K,J)+FUNVAL(22-K,J)
DIFF=FUNVAL(K,J)-FUNVAL(22-K,J)
DO 70 I=1,NUMNUL/2
I2=2*I
I1=I2-1
NULL(I1,J)=NULL(I1,J)+NULLW(I1,K)*SUM
NULL(I2,J)=NULL(I2,J)+NULLW(I2,K)*DIFF
70 CONTINUE
80 CONTINUE
90 CONTINUE
C
C Scale the results with the length of the interval
C
DO 100 J=1,NUMFUN
BASVAL(J)=HLEN*BASVAL(J)
BASABS(J)=HLEN*BASABS(J)
100 CONTINUE
C
C Compute errors.
C
GREATE=0
DO 130 J=1,NUMFUN
EMAX=0
DO 110 I=1,NUMNUL/2
I2=2*I
I1=I2-1
170
E(I)=HLEN*DSQRT(NULL(I1,J)**2+NULL(I2,J)**2)
EMAX=MAX(EMAX,E(I))
110 CONTINUE
R=0
DO 120 I=1,NUMNUL/2-1
IF (E(I+1).NE.0) THEN
RR(I)=E(I)/E(I+1)
ELSE
RR(I)=1
END IF
R=MAX(R,RR(I))
120 CONTINUE
IF (R.GE.1) THEN
RGNERR(J)=SAFETY*EMAX
ELSE IF (R.GE.RCRIT) THEN
RGNERR(J)=SAFETY*R*E(1)
ELSE
RGNERR(J)=CONST*(R**ALPHA)*E(1)
END IF
C
C Check the noise level.
C
NOISE=50*EPMACH*BASABS(J)
IF ((E(1).LT.NOISE).AND.(E(2).LT.NOISE)) THEN
RGNERR(J)=NOISE
ELSE
RGNERR(J)=MAX(NOISE,RGNERR(J))
END IF
C
IF (.NOT.(MORE.AND.DONE(J)).AND.(RGNERR(J).GT.GREATE)) THEN
GREATE=RGNERR(J)
WORST=J
END IF
130 CONTINUE
C
C Compute the new subdivision points.
C
C=A+(B-A)/3
D=B-(B-A)/3
RETURN
END
C ****************************************
SUBROUTINE DGAINF (NUMFUN,F,B,PERIOD,GAMMA,NEVAL,IFAIL)
C
C
C Global variables.
C
EXTERNAL F,D1MACH
DOUBLE PRECISION PERIOD,GAMMA,B,D1MACH
INTEGER NUMFUN,IFAIL,NEVAL
C
C Local variables.
C
DOUBLE PRECISION P,X,H,XB,FX,FXP,XR,FR,EPMACH,FOLD,FNEW,SIGN,XOLD,
+XNEW
INTEGER I,IMX,M
PARAMETER (IMX=15)
171
DOUBLE PRECISION A(0:IMX),DA,FACTOR(0:IMX),FI(0:IMX),DF,LAST
EPMACH=D1MACH(4)
C
XOLD=B
CALL F (XOLD,1,1,FOLD)
P=PERIOD/4
XNEW=XOLD
DO 10 I=1,3,1
XNEW=XNEW+P
CALL F (XNEW,1,1,FNEW)
IF (ABS(FOLD).LT.ABS(FNEW)) THEN
XOLD=XNEW
IF ((FOLD*FNEW).GT.0) THEN
P=-P
END IF
FOLD=FNEW
ELSE
IF ((FOLD*FNEW).GT.0) THEN
XNEW=XOLD
END IF
END IF
P=P/2
10 CONTINUE
X=XOLD
FX=FOLD
C
CALL F (X+PERIOD/2,1,1,FXP)
NEVAL=5
LAST=FXP
SIGN=FX
FI(0)=X
IF (((SIGN*FXP).LT.0).AND.(ABS(FX).GT.ABS(FXP))) THEN
FR=-FXP/FX
XR=X/(X+PERIOD/2)
A(0)=LOG(FR)/LOG(XR)
FACTOR(0)=A(0)*(1/ABS(FR*LOG(FR))+1/ABS(XR*LOG(XR)))*EPMACH
ELSE
IFAIL=12
RETURN
END IF
C
H=PERIOD/2
XB=X
DO 30 I=1,IMX
H=2*H
X=XB+H
FI(I)=X
CALL F (X+PERIOD/2,1,1,FXP)
CALL F (X,1,1,FX)
NEVAL=NEVAL+2
C
C Check the computed function values.
C
IF (((SIGN*FX).GT.0).AND.((SIGN*FXP).LT.0).AND.(ABS(FX).GT.
+ ABS(FXP)).AND.(ABS(LAST).GT.ABS(FX))) THEN
LAST=FXP
172
FR=-FXP/FX
XR=X/(X+PERIOD/2)
A(I)=LOG(FR)/LOG(XR)
FACTOR(I)=A(0)*(2**I/ABS(FR*LOG(FR))+1/ABS(XR*LOG(XR)))*
+ EPMACH
ELSE
C
C Either PERIOD is wrong or B is too small.
C
IFAIL=12
RETURN
END IF
DO 20 M=I-1,0,-1
DA=(A(M+1)-A(M))/((FI(I)-FI(M))/FI(M))
DF=(FACTOR(M+1)+FACTOR(M))/((FI(I)-FI(M))/FI(M))
A(M)=A(M+1)+DA
FACTOR(M)=FACTOR(M+1)+DF
20 CONTINUE
IF (ABS(DA).LE.FACTOR(0)) THEN
GAMMA=A(0)
C
C Normal return
C
RETURN
END IF
30 CONTINUE
C
C We did not succeed in estimating GAMMA to sufficient precision.
C
GAMMA=A(0)
IFAIL=13
RETURN
END
C *************************************
SUBROUTINE DEXINF (NUMFUN,GAMMA,C,KEY,N,T,CT,UPDATE,EMAX,U,E,
+RESULT,ABSERR,EXTERR,BETA,MAXEER)
C
C Global variables.
C
INTEGER N,EMAX,UPDATE,NUMFUN,KEY
DOUBLE PRECISION GAMMA,T(NUMFUN),CT(NUMFUN),C,U(NUMFUN,0:N),
+E(NUMFUN,0:N),BETA(0:EMAX,0:EMAX),MAXEER,RESULT(NUMFUN),
+ABSERR(NUMFUN),EXTERR(NUMFUN)
C
C Local variables.
C
INTEGER I,J
DOUBLE PRECISION SAVE1,SAVE2,CONST,BASE1,BASE2,P,Q
PARAMETER (CONST=1.D0)
BETA(0,0)=1
IF (KEY.EQ.2) THEN
BASE1=GAMMA
BASE2=MAX(GAMMA-1,4*C)
P=2.D0
ELSE
BASE1=0.D0
173
BASE2=MAX(0.D0,4*C)
P=1.D0
END IF
C
C Compute the new extrapolation coefficients.
C
UPDATE=-1
IF (UPDATE.LT.0) THEN
BETA(0,N)=1
BETA(N,0)=1
DO 20 I=1,N-1
SAVE1=1
DO 10 J=N-I,N-1
IF (KEY.EQ.0) THEN
Q=0.5D0
ELSE
Q=(1+(-BASE1-P*J)/(2*N+BASE2))/2.D0
END IF
SAVE2=SAVE1-(SAVE1-BETA(I,J))*Q
BETA(I,J)=SAVE1
10 SAVE1=SAVE2
BETA(I,N)=SAVE1
20 CONTINUE
DO 30 J=1,N
IF (KEY.EQ.0) THEN
Q=0.5D0
ELSE
Q=(1+(-BASE1-P*(J-1))/(2*N+BASE2))/2.D0
END IF
30 BETA(N,J)=(1-Q)*BETA(N,J-1)
END IF
C
C Extrapolate; use the extrapolation coefficients.
C
DO 40 J=1,NUMFUN
T(J)=0
CT(J)=0
DO 40 I=N,0,-1
CT(J)=CT(J)+BETA(I,N-1)*U(J,I)
40 T(J)=T(J)+BETA(I,N)*U(J,I)
C
MAXEER=0
DO 60 J=1,NUMFUN
EXTERR(J)=CONST*ABS(T(J)-CT(J))/Q
MAXEER=MAX(EXTERR(J),MAXEER)
C
C Note: The last U-errors are effected by the extrapolation-process
C Accumulate the total error in exterr.
C
DO 50 I=0,N
50 EXTERR(J)=EXTERR(J)+BETA(I,N)*E(J,I)
60 CONTINUE
C
C Define the results: update the results only if the error is improved
C
DO 70 J=1,NUMFUN
174
C IF ((EXTERR(J).LE.ABSERR(J)).OR.(ABSERR(J).EQ.0)) THEN
RESULT(J)=T(J)
ABSERR(J)=EXTERR(J)
C END IF
70 CONTINUE
RETURN
END
C ***************************************
DOUBLE PRECISION FUNCTION D1MACH(I)
C
C
INTEGER SMALL(4)
INTEGER LARGE(4)
INTEGER RIGHT(4)
INTEGER DIVER(4)
INTEGER LOG10(4)
INTEGER SC,I
C
DOUBLE PRECISION DMACH(5)
C
EQUIVALENCE (DMACH(1),SMALL(1))
EQUIVALENCE (DMACH(2),LARGE(1))
EQUIVALENCE (DMACH(3),RIGHT(1))
EQUIVALENCE (DMACH(4),DIVER(1))
EQUIVALENCE (DMACH(5),LOG10(1))
C
DATA SMALL(1),SMALL(2) / 1048576, 0 /
DATA LARGE(1),LARGE(2) / 2146435071, -1 /
DATA RIGHT(1),RIGHT(2) / 1017118720, 0 /
DATA DIVER(1),DIVER(2) / 1018167296, 0 /
DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/
C
C *** ISSUE STOP 779 IF ALL DATA STATEMENTS ARE COMMENTED...
IF (SC .NE. 987) STOP 779
IF (I .LT. 1 .OR. I .GT. 5) GOTO 999
D1MACH = DMACH(I)
RETURN
999 WRITE(*,1999) I
1999 FORMAT(' D1MACH - I OUT OF BOUNDS',I10)
STOP
END
175
LÝ LỊCH KHOA HỌC
I. LÝ LỊCH SƠ LƯỢC:
Họ và tên: Trần Minh Tuấn Giới tính: Nam
Ngày sinh: 20 tháng 9 năm 1970 Nơi sinh: Hà Nội
Quê quán: Xã Tây Mỗ, huyện Từ Liêm, Hà Nội
Dân tộc: Kinh
Chức vụ, đơn vị công tác: Trưởng Ban Cơ sở hạ tầng thông tin và
Công nghiệp, Viện Chiến lược BCVT và CNTT, Bộ BCVT.
Nơi ở hiện nay: 93 phố Trúc Bạch, quận Ba Đình, Hà Nội
Điện thoại cơ quan: + 84 4 556 53; Fax:+ 84 4 556 73 99
E-mail: tm_tuan@mpt.gov.vn
II. QUÁ TRÌNH ĐÀO TẠO:
1. Đại học:
Hệ đào tạo: Chính quy. Thời gian đào tạo: Từ tháng 10/1989 đến 10/1994.
Nơi học: Trường Đại học Công nghệ Moscow, Moscow, CHLB Nga.
Ngành học: Kỹ thuật vô tuyến và hệ thống vô tuyến
Các môn thi quốc gia chủ yếu: 4 (bốn) môn: Mạch điện tử và Xử lý tín hiệu; Lý thuyết điện
động học và truyền sóng; Thông tin viba, vệ tinh và cáp quang; Truyền hình.
Ngày và nơi bảo vệ luận án: Tháng 8/1994 tại Moscow, CHLB Nga.
Người hướng dẫn: TS. Atiyshenko A. A.
2. Thạc sĩ:
Thời gian đào tạo: Từ tháng 10/1994 đến 7/1995.
Nơi học: Trường Đại học Công nghệ Moscow, Moscow, CHLB Nga.
Ngành học: Điện tử - Viễn thông
Luận án: Thiết kế, xây dựng và khai thác hệ thống truyền hình cáp cho Việt Nam
Ngày và nơi bảo vệ luận án: Tháng 6/1995 tại Moscow, CHLB Nga.
Người hướng dẫn: TS. Atiyshenko A. A.
3. Tiến sĩ:
Thời gian đào tạo: Từ tháng 10/1999 đến 01/2004.
Nơi học: Trường Đại học Bách Khoa Hà Nội, Hà Nội, Việt Nam.
Ngành học: Điện tử - Viễn thông
Luận văn: Ứng dụng phương pháp moment trong bài toán phân tích các cấu trúc điện từ
phẳng được kích thích bởi sóng chạy.
Ngày và nơi bảo vệ luận văn: 08/01/2004 tại Đại học Bách Khoa Hà Nội.
Người hướng dẫn: GS. TSKH. Phan Anh
4. Trình độ ngoại ngữ:
- Tiếng Nga: Tốt nghiệp Đại học và Thạc sĩ tại CHLB Nga
- Tiếng Anh: Hai Chứng chỉ tiếng Anh trình độ trên C (Upper-Intermediate) cấp tại Học
viện Công nghệ Hoàng gia Melbourne (Úc) các năm 1996 và 1997.
- Tiếng Pháp: Trình độ B
5. Học vị, học hàm, chức vụ kỹ thuật được chính thức cấp, số bằng, ngày và nơi cấp:
- Chứng chỉ hành nghề kỹ sư bảo dưỡng hệ thống truyền hình cáp (bao gồm hệ thống
cáp quang, cáp đồng và truyền hình cáp) cấp tại Nhà máy Kỹ thuật Điện Mitishi,
Moscow, CHLB Nga ngày 10/11/1993.
- Chứng chỉ hành nghề kỹ sư khai thác, bảo dưỡng các thiết bị thông tin vệ tinh cấp tại
Trung tâm Thông tin Vệ tinh Dubna (Moscow, CHLB Nga) ngày 25/11/1994.
- Chứng chỉ cho phép hành nghề lắp đặt, điều chỉnh, đo đạc, bảo dưỡng và sửa chữa các
thiết bị thu phát vệ tinh phục vụ tập thể và cá nhân cấp tại Trung tâm Công nghiệp -
Khoa học Hàng không Vũ trụ Quân sự CROSNA (Moscow, CHLB Nga) ngày
12/3/1995.
176
- Hai Chứng chỉ tiếng Anh trình độ trên C (Upper-Intermediate) cấp tại Học viện Công
nghệ Hoàng gia Melbourne (Úc) các năm 1996 và 1997.
- Học vị: Tiến sỹ Kỹ thuật từ năm 2004.
III. QUÁ TRÌNH CÔNG TÁC CHUYÊN MÔN KỂ TỪ KHI TỐT NGHIỆP ĐẠI HỌC
Thời gian Nơi công tác Công việc đảm nhiệm
Từ 11/1995
đến 02/1996:
Công ty Siemens
Việt Nam. Ngõ Mai
Hương, Quận Hai
Bà Trưng, Hà Nội.
Chuyên viên kỹ thuật
- Lắp đặt, nghiệm thu và vận hành hệ thống thiết bị BTS của
Motorola và BSC của Siemens thuộc mạng Vinaphone tại Hà
Nội và một số tỉnh phía Bắc.
Từ 02/1996
đến 10/1999:
Công ty Viễn thông
Quốc tế, Tổng Công
ty Bưu chính - Viễn
thông Việt Nam
Chuyên viên kỹ thuật
- Tham gia nghiên cứu và lập dự án tiền khả thi, dự án khả thi,
thiết kế, lập dự toán, chấm thầu dự án: Tuyến cáp quang biển
quốc tế SEA-ME-WE-3 nối liền Đông Bắc Á - Đông Nam Á -
Úc - Ấn độ, Trung Đông và Tây Âu, trong đó phần tham gia
của phía Việt nam là 31,4 triệu USD trong tổng giá thành dự
án là 1,2 tỷ USD với sự tham gia của 90 hãng viễn thông thuộc
41 nước. Đã đưa vào khai thác từ năm 1998.
- Tham gia nghiên cứu và lập dự án tiền khả thi triển khai hệ
thống thông tin di động sử dụng vệ tinh tầm thấp
GLOBALSTARS tại Việt Nam với tổng giá trị đầu tư là 21
triệu USD.
- Tham gia nghiên cứu và lập dự án khả thi đầu tư tuyến cáp
quang CSC nối liền 6 nước: Trung Quốc - Việt Nam - Lào -
Thái Lan - Malaysia - Singapore với tổng giá trị dự án khoảng
13 triệu USD, gồm tuyến cáp và 3 trạm cáp lớn tại Lạng Sơn
(nối sang Trung Quốc), Hà Nội (hệ thống theo dõi và điều
khiển toàn tuyến), Vinh (nối sang Lào). Đã hoàn thành xong và
đưa vào khai thác giai đoạn 1 chuyển lưu lượng từ Việt Nam đi
Trung Quốc và Lào, đang triển khai giai đoạn 2 là các tuyến vu
hồi từ Hà Nội đi Vinh và Đông Hà, Quảng Trị.
- Tham gia lập dự án khả thi, lập thiết kế dự toán và chấm thầu
dự án: Xây dựng hệ thống tính cước tập trung cho Công ty
Viễn thông Quốc tế với tổng giá trị dự án là 2,4 triệu USD
trong đó: 500.000 USD đầu tư thiết bị phần cứng bao gồm các
máy chủ, router, switch... và 1,9 triệu USD đầu tư phần mềm
xử lý, thu thập số liệu cước và xây dựng các báo cáo quản lý.
Hệ thống đang hoạt động hiệu quả tại Công ty Viễn thông
Quốc tế.
- Tham gia lập dự án khả thi, lập thiết kế dự toán, chấm thầu,
theo dõi thi công và nghiệm thu dự án: Xây dựng mạng quản
lý nội bộ LAN cho Văn phòng Công ty Viễn thông Quốc tế với
tổng giá trị dự án là 80.000 - 100.000 USD. Hệ thống mạng đã
hoạt động tốt.
- Tham gia lập dự án khả thi, thiết kế dự toán với Công ty
VDC trong dự án: Xây dựng mạng Frame Relay. Dịch vụ đã đi
vào hoạt động hiệu quả.
- Tham gia tổ chức thông tin vệ tinh quốc tế INTELSAT,
INTERSPUTNIK với chức danh quản lý khai thác của phía
Việt Nam.
Từ 10/1999
đến 10/2000:
Công ty Viễn thông
Quốc tế, Tổng Công
ty Bưu chính - Viễn
thông Việt Nam
Tổ trưởng Tổ Tổng hợp
- Tham gia tổ chức thông tin vệ tinh quốc tế INTELSAT,
INTERSPUTNIK với chức danh quản lý khai thác của phía
Việt Nam.
Từ 10/2000
đến 3/2002:
Công ty Viễn thông
Quốc tế, Tổng Công
ty Bưu chính - Viễn
Phó phòng phụ trách Phòng Tổng hợp
- Tham gia tổ chức thông tin vệ tinh quốc tế INTELSAT,
INTERSPUTNIK với chức danh quản lý khai thác của phía
177
thông Việt Nam Việt Nam.
- Tham gia lập dự án khả thi, lập thiết kế dự toán và tổ chức
đấu thầu dự án: Xây dựng hệ thống cung cấp dịch vụ VoIP
quốc tế, với tổng giá trị hệ thống ban đầu là 900.000 USD kết
nối với hãng MCI-Worldcom. Đã đưa vào khai thác từ tháng
10/2001.
- Tham gia lập dự án tiền khả thi đầu tư vào hệ thống cáp
quang biển C2C với tổng số vốn tham gia của phía Việt Nam
là 38 triệu USD.
Từ 3/2002
đến 3/2003:
Công ty Viễn thông
Quốc tế, Tổng Công
ty Bưu chính - Viễn
thông Việt Nam
Trưởng Phòng Tổng hợp
- Tham gia tổ chức thông tin vệ tinh quốc tế INTELSAT,
INTERSPUTNIK với chức danh quản lý khai thác của phía
Việt Nam.
- Tham gia lập dự án khả thi đầu tư vào hệ thống cáp quang
biển C2C với tổng số vốn tham gia của phía Việt Nam là 38
triệu USD.
Từ 3/2003
đến nay:
Viện Chiến lược
Bưu chính, Viễn
thông và CNTT, Bộ
Bưu chính, Viễn
thông, 115 Trần Duy
Hưng, Hà Nội
Trưởng Ban Cơ sở hạ tầng thông tin và công nghiệp,
- Xây dựng quy hoạch phát triển viễn thông và Internet Việt
Nam đến 2010.
- Xây dựng Chiến lược phát triển Công nghệ Thông tin Việt
Nam đến năm 2010.
- Chủ trì nhiều đề tài cấp Bộ và xây dựng các quy hoạch phát
triển BCVT và CNTT cho các địa phương trong cả nước.
IV. CÁC CÔNG TRÌNH KHOA HỌC ĐÃ CÔNG BỐ:
1. Đề tài nghiên cứu khoa học cấp Công ty:
Nội dung: Nghiên cứu công nghệ Internet Telephony và ảnh hưởng của công nghệ này tới
doanh thu viễn thông quốc tế.
Hoàn thành: Tháng 5/1998. Nhận xét: Đề tài được đánh giá kết quả tốt, mang tính cấp thiết,
đáp ứng được nhu cầu sản xuất kinh doanh của Công ty.
2. Đề tài nghiên cứu khoa học cấp Công ty:
Nội dung: Xây dựng chiến lược con người viễn thông quốc tế trong thập kỷ đầu thế kỷ 21.
Hoàn thành: Tháng 10/2002. Nhận xét: Đề tài được đánh giá kết quả tốt, mang tính cấp thiết,
đáp ứng được nhu cầu xây dựng chiến lược và đào tạo con người của Công ty.
3. Công trình nghiên cứu được đưa trên tạp chí IEEE của Hoa Kỳ:
Tên bài báo: On a method of formation of surface wave structures using for strip antennas.
Tác giả: GS. TSKH. Phan Anh và Th.S. Trần Minh Tuấn, Tạp chí: Antennas and Propagation
Society, 2001 IEEE International Sym, Năm: 2001, Vol: 3 pages: 720, Nhà xuất bản
(publisher): IEEE provider. Internet: IEEE.
/7598/20721/00960197.pdf
4. Công trình nghiên cứu được báo cáo tại Hội nghị về sóng siêu cao tần châu Á - Thái
Bình Dương (Asia-Pacific Microwave Conference - APMC) tổ chức tại Hàn Quốc năm
2003:
a) Tên công trình: Moment Method Formulation for Traveling-Wave Excitation of
Cavity-Backed Arbitrary Narrow Slot Structures.
Tác giả: GS. TSKH. Phan Anh và Th.S. Trần Minh Tuấn
b) Tên công trình: Some Studies on Scattering Properties of Traveling-Wave Excitation
Microstrip Structures by Moment Method.
Tác giả: GS. TSKH. Phan Anh và Th.S. Trần Minh Tuấn
5. H i ngh Vô tuy n & i n t l n th 9 – REV 2004, Hà N i, 26 - 27/11/2004
Tên công trình: Some studies on scattering properties of plannar structures to simulate
arbitrary ones in building antennas with designed parameters.
Tác giả: TS. Trần Minh Tuấn, GS. TSKH Phan Anh
5. Xây dựng quy hoạch phát triển viễn thông và Internet Việt Nam đến năm 2010:
178
Ngày 07/2/2006 của Thủ tướng chính phủ đã ký quyết định số 32/2006/QĐ-TTg phê duyệt
quy hoạch phát triển Viễn thông và Internet Việt Nam đến năm 2010.
6. Các đề tài khoa học công nghệ:
Chủ trì hoàn thành 3 đề tài nghiên cứu KHCN cấp Bộ trong giai đoạn 2003 – 2006:
a) Tên đề tài: “Nghiên cứu định hướng xây dựng chính sách và khuôn khổ pháp lý về Bưu
chính, Viễn thông và CNTT”. Mã số: 77-03-KHKT-QL. Kết quả: Khá.
b) Tên đề tài: Nghiên cứu đề xuất các quy định pháp lý về cạnh tranh, chống độc quyền,
chống bán phá giá trong lĩnh vực Bưu chính, Viễn thông. Mã số: 41-04-KHKT-RD. Kết
quả: Tốt.
c) Tên đề tài: “Nghiên cứu các hình thức phân tách mạch vòng nội hạt và đề xuất cơ chế cho
Việt Nam”. Mã số: 60-05-KHKT-RD. Kết quả: Tốt.
7. Xây dựng quy hoạch BCVT và CNTT cho các địa phương:
Giai đoạn 2005 - 2006 trực tiếp chỉ đạo và thực hiện quy hoạch về Bưu chính, Viễn thông và
CNTT cho 23 tỉnh và thành phố sau: Vĩnh Phúc, Phú Thọ, Hà Tĩnh, Thanh Hóa, Huế, Bình
Thuận, Đà Nẵng, Đồng Nai, Hà Nam, Khánh Hòa, Bình Định, Bắc Giang, Lâm Đồng, Long
An, Quảng Ngãi, Nam Định, Hà Tây, Lào Cai, Cao Bằng, Nghệ An, Đăk Lăk, Ninh Bình,
Hòa Bình.
Trên đây là tóm tắt lý lịch khoa học của tôi. Tôi xin cam đoan những lời khai trên là đúng.
Nếu sai tôi hoàn toàn chịu trách nhiệm.
Hà Nội, ngày 20 tháng 7 năm 2006
Người khai
Trần Minh Tuấn
Các file đính kèm theo tài liệu này:
Ứng dụng phương pháp moment trong bài toán phân tích các kết cấu điện từ phẳng được kích thích bởi sóng chạy.pdf
