长安大学全球地磁场IGRF报告.docVIP

长安大学全球地磁参考场IGRF作业

姓名xxx学号xxxxxx

1.实验原理

N根据以前的学习我们知道，假设假设地球为均匀磁化球体，球体半径为R，N为地理北极(其对应地磁偶极子S极性)。假设采用球极坐标系〔图六〕，以球心为坐标原点，球外任意一点P的地心距离为r，余纬度为θ〔θ=90°—φ，φ为纬度〕，经度为λ。那么在地磁场源区之外空间域坐标系〔r,θ,λ〕中，磁位u的拉普拉斯方程可以写成如下形式：

N

对上式采用别离变量法，即可解得拉普拉斯方程的一般解，从而分别获得其内源场和外源场的磁位球谐表达式。假设设外源场磁位为零，那么内源场的磁位球谐一般表达式为

①

其中为施密特准归一化的缔合勒让德函数。

②

其中为勒让德多项式，

③

为n阶m次缔合〔或伴随〕勒让德多项式。当m=0，Cm=1；当m≠0，Cm=2.

而、为内源场磁位的球谐级数系数，它与球体内任意体积元的磁荷量dm0有关。假设小体积元中心点坐标为〔r0，θ0，λ0〕，那么计算①式沿轴向的微商位，可得其北向水平分量X、东向水平分量Y、垂直分量Z如下：

其中，R为国际参考球半径，即地球平均半径，R=6371.2Km；θ为余纬；λ为地理经度；g，h为球谐系数。

通过这三个分量可以计算出其他要素值：

标量总强度F=〔X2+Y2+Z2〕1/2，水平强度H=〔X2+Y2〕1/2，磁偏角D=arctan〔Y/X〕，磁倾角I=arctan〔Z/H〕。

2.编程思路

通过以上原理，我确定了如下编程思路：

输入IGRF→对经度，余纬，以及半径进行插值→求取相应勒让德多项式→求取n阶m次缔合勒让德多项式→求取施密特归一化的缔合勒让德函数→求取X，Y，Z分量→求取其他分量。

3.程序代码

programglobal

realn,m,dam,xi,r,na,xi2,RE,x,Y,Z,cosxi,xinxi,lin,ce1,ce2,ce3,pn1,h,sum,D,I,an

!定义n,m同书上定义，dam为n次勒让德函数的连加上限,xi(余纬度),r〔半径〕,na〔经度〕为地磁三要素，RE为地球半径,XYZ为三坐标轴分量

!cosxi为cosξ，lin为临时变量,xi2为纬度,ce系列为测试所用临时变量,h为水平分量，sum为总场强,D,I分别为磁偏角，磁倾角，an为年份

real,dimension(:,:),allocatable::gn,hn,gna,hna

!各变量定义均以书上名称为准

1format(f7.2,2x,f7.2,2x,f15.5,2x,f15.5,2x,f15.5,2x,f15.5)

2format(f30.5,\tc,f30.5,\tc,f30.5)

allocate(gn(1:10,0:10))

allocate(hn(1:10,0:10))

allocate(gna(1:10,0:10))

allocate(hna(1:10,0:10))

open(1,file=input1.txt)

open(2,file=check1.txt)

open(3,file=x-axis.dat)

open(4,file=y-axis.dat)

open(5,file=z-axis.dat)

open(6,file=horizontal.dat)

open(7,file=sum.dat)

open(8,file=d.dat)

open(9,file=I.dat)

don=1,10

dom=0,n,1

read(1,*)gn(n,m),hn(n,m),gna(n,m),hna(n,m)

enddo

enddo

!输入IGRF球谐系数

don=1,10

dom=0,n,1

write(2,1)n,m,gn(n,m),hn(n,m),gna(n,m),hna(n,m)

enddo

write(2,(**********))

enddo

!以上检查输入是否正确

pause

ce1=15.0

ce2=jie(ce1)

print*,ce2

ce1=5.0

ce2=jie(ce1)

print*,ce2

pause

!以上为杂项检查

an=1991-1990

!以上规定年份

z=0

x=0

y=0

d=0

i=0

sum=0

dona=180,-180,-1

doxi=80,-80,-1

xi2=90-xi

sinxi=sind(xi2)

cosxi=cosd(xi2)

don=1,10

dom=0,n,1

gn(n,m)=gn(n,m)+an*gna(n,m)

hn(n,m)=hn(n,m)+an*hna(n,m)

z=z+(n+1)*(gn(n,m)*cosd(m*na)+hn(n,m)*sind(m*na))*pnm(n,m,cosxi)