求谐波
本帖最后由 wawewen 于 2020-1-15 11:24 编辑想写一个求谐波的一个子函数。k表示谐波次数，i是行数，，j是列数，Bz(m,n-1)是纵向沿角度分布，横向沿半径分布的一组数据，现在就想每个半径下沿角度求积分，算出来结果不对，有没有大佬纠正一下
program harmonic
implicit none
integer i,j,jj
integer,parameter::m=721,n=604,l=10
real B(m,n),Bz(m,n-1),r(n-1),theta(m)
open(unit=10,file="C:\Users\zklz\Desktop\Bz_T73&T72.txt")
!打开磁场MAP文件tosca格式
do i=1,m
read(10,*)(B(i,j),j=1,n) !读入全部数据
end do
close(10)
!读入theta并输出
do i=1,m
if(i<m)then
theta(i)=B(i+1,1)
else
theta(i)=360.0
end if
end do
open(unit=20,file="C:\Users\zklz\Desktop\theta.txt")
do i=1,m
write(20,"(f6.1)")theta(i)
end do
close(20)
!读入r并输出
do j=1,n
if(j==1)cycle
r(j-1)=B(1,j)
end do
open(unit=30,file="C:\Users\zklz\Desktop\r.txt")
do jj=1,n-1
write(30,"(f4.1)")r(jj)
end do
close(30)
!读入Bz,输出BMAP
do i=1,m
do jj=1,n-1
if(i<m)then
Bz(i,jj)=B(i+1,jj+1)
else
Bz(i,jj)=B(2,jj+1)
end if
end do
end do
open(unit=40,file="C:\Users\zklz\Desktop\BMAP.txt")
do jj=1,n-1
write(40,"(12X,f4.1)")r(jj)
write(40,"(10f12.4)")(Bz(i,jj),i=1,m-1) !不输出360度的值
end do
close(40)
call Bnharm(Bz,theta)
end
subroutine Bnharm(Bz,theta)
implicit none
integer i,j,k
integer,parameter::m=721,n=604,l=10
real an(l,n-1),bn(l,n-1),Bzn(l,n-1),Bz(m,n-1),theta(m),
&sa(l,m-1,n-1),sb(l,m-1,n-1)
real,parameter::pi=3.141592654
real::sum1=0.0,sum2=0.0
do k=1,l
do j=1,n-1
do i=1,m-1
sa(k,i,j)=(Bz(i,j)*cosd(k*theta(i))+Bz(i+1,j)
&*cosd(k*theta(i+1)))*0.25
sb(k,i,j)=(Bz(i,j)*sind(k*theta(i))+Bz(i+1,j)
&*sind(k*theta(i+1)))*0.25
end do
end do
end do
do k=1,l
do j=1,n-1
do i=1,m-1
sum1=sum1+sa(k,i,j)
sum2=sum2+sb(k,i,j)
end do
an(k,j)=sum1/180.0
bn(k,j)=sum2/180.0
Bzn(k,j)=sqrt(an(k,j)**2+bn(k,j)**2)
end do
end do
open(unit=50,file="C:\Users\zklz\Desktop\Bn.txt")
do k=1,l
!write(50,"(603f18.5)" )(an(k,j),j=1,n-1)
write(50,"(603f18.5)" )(Bzn(k,j),j=1,n-1)
end do
close(50)
end 各位大佬帮帮忙啊:'( 这属于逻辑问题，不清楚原理的人没法帮忙。 necrohan 发表于 2020-1-15 14:16
这属于逻辑问题，不清楚原理的人没法帮忙。
唉，问题找到了，看了另一个帖子学了一下单步调试，sum1和sum2没有置零，这么简单的错误找了好久，果然我还是太菜了啊 wawewen 发表于 2020-1-16 11:02
唉，问题找到了，看了另一个帖子学了一下单步调试，sum1和sum2没有置零，这么简单的错误找了好久，果然我 ...
学了单步调试，你进步就很大了。
此后，很多问题你都可以自己排查了。 fcode 发表于 2020-1-17 08:08
学了单步调试，你进步就很大了。
此后，很多问题你都可以自己排查了。
实不相瞒，看的正是你的某一个回帖的链接，你们几个大佬对我的学习还是有很大的帮助的:-lol
页:
[1]