地磁要素程序.doc

program diciyaosuparameter(n=10,m=10,r=6371200)integer i,j!r=R!D磁偏角，C磁倾角,H1水平强度，z曲直强度，T总场强real x,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)do i=1,10do j=0,iread(1,*)g(i,j),h(i,j)end doenddoclose(1)do q=5,175,5call lerangde(p,p1,q,n,m)do v=-180,180,5x=0.0y=0.0z=0.0do i=1,10do j=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,file=2.dat)write(2,*)90-q,v,z,d,h1,c,t end doend doend !勒让德函数递推公式subroutine lerangde(p,p1,q,n,m)integer nreal i,j,k,zh,p(0:n+1,0:m),p1(0:n+1,0:m)p(0,0)=1p(1,0)=cosd(q)do i=2,np(i,0)=(2*i-1)*cosd(q)*p(i-1,0)-(i-1)*p(i-2,0)/i end dodo j=1,10do i=j,j+1call han(zh,i,j)a=(2*i+1)*(1-cosd(q)*2)*0.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*zh end dodo k=i+1,ncall han(zh,k,j)p(k,j)=(2*k-1)*cosd(q)*p(k-1,j)-(k-1+j)*p(k-2,j)/(k-j)*zh end doend dodo i=1,10do j=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)&*(-sind(q)end do end doend subroutinesubroutine han(zh,i1,j1)real zh,k,f,i1,j1if(j1=0)thenxm=1elsexm=2end ifk=jch(i1-j1)f=jch(i1+j1)zh=(xm*k/f)*0.5contains!阶乘function jch(e)real e,jchjch=1.0if(e=0.or.e=1)thenjch=1.0elsedo i=1,ejch=jch*iend doend ifend functio