c c main program to test IGRF.FOR c real*8 xlat, xlong, elev, year, xn, xe, xz, xt, decl, incl, r2d c c get location and time to compute potential c r2d = 57.29577951d00 write (*,1276) 1276 format (' Enter geographic latitude (degrees N): ') read (*,*) xlat 1277 format (f10.5) write (*,1278) 1278 format (' Enter geographic longitude (degrees E): ') read (*,*) xlong write (*,1280) 1280 format (' Enter elevation above sea level (meters): ') read (*,*) elev write (*,1281) 1281 format (' Enter decimal year for which field is required: ') read (*,*) year call igrf(xlat, xlong, elev, year, xn, xe, xz) xt = dsqrt(xn**2 + xe**2 + xz**2) decl = datan2(xe,xn)*r2d incl = datan2(xz, dsqrt(xe*xe + xn*xn))*r2d write (*,1938) xn, xe, xz, xt, decl, incl write (*,*) xlat, xlong, year, elev 1938 format (' xn, xe, xz, xt, decl, incl = ', 4(1x,f7.0),2(1x,f7.1)) stop end c c IGRF.FOR computes International Geomagnetic Reference Field for c any time between 1975.0 to 1990.0, using the spherical c harmonic expansion of Earth's magnetic field for 1975.0, c 1980.0, 1985.0, and the secular variation between 1985.0 c and 1990.0. c c REFERENCES IAGA Division I, Working Group 1, 1986, International c Geomagnetic Reference Field Revision 1985. EOS, June 17, c 1986, p. 523-524. c Peddie, Norman W., 1982, International Geomagnetic Refer- c ence Field: the Third Generation. J. Geomag. Geoelectr., c Vol. 34, p. 309-326. c EOS Transactions, American Geophysical Union, vol. 69, c no. 17, April 26, 1988, pages 557-558. c Chapman, S. and J. Bartels, 1940, Geomagnetism. Oxford c University Press, New York, p. 611-612. c c NOTE: INPUT and OUTPUT should be DOUBLE PRECISION (REAL*8) c c INPUT xlat: latitude of point of interest - decimal degrees c xlong: longitude " " " " - " " c elev: elevation " " " " - meters a.m.s.l. c year: time desired - decimal year, 1985.0 to 1990.0 c c OUTPUT xn: north component of magnetic field, nanoteslas (gammas) c xe: east " " " " " " c xz: vertical " " " " " " c c TESTS: Results given here are from NEIS magnetic service. Values from this c program, as of 4/19/89, differ by more than 1 nT in any component. c c D I H X Y Z F c deg min deg min nT nT nT nT nT c ---------- ---------- ------- ------- ------- ------- ------- c Latitude: 35.0 N Date: 1/1/75 Longitude: 90.0 W Elevation: 0.0 c 3 41.0 65 31.6 22647 22600 1455 49757 54668 c Latitude: 35.0 N Date: 1/1/75 Longitude: 90.0 W Elevation: 1000.0 m c 3 41.0 65 31.6 22635 22589 1453 49731 54640 c Latitude: 35.0 N Date: 1/1/78 Longitude: 90.0 W Elevation: 0.0 c 3 13.6 65 24.4 22620 22584 1273 49420 54351 c Latitude: 35.0 N Date: 1/1/80 Longitude: 90.0 W Elevation: 0.0 c 2 55.4 65 19.5 22602 22573 1152 49196 54140 c Latitude: 35.0 N Date: 1/1/83 Longitude: 90.0 W Elevation: 0.0 c 2 31.5 65 14.7 22568 22546 994 48942 53895 c Latitude: 35.0 N Date: 1/1/85 Longitude: 90.0 W Elevation: 0.0 c 2 15.5 65 11.5 22545 22528 888 48773 53732 c Latitude: 35.0 N Date: 1/1/88 Longitude: 90.0 W Elevation: 0.0 c 1 51.3 65 7.9 22504 22492 728 48549 53510 c Latitude: 35.0 N Date: 1/1/90 Longitude: 90.0 W Elevation: 0.0 c 1 35.1 65 5.4 22476 22468 621 48399 53363 c Latitude: 35.0 N Date: 1/1/93 Longitude: 90.0 W Elevation: 0.0 c 1 10.8 65 1.6 22437 22432 461 48174 53142 c Latitude: 35.0 S Date: 1/1/88 Longitude: 90 E Elevation: 100.0 m c -27 4.0 -68 40.9 19760 17596 -8991 -50634 54353 c subroutine igrf(xlt, xlg, el, year, xn, xe, xz) implicit real*8 (a-h, o-z) dimension g75(11,11), h75(11,11), g80(11,11), h80(11,11), $ g85(11,11), h85(11,11), g90(11,11), h90(11,11), $ p(13,13), v(7), vs(7) data g75 /0.,-30100.,-1902.,1276.,946.,-218.,45.,71.,14.,7.,-3., $0.,-2013.,3010.,-2144.,791.,356.,66.,-56.,6.,10.,-3.,2*0.,1632., $1260.,438.,264.,28.,1.,-1.,2.,2.,3*0.,830.,-405.,-59.,-198.,16., $-12.,-12.,-5.,4*0.,216.,-159.,1.,-14.,-8.,10.,-2.,5*0.,-49.,6., $0.,4.,-1.,5.,6*0.,-111.,12.,0.,-1.,4.,7*0.,-5.,10.,4.,1.,8*0.,1., $1.,10*0.,-2.,3.,10*0.,-1./ data h75 /12*0.,5675.,-2067.,-333.,191.,31.,-13.,-77.,6.,-21.,1., $2*0.,-68.,262.,-265.,148.,99.,-26.,-16.,16.,1.,3*0.,-223.,39., $-152.,75.,-5.,4.,7.,3.,4*0.,-288.,-83.,-41.,10.,-19.,-4.,4.,5*0., $88.,-4.,22.,6.,-5.,-4.,6*0.,11.,-23.,18.,10.,-1.,7*0.,-12.,-10., $11.,-1.,8*0.,-17.,-3.,3.,9*0.,1.,1.,10*0.,-5./ data g80 /0.,-29992.,-1997.,1281.,938.,-218.,48.,72.,18.,5.,-4., $0.,-1956.,3027.,-2180.,782.,357.,66.,-59.,6.,10.,-4.,2*0.,1663., $1251.,398.,261.,42.,2.,0.,1.,2.,3*0.,833.,-419.,-74.,-192.,21., $-11.,-12.,-5.,4*0.,199.,-162.,4.,-12.,-7.,9.,-2.,5*0.,-48.,14., $1.,4.,-3.,5.,6*0.,-108.,11.,3.,-1.,3.,7*0.,-2.,6.,7.,1.,8*0.,-1., $2.,2.,9*0.,-5.,3.,11*0./ data h80 /12*0.,5604.,-2129.,-336.,212.,46.,-15.,-82.,7.,-21.,1., $2*0.,-200.,271.,-257.,150.,93.,-27.,-18.,16.,4*0.,-252.,53.,-151., $71.,-5.,4.,9.,3.,4*0.,-297.,-78.,-43.,16.,-22.,-5.,6.,5*0.,92., $-2.,18.,9.,-6.,-4.,6*0.,17.,-23.,16.,9.,8*0.,-10.,-13.,10.,-1., $8*0.,-15.,-6.,4.,9*0.,2.,11*0.,-6./ data g85 /0.,-29877.,-2073.,1300.,937.,-215.,52.,75.,21.,5.,-4., $0.,-1903.,3045.,-2208.,780.,356.,65.,-61.,6.,10.,-4.,2*0.,1691., $1244.,363.,253.,50.,2.,0.,1.,2.,3*0.,835.,-426.,-94.,-186.,24., $-11.,-12.,-5.,4*0.,169.,-161.,4.,-6.,-9.,9.,-2.,5*0.,-48.,17.,4., $2.,-3.,5.,6*0.,-102.,9.,4.,-1.,3.,8*0.,4.,7.,1.,8*0.,-6.,2.,2., $9*0.,-5.,3.,11*0./ data h85 /12*0.,5497.,-2191.,-312.,233.,47.,-16.,-82.,7.,-21.,1., $2*0.,-309.,284.,-250.,148.,90.,-26.,-21.,16.,4*0.,-296.,68.,-155., $69.,-1.,5.,9.,3.,4*0.,-298.,-75.,-50.,23.,-25.,-5.,6.,5*0.,95., $-4.,17.,11.,-6.,-4.,6*0.,20.,-21.,12.,9.,8*0.,-6.,-16.,10.,-1., $8*0.,-10.,-6.,4.,9*0.,2.,11*0.,-6./ data g90 /0.,-29761.,-2141.5,1325.5,937.5,-208.5,59.,76.,24.5, $5.,-4.,0.,-1853.,3062.,-2231.,777.,356.5,63.5,-64.,6.,10.,-4., $2*0.,1726.,1241.,324.,245.5,58.5,-.5,1.5,1.,2.,3*0.,835.5,-433., $-110.,-183.,28.,-9.,-12.,-5.,4*0.,135.,-160.5,4.,-1.,-10.5,9.,-2., $5*0.,-48.5,21.5,6.,.5,-3.,5.,6*0.,-96.,6.5,4.5,-1.,3.,7*0.,-.5, $1.5,7.,1.,8*0.,-10.,2.,2.,9*0.,-5.,3.,11*0./ data h90 /12*0.,5374.5,-2248.5,-285.5,252.,47.5,-18.,-81.,7.5, $-21.,1.,2*0.,-410.,295.5,-239.,147.,84.5,-21.,-26.,16.,4*0.,-350., $80.5,-155.5,65.,4.5,5.5,9.,3.,4*0.,-293.5,-72.,-61.5,32.5,-29., $-5.,6.,5*0.,95.,-6.5,18.5,12.,-6.,-4.,6*0.,19.5,-20.,8.,9.,8*0., $-1.5,-16.5,10.,-1.,8*0.,-3.5,-6.,4.,9*0.,2.,11*0.,-6./ c c constants (d2r: degrees to radians) c d2r = 1.745329252d-02 pi = 3.141592654 a = 6371.2d05 nl = 10 c c convert input to radians, cm, years since 1985.0 c xlat = xlt*d2r xlong = xlg*d2r elev = el*100.0 c c increments for theta, lamda and r c dtheta = 0.01*d2r dlamda = 0.01*d2r dr = 10.0*100.0 c c compute coordinates in spherical (NOT ellipsoidal) coord. system c call oblate(xlat, elev, dcosth, r1) theta1 = dacos(dcosth) c c compute potential at 7 points (point of interest, plus a point N, S, c E, W, Up and Down from point of interest, to get gradients in those c directions c do 99 ii = 1, 7 c c N of point c if (ii .eq. 1) then xln = xlat + dtheta call oblate(xln, elev, dcosth, r1) call legendre(nl,dcosth,p) thetan = dacos(dcosth) xlamda = xlong r = r1 go to 23 endif c c S of point c if (ii .eq. 2) then xls = xlat - dtheta call oblate(xls, elev, dcosth, r1) call legendre(nl,dcosth,p) thetas = dacos(dcosth) xlamda = xlong r = r1 go to 23 endif c c E of point c if (ii .eq. 3) then xl = xlat call oblate(xl, elev, dcosth, r1) call legendre(nl,dcosth,p) xlamda = xlong + dlamda r = r1 go to 23 endif c c W of point c if (ii .eq. 4) then xlamda = xlong - dlamda r = r1 go to 23 endif c c Up from point c if (ii .eq. 5) then xlamda = xlong elu = elev + dr call oblate(xl, elu, dcosth, r1) call legendre(nl,dcosth,p) r = r1 ru = r1 go to 23 endif c c Down from point c if (ii .eq. 6) then eld = elev - dr call oblate(xl, eld, dcosth, r1) call legendre(nl,dcosth,p) r = r1 rd = r1 go to 23 endif c c At point of interest (not actually used for vector field calculation) c if (ii .eq. 7) then el0 = elev call oblate(xl, el0, dcosth, r1) call legendre(nl,dcosth,p) r = r1 go to 23 endif c c begin summing terms in the expansion for field and secular change c 23 if (year .lt. 1980.0) then yearb = 1975.0 call summer (a,r,xlamda,g75,g80,h75,h80,v(ii),vs(ii),p,nl) go to 99 end if if (year .lt. 1985.0) then yearb = 1980.0 call summer (a,r,xlamda,g80,g85,h80,h85,v(ii),vs(ii),p,nl) go to 99 end if c or else... yearb = 1985.0 call summer (a,r,xlamda,g85,g90,h85,h90,v(ii),vs(ii),p,nl) 99 continue c c field components at early time (nT) c xn1 = (v(1)-v(2))/(thetan-thetas)/r1 xe1 = (v(4)-v(3))/2.0d0/dlamda/r1/dsin(theta1) xz1 = (v(5)-v(6))/(ru-rd) c c field components at late time (nT) c xn2 = (vs(1)-vs(2))/(thetan-thetas)/r1 xe2 = (vs(4)-vs(3))/2.0d0/dlamda/r1/dsin(theta1) xz2 = (vs(5)-vs(6))/(ru-rd) c c field at time of interest c xn = xn1 + (xn2-xn1)*(year-yearb)/5.0d0 xe = xe1 + (xe2-xe1)*(year-yearb)/5.0d0 xz = xz1 + (xz2-xz1)*(year-yearb)/5.0d0 return end c c Subroutine Summer performs summation part of spherical harmonic eq. c subroutine summer (a,r,xlamda,g1,g2,h1,h2,v1,v2,p,nl) implicit real*8 (a-h, o-z) dimension g1(11,11), g2(11,11), h1(11,11), h2(11,11), p(13,13) 23 sum = 0.0 sums = 0.0 do 20 n = 1, nl np1 = n + 1 arg0 = (a/r)**(n+1) do 20 m = 0, n mp1 = m + 1 c c terms for early time field components c arg1 = g1(np1,mp1)*dcos(m*xlamda) arg2 = 0.0 if (m .eq. 0) goto 30 arg2 = h1(np1,mp1)*dsin(m*xlamda) 30 sum = sum + arg0*(arg1 + arg2)*p(np1,mp1) c c terms for late time field components c args1 = g2(np1,mp1)*dcos(m*xlamda) args2 = 0.0 if (m .eq. 0) goto 35 args2 = h2(np1,mp1)*dsin(m*xlamda) 35 sums = sums + arg0*(args1 + args2)*p(np1,mp1) 20 continue v1 = a*sum v2 = a*sums return end c c legendre(lmax,x,p) c c computes values of all associated legendre polynomials up to c degree and order lmax, with argument x. values are returned in c array p(lmax+1,lmax+1). lmax must not exceed 12. Because FORTRAN c does not allow 0 as an array index, Po,o is stored as P(1,1), etc. c subroutine legendre(lmax,x,p) implicit real*8 (a-h, o-z) dimension a(13), c(13,13), p(13,13) data a /1., 1., 2., 2., 8., 8., 16., 16., 128., 128., 256., 256., $ 1024. / data c /1., 12*0., 0., 1., 11*0., -1., 0., 3., 10*0., 0., -3., 0., 1 5., 9*0., 3., 0., -30., 0., 35., 8*0., 0., 15., 0., -70., 2 0., 63., 7*0., -5., 0., 105., 0., -315., 0., 231., 3 6*0., 0., -35., 0., 315., 0., -693., 0., 429., 5*0., 4 35., 0., -1260., 0., 6930., 0., -12012., 0., 6435., 4*0., 5 0., 315., 0., -4620., 0., 18018., 0., -25740., 0., 12155., 6 3*0., -63., 0., 3465., 0., -30030., 0., 90090., 0., 7 -109395., 0., 46189., 2*0., 0., -693., 0., 15015., 0., 8 -90090., 0., 218790., 0., -230945., 0., 88179., 0., 231., 9 0., -18018., 0., 225225., 0., -1021020., 0., 2078505., 0., 1 -1939938., 0., 676039./ c c compute polynomial for l=0, m=0: (NOTE: since FORTRAN does not allow 0 c as an array index, Pn,m must be represented as P(n+1,m+1) c p(1,1) = 1.0 c c for each higher value of l: c do 10 l = 1, lmax lp1 = l + 1 sum = 0.0 c c compute m=0 term: c do 20 n = 0, l np1 = n + 1 sum = sum + c(np1,lp1)*x**n 20 continue p(lp1,1) = sum/a(l+1) c c compute higher terms with recursion relation c do 30 m = 0, l - 1 mp1 = m + 1 p(lp1,mp1+1) = -((l-m)*x*p(lp1,mp1)-(l+m)*p(lp1-1,mp1)) $ /dsqrt(1.0-x**2) 30 continue 10 continue c c normalize terms (only for m>0) c do 40 l = 1, lmax lp1 = l + 1 do 40 m = 1, l mp1 = m + 1 p(lp1,mp1) = p(lp1,mp1)*xnorm(l,m) 40 continue return end c c function xnorm(l,m) normalizes legendre functions to Schmidt c quasi-normalized form (Chapman and Bartels, 1940, p. 610-611) c function xnorm(l,m) implicit real*8 (a-h, o-z) xnorm = dsqrt(2.0*fac(l - m)/fac(l + m)) return end c c function fac(i) finds i! (as a real*8, to prevent overflow) c function fac(i) implicit real*8 (a-h, o-z) fac = 1.0 if ((i .eq. 0) .or. (i .eq. 1)) return do 10 j = 2,i fac = fac*j 10 continue return end c c subroutine oblate computes spherical coordinates given geographic c coordinates on an oblate spheroid - uses IAU (1966) spheroid c subroutine oblate(phi, h, dcosth, r) implicit real*8 (a-h, o-z) a = 6378.16d05 b = 6356.78d05 c c compute f c arg1 = dsqrt(a*a - (a*a - b*b)*(dsin(phi))**2) top = (h*arg1 + a*a)**2 bot = (h*arg1 + b*b)**2 f = top/bot c c compute r2 (r-squared) and r c arg2 = a**4 - (a**4 - b**4)*(dsin(phi))**2 r2 = h*h + 2.0*h*arg1 + arg2/(arg1*arg1) r = dsqrt(r2) c c compute costh (cos(theta)) c dcosth = dsin(phi)/dsqrt(f*(dcos(phi))**2 + (dsin(phi))**2) return end