[其他行业算法] 求谐波

      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