谁能帮我看个矩量法的程序

fotran写的，算一个圆的2维散射体，先解出电流，再算散射截面。没有语法错误了，但是运行以后有错误，大概是数值范围益处，我查不出来了，急用！感谢死了！
program homework1
use imsl
implicit none
integer::maxnodes,maxeles
real r,FIi,k,z,deltaCn,aii,pi
integer i, j, m
real,allocatable::xy(:,:)
integer,allocatable::ln(:,:)
complex,allocatable::x(:),a(:,:),b(:)
real delta, H, L
real gama, e
real FI
real Taw
complex jj
complex q
write(*,*) "please input r:"
read(*,*) r
write(*,*) "please input FIi:"
read(*,*) FIi
write(*,*) "please inpute maxnodes:"
read(*,*) maxnodes
write(*,*) "please input k:"
read(*,*) k
write(*,*) "please input z:"
read(*,*) z
allocate(xy(2,maxnodes),ln(2,maxeles),a(maxeles, maxeles),x(maxeles),b(maxeles))
jj=cmplx(0.,1.)
maxeles = maxnodes
PI = 3.14
deltaCn = 2 * PI * r / maxeles
do i = 1, maxnodes
    aii = 2 * PI * i / maxnodes
    xy(1, i) = r * cos(aii)
    xy(2, i) = r * sin(aii)
end do
do i = 1, maxeles-1
ln(1, i) = i
ln(2, i) = i+1
end do
ln(1, maxeles)=maxeles
ln(2, maxeles)=1
do i = 1, maxeles
    do j = 1, maxeles
        if(i .ne. j)  then 
            delta = k * sqrt( (xy(1,i)-xy(1,j))**2 + (xy(2,i)-xy(2,j))**2 )
            call bsj0(delta)
            H=bsj0(delta)
            call bsy0(delta)
            L= bsy0(delta)
            a(i,j) = z*k*deltaCn*(H+jj*L)/4
        else 
            gama = 1.718
            e = 2.718
            m = log10(gama*k*deltaCn/4/e)
            a(i,j) =z*k*deltaCn/4*(1-jj*2*m/PI)
        end if
    end do
end do
do j = 1, maxeles
    b(i) = exp(-jj*k*(xy(1,j)*cos(FIi)+xy(2,j)*sin(FIi)))
end do
call lsacg(maxeles, a, maxeles, b, 1, x)
q = 0
write(*,*) "please input FI:"
read(*,*) FI
do i = 1, maxeles
    q = q + x(i) * exp(jj*( xy(1,i)*FI + xy(2,i)*FI) )
end do
Taw = k * z * z * ( q * conjg(q) )
write(*,*) "Taw="
print*, Taw 
pause
end program

申明：网友回复良莠不齐，仅供参考。如需专业解答，请学习本站推出的微波射频专业培训课程。