Fortran Coder

标题: 求谐波 [打印本页]

作者: wawewen    时间: 2020-1-10 15:09
标题: 求谐波
本帖最后由 wawewen 于 2020-1-15 11:24 编辑

想写一个求谐波的一个子函数。k表示谐波次数,i是行数,,j是列数,Bz(m,n-1)是纵向沿角度分布,横向沿半径分布的一组数据,现在就想每个半径下沿角度求积分,算出来结果不对,有没有大佬纠正一下
[Fortran] 纯文本查看 复制代码
      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  

作者: wawewen    时间: 2020-1-15 10:07
各位大佬帮帮忙啊
作者: necrohan    时间: 2020-1-15 14:16
这属于逻辑问题,不清楚原理的人没法帮忙。
作者: wawewen    时间: 2020-1-16 11:02
necrohan 发表于 2020-1-15 14:16
这属于逻辑问题,不清楚原理的人没法帮忙。

唉,问题找到了,看了另一个帖子学了一下单步调试,sum1和sum2没有置零,这么简单的错误找了好久,果然我还是太菜了啊
作者: fcode    时间: 2020-1-17 08:08
wawewen 发表于 2020-1-16 11:02
唉,问题找到了,看了另一个帖子学了一下单步调试,sum1和sum2没有置零,这么简单的错误找了好久,果然我 ...

学了单步调试,你进步就很大了。
此后,很多问题你都可以自己排查了。
作者: wawewen    时间: 2020-1-17 10:19
fcode 发表于 2020-1-17 08:08
学了单步调试,你进步就很大了。
此后,很多问题你都可以自己排查了。

实不相瞒,看的正是你的某一个回帖的链接,你们几个大佬对我的学习还是有很大的帮助的




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2