[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
