Fortran Coder

查看: 588|回复: 5

[其他行业算法] 求谐波

[复制链接]

21

帖子

7

主题

0

精华

熟手

F 币
169 元
贡献
93 点

新人勋章

发表于 2020-1-10 15:09:14 | 显示全部楼层 |阅读模式
本帖最后由 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  
心怀不惧,方能翱翔于天际
回复

使用道具 举报

21

帖子

7

主题

0

精华

熟手

F 币
169 元
贡献
93 点

新人勋章

 楼主| 发表于 2020-1-15 10:07:20 | 显示全部楼层
各位大佬帮帮忙啊
心怀不惧,方能翱翔于天际

135

帖子

2

主题

0

精华

大师

F 币
885 元
贡献
461 点

规矩勋章

发表于 2020-1-15 14:16:53 | 显示全部楼层
这属于逻辑问题,不清楚原理的人没法帮忙。

21

帖子

7

主题

0

精华

熟手

F 币
169 元
贡献
93 点

新人勋章

 楼主| 发表于 2020-1-16 11:02:09 | 显示全部楼层
necrohan 发表于 2020-1-15 14:16
这属于逻辑问题,不清楚原理的人没法帮忙。

唉,问题找到了,看了另一个帖子学了一下单步调试,sum1和sum2没有置零,这么简单的错误找了好久,果然我还是太菜了啊
心怀不惧,方能翱翔于天际

1543

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1454 元
贡献
956 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2020-1-17 08:08:31 | 显示全部楼层
wawewen 发表于 2020-1-16 11:02
唉,问题找到了,看了另一个帖子学了一下单步调试,sum1和sum2没有置零,这么简单的错误找了好久,果然我 ...

学了单步调试,你进步就很大了。
此后,很多问题你都可以自己排查了。

21

帖子

7

主题

0

精华

熟手

F 币
169 元
贡献
93 点

新人勋章

 楼主| 发表于 2020-1-17 10:19:29 | 显示全部楼层
fcode 发表于 2020-1-17 08:08
学了单步调试,你进步就很大了。
此后,很多问题你都可以自己排查了。

实不相瞒,看的正是你的某一个回帖的链接,你们几个大佬对我的学习还是有很大的帮助的
心怀不惧,方能翱翔于天际
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-9-19 08:35

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表