资源描述：
《地磁要素程序》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、programdiciyaosuparameter(n=10,m=10,r=6371200)integeri,j!r=R!D磁偏角，C磁倾角,H1水平强度，z曲直强度，T总场强realx,y,z,q,g(1:10,0:10),h(1:10,0:10),v,d,c,h1,t,zh,p(0:n+1,0:m),p1(0:n+1,0:m)open(1,file='1.dat')doi=1,10doj=0,iread(1,*)g(i,j),h(i,j)enddoenddoclose(1)doq=5,175,5calllerangde(p,p1,q,n,m)dov=-
2、180,180,5x=0.0y=0.0z=0.0doi=1,10doj=0,ix=x+(g(i,j)*cosd(j*v)+h(i,j)*sind(j*v))*p1(i,j)y=y+j/sind(q)*(g(i,j)*sind(j*v)-h(i,j)*cosd(j*v))*p(i,j)z=z-(i+1)*(g(i,j)*cosd(j*v)+h(i,j)*sind(j*v))*p(i,j)enddoenddod=atan2d(y,x)h1=sqrt(x**2+y**2)c=atan2d(z,h1)t=sqrt(x**2+y**2+z**2)open(2,fil
3、e='2.dat')write(2,*)90-q,v,z,d,h1,c,tenddoenddoend!勒让德函数递推公式subroutinelerangde(p,p1,q,n,m)integernreali,j,k,zh,p(0:n+1,0:m),p1(0:n+1,0:m)p(0,0)=1p(1,0)=cosd(q)doi=2,np(i,0)=((2*i-1)*cosd(q)*p(i-1,0)-(i-1)*p(i-2,0))/ienddodoj=1,10doi=j,j+1callhan(zh,i,j)a=(2*i+1)*(1-cosd(q)**2)**0.
4、5p(i,j)=((i+j)*(i+j-1)*p(i-1,j-1)-(i-j+2)*(i-j+1)*p(i+1,j-1))/a*zhenddodok=i+1,ncallhan(zh,k,j)p(k,j)=((2*k-1)*cosd(q)*p(k-1,j)-(k-1+j)*p(k-2,j))/(k-j)*zhenddoenddodoi=1,10doj=0,ip1(i,j)=((i+1)*(2*i+1)*cosd(q)*p(i,j)-(i-j+1)*(i+1)*p(i+1,j)-i*(i-j+1)*p(i+1,j))/(2*i+1)*(1-cosd(q)**2
5、)&&*(-sind(q))enddoenddoendsubroutinesubroutinehan(zh,i1,j1)realzh,k,f,i1,j1if(j1==0)thenxm=1elsexm=2endifk=jch(i1-j1)f=jch(i1+j1)zh=(xm*k/f)**0.5contains!阶乘functionjch(e)reale,jchjch=1.0if(e==0.or.e==1)thenjch=1.0elsedoi=1,ejch=jch*ienddoendifendfunctionendsubroutine等偏线图（°）等倾线图（°
6、）水平强度等值线图（103nT)垂直强度等值线图（103nT)总场强度等值线图（103nT)
