!*************************************************************************************************** ! Copyright 1976, K. Pfitzer ! ! This file is part of IRBEM-LIB. ! ! IRBEM-LIB is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! IRBEM-LIB is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public License ! along with IRBEM-LIB. If not, see . ! c C----------------------------------------------------------------------- SUBROUTINE BXYZMU(x, y , z, abx, aby, abz) C C VERSION 11/01/76 C C PURPOSE C TO CALCULATE THE CONTRIBUTION TO THE EARTHS MAGNETIC FIELD BY SOURCES C EXTERNAL TO THE EARTH. NO INTERNAL FIELD IS INCLUDED IN THIS ROUTINE. C C METHOD C THE ROUTINE INCLUDES THE FIELD CONTRIBUTIONS FROM THE MAGNETOPAUSE C CURRENTS, AND CURRENTS DISTRIBUTED THROUGHOUT THE MAGNETOSPHERE C (THE TAIL AND RING CURRENTS). IT IS VALID FOR ALL TILTS OF THEEARTH'S C DIPOLE AXIS AND IS VALID DURING QUIET MAGNETIC CONDITIONS. C A GENERALIZED ORTHONORMAL LEAST-SQUARES PROGRAM WAS USED TO FIT THE C COEFFICIENTS OF A POWER SERIES (INCLUDING EXPONENTIAL TERMS) THROUGH C FOURTH ORDER IN SPACE AND THIRD ORDER IN TILT. THIS EXPANSION HAS C BEEN OPTIMIZED FOR THE NEAR EARTH REGION AND IS VALID TO 15 EARTH C RADII. OUTSIDE OF THIS REGION THE FIELD DIVERGES RAPIDLY AND A C TEMPLATE SETS THE FIELD TO ZERO. IN ORDER TO IMPROVE COMPUTATIONAL C SPEED THE FIELD IS SET TO ZERO BELOW 2 EARTH RADII. (IN THIS REGION C THE EARTH'S INTERNAL FIELD DOMINATES AND THE VARIATIONS EXPRESSED C BY THIS EXPANSION ARE NOT SUFFICIENTLY ACCURATE TO PREDICT VARIATIONS C ON THE EARTH'S SURFACE). C C THE POWER SERIES REPRESENTING THE MAGNETIC FIELD IS C BX=SUM OVER I,J,K OF ( A(I,J,K)*X**(I-1)*Y**(2*J-2)*Z**(K-1) C + B(I,J,K)*X**(I-1)*Y**(2*J-2)*Z**(K-1)*EXP(-.06*R**2)) C I GOES FROM 1 TO 5,J FROM 1 TO 3, K FROM 1 TO 5 C THE SUM OF I + 2*J + K IS LESS THAN OR EQUAL TO 9 C BY=SUM OVER I,J,K OF ( C(I,J,K)*X**(I-1)*Y**(2*J-1)*Z**(K-1) C + D(I,J,K)*X**(I-1)*Y**(2*J-1)*Z**(K-1)*EXP(-.06*R**2)) C I GOES FROM 1 TO 5, J FROM 1 TO 3, K FROM 1 TO 5 C THE SUM OF I + 2*J+1 + K IS LESS THAN OR EQUAL TO 9 C BZ=SUM OVER I,J,K OF ( E(I,J,K)*X**(I-1)*Y**(2*J-2)*Z**(K-1) C + F(I,J,K)*X**(I-1)*Y**(2*J-2)*Z**(K-1)*EXP(-.06*R**2)) C I GOES FROM 1 TO 5, J FROM 1 TO 3, K FROM 1 TO 5 C THE SUM OF I + 2*J + K IS LESS THAN OR EQUAL TO 9 C THE COEFFICIENTS A-F ARE DEPENDENT ONLY ON POSITION AND ARE C RECALCULATED EACH TIME THE TILT OF THE DIPOLE IS CHANGED. C THE COEEFICIENTS A-F ARE DETERMINED FROM THE TILT DEPENDENT C CONSTANTS AA-FF BY THE FOLLOWING EXPRESSIONS C A(I,J,K)=AA(I,J,K,1)*TILT**(K-1-(K-1)/2*2) C +AA(I,J,K,2)*TILT**(K+1-(K-1)/2*2) C B(I,J,K)=BB...... C C(I,J,K)=CC...... C D(I,J,K)=DD...... C E(I,J,K)=EE(I,J,K,1)*TILT**(K-(K)/2*2) C +EE(I,J,K,2)*TILT**(K+2-(K)/2*2) C F(I,J,K)=FF...... C C INPUT -- CALLING SEQUENCE C X,Y,Z REALS GIVING THE POSITION WHERE THE MAGNETIC FIELD C IS TO BE EVALUATED. XX(1), XX(2), XX(3) ARE RESPECTIVELY C THE X, Y, AND Z SOLAR MAGNETIC COORDINATES IN EARTH RADII. C Z IS ALONG THE EARTH'S NORTH DIPOLE AXIS, X IS PERPENDICULAR C TO Z AND IN THE PLANE CONTAINING THE Z AXIS AND THE SUN-EARTH C LINE, Y IS PERPENDICULAR TO X AND Z FORMING A RIGHT HANDED C COORDINATE SYSTEM. X IS POSITIVE IN THE SOLAR DIRECTION. C TILT IS THE TILT OF THE DIPOLE AXIS IN DEGREES. IT IS THE COMPLEMENT C OF THE ANGLE BETWEEN THE NORTH DIPOLE AXIS AND THE SOLAR C DIRECTION (PSI). TILT=90-PSI. C C OUTPUT -- CALLING SEQUENCE C ABX,ABY,ABZ THE X, Y, AND Z COMPONENTS OF THE MAGNETOSPHERIC C MAGNETIC FIELD IN GAMMA. C C CONSTANTS C AA,BB,CC,DD,EE,FF ARE REAL ARRAYS CONTAINING THE TILT DEPENDENT C COEFFICIENTS. AA(I,J,K,L) ARE STORED SUCH THAT L VARIES MOST C RAPIDLY, FOLLOWED IN ORDER BY K, J AND I. I VARIES THE SLOWEST. C THE ARRAY IS CLOSE PACKED AND ALL COEFFICIENTS THAT ARE ZERO C BECAUSE OF SYMMETRY OR BECAUSE THE CROSS TERM POWER IS TOO C LARGE ARE DELETED. C C VARIABLES C A,B,C,D,E,F THE TILT INDEPENDENT COEFFICIENTS. THEIR USE IS C DESCRIBED UNDER METHOD. C ITA A REAL ARRAY WHICH CONTAINS THE SYMMETRY OF THE TILT C DEPENDENCE FOR EACH OF THE A AND B COEFFICIENTS C ITA(1) HAS THE SYMMETRY INFORMATION FOR A(1,1,1,1) C AND A(1,1,1,2) C ITA(2) HAS THE SYMMETRY INFORMATION FOR A(1,1,2,1) C AND A(1,1,2,2) ETC. C IF ITA = 1 TILT SYMMETRY IS EVEN WITH RESPECT TO Z SYM. C IF ITA = 2 TILT SYMMETRY IS ODD WITH RESPECT TO Z SYM. C ITB SYMMETRY POINTER FOR C AND D ARRAYS C ITC SYMMETRY POINTER FOR E AND F ARRAYS C X X COMPONENT OF POSITION C Y Y COMPONENT OF POSITION C Z Z COMPONENT OF POSITION C Y2 Y**2 C Z2 Z**2 C R2 X**2 + Y**2 + Z**2 C R SQRT(R2) C I DO LOOP VARIABLE. IN THE FIELD EXPANSION LOOP IT REPRESENTS C THE POWER TO WHICH X IS CURRENTLY RAISED, I.E. X**(I-1) C J DO LOOP VARIABLE. ALSO Y**(2*J-2) C K DO LOOP VARIABLE. ALSO Z**(K-1) C XB X**(I-1) C YEXB X**(I-1)*Y**(2*J-2) C ZEYEXB X**(I-1)*Y**(2*J-2)*Z**(K-1) C IJK I + 2*J + K C II POINTS TO THE ARRAY LOCATION WHERE THE CURRENT POWER SERIES C COEFFICIENT FOR BX IS LOCATED. C JJ BY COEFFICIENT LOCATION POINTER C KK BZ COEFFICIENT LOCATION POINTER C BX,BY,BZ ARE USED TO CONSTRUCT THE MAGNETIC FIELD WITHIN THE POWER C SERIES LOOP. C EXPR EXP(-.06*R2) C TILTL HOLDS THE LAST VALUE OF THE TILT FOR WHICH THE TILT C INDEPENDENT COEFFICIENTS A-F WERE CALCULATED. C TT A REAL ARRAY HOLDING THE POWERS OF THE TILT. C TT(1)=TILT**0, TT(2)=TILT**1, ETC. C CON =0 FOR R LESS THAN 2 C =1 FOR R GREATER THAN 2.5 C GOES FROM 0 TO 1 IN THE REGION 2 TO 2.5. C C FOR MORE INFORMATION CALL OR WRITE K.A. PFITZER OR W.P. OLSON AT C MCDONNEL DOUGLAS ASTRONAUTICS CO. 5301 BOLSA AVE, HUNTINGTON C CALIF., PHONE (714) 896-3231. C IMPLICIT NONE C Input REAL*8 X, Y, Z, TILT C Output REAL*8 ABX, ABY, ABZ C Local REAL*8 AA(64), BB(64), CC(44), DD(44), EE(64), FF(64), A(32) REAL*8 B(32), C(22), D(22), E(32), F(32), TT(4), TILTL REAL*8 Y2, Z2, R2, BX, BY, BZ, CON, EXPR, XB, YEXB, ZEYEXB INTEGER*4 ITA(32), ITB(22), ITC(32), I, J, K, II, JJ, KK, IJK SAVE TILTL,TT,A,B,E,F,C,D COMMON /dip_ang/TILT DATA ITA /2,1,2,1,2,2,1,2,1,2,1,2,1,2,1,2,2,1,2,2,2,1, * 2,1,2,1,2,1,2,2,2,1/ DATA ITB /2,1,2,1,2,2,1,2,2,2,1,2,1,2,1,2,1,2,2,2,1,2/ DATA ITC /1,2,1,2,1,1,2,1,2,1,2,1,2,1,2,1,1,2,1,1,1,2, * 1,2,1,2,1,2,1,1,1,2/ DATA AA/-2.26836d-02,-1.01863d-04,3.42986d+00, *-3.12195d-04, 9.50629d-03,-2.91512d-06,-1.57317d-03, 8.62856d-08, *-4.26478d-05, 1.62924d-08,-1.27549d-04, 1.90732d-06,-1.65983d-02, * 8.46680d-09,-5.55850d-05, 1.37404d-08, 9.91815d-05, 1.59296d-08, * 4.52864d-07,-7.17669d-09, 4.98627d-05, 3.33662d-10,-5.97620d-02, * 1.60669d-05,-2.29457d-01,-1.43777d-04, 1.09403d-03,-9.15606d-07, * 1.60658d-03,-4.01198d-07,-3.15064d-06, 2.03125d-09, 4.92887d-04, *-1.80676d-07,-1.12022d-03, 5.98568d-07,-5.90009d-06, 5.16504d-09, *-1.48737d-06, 4.83477d-10,-7.44379d-04, 3.82472d-06, 7.41737d-04, *-1.31468d-05,-1.24729d-04, 1.92930d-08,-1.91764d-04,-5.30371d-08, * 1.38186d-05,-2.81594d-08, 7.46386d-06, 2.64404d-08, 2.45049d-04, *-1.81802d-07,-1.00278d-03, 1.98742d-06,-1.16425d-05, 1.17556d-08, *-2.46079d-06,-3.45831d-10, 1.02440d-05,-1.90716d-08,-4.00855d-05, * 1.25818d-07/ DATA BB/ 9.47753d-02, 1.45981d-04,-1.82933d+00, * 5.54882d-04, 5.03665d-03,-2.07698d-06, 1.10959d-01,-3.45837d-05, *-4.40075d-05, 5.06464d-07,-1.20112d-03, 3.64911d-06, 1.49849d-01, *-7.44929d-05, 2.46382d-04, 9.65870d-07,-9.54881d-04, 2.43647d-07, * 3.06520d-04, 3.07836d-07, 6.48301d-03, 1.26251d-06,-7.09548d-03, *-1.55596d-05, 3.06465d+00,-7.84893d-05, 4.95145d-03, 3.71921d-06, *-1.52002d-01, 6.81988d-06,-8.55686d-05,-9.01230d-08,-3.71458d-04, * 1.30476d-07,-1.82971d-01, 1.51390d-05,-1.45912d-04,-2.22778d-07, * 6.49278d-05,-3.72758d-08,-1.59932d-03, 8.04921d-06, 5.38012d-01, *-1.43182d-04, 1.50000d-04, 5.88020d-07,-1.59000d-02, 1.60744d-06, * 3.17837d-04, 1.78959d-07,-8.93794d-03, 6.37549d-06, 1.27887d-03, *-2.45878d-07,-1.93210d-01, 6.91233d-06,-2.80637d-04,-2.57073d-07, * 5.78343d-05, 4.52128d-10, 1.89621d-04,-4.84911d-08,-1.50058d-02, * 6.21772d-06/ DATA CC/-1.88177d-02,-1.92493d-06,-2.89064d-01, *-8.49439d-05,-4.76380d-04,-4.52998d-08, 1.61086d-03, 3.18728d-07, * 1.29159d-06, 5.52259d-10, 3.95543d-05, 5.61209d-08, 1.38287d-03, * 5.74237d-07, 1.86489d-06, 7.10175d-10, 1.45243d-07,-2.97591d-10, *-2.43029d-03,-6.70000d-07,-2.30624d-02,-6.22193d-06,-2.40815d-05, * 2.01689d-08, 1.76721d-04, 3.78689d-08, 9.88496d-06, 7.33820d-09, * 7.32126d-05, 8.43986d-08, 8.82449d-06,-6.11708d-08, 1.78881d-04, * 8.62589d-07, 3.43724d-06, 2.53783d-09,-2.04239d-07, 8.16641d-10, * 1.68075d-05, 7.62815d-09, 2.26026d-04, 3.66341d-08, 3.44637d-07, * 2.25531d-10/ DATA DD/ 2.50143d-03, 1.01200d-06, 3.23821d+00, * 1.08589d-05,-3.39199d-05,-5.27052d-07,-9.46161d-02,-1.95413d-09, *-4.23614d-06, 1.43153d-08,-2.62948d-04, 1.05138d-07,-2.15784d-01, *-2.20717d-07,-2.65687d-05, 1.26370d-08, 5.88917d-07,-1.13658d-08, * 1.64385d-03, 1.44263d-06,-1.66045d-01,-1.46096d-05, 1.22811d-04, * 3.43922d-08, 9.66760d-05,-6.32150d-07,-4.97400d-05,-2.78578d-08, * 1.77366d-02, 2.05401d-07,-1.91756d-03,-9.49392d-07,-1.99488d-01, *-2.07170d-06,-5.40443d-05, 1.59289d-08, 7.30914d-05, 3.38786d-08, *-1.59537d-04,-1.65504d-07, 1.90940d-02, 2.03238d-06, 1.01148d-04, * 5.20815d-08/ DATA EE/-2.77924d+01,-1.01457d-03, 9.21436d-02, *-8.52177d-06, 5.19106d-01, 8.28881d-05,-5.59651d-04, 1.16736d-07, *-2.11206d-03,-5.35469d-07, 4.41990d-01,-1.33679d-05,-7.18642d-04, * 6.17358d-08,-3.51990d-03,-5.29070d-07, 1.88443d-06,-6.60696d-10, *-1.34708d-03, 1.02160d-07, 1.58219d-06, 2.05040d-10, 1.18039d+00, * 1.58903d-04, 1.86944d-02,-4.46477d-06, 5.49869d-02, 4.94690d-06, *-1.18335d-04, 6.95684d-09,-2.73839d-04,-9.17883d-08, 2.79126d-02, *-1.02567d-05,-1.25427d-04, 3.07143d-08,-5.31826d-04,-2.98476d-08, *-4.89899d-05, 4.91480d-08, 3.85563d-01, 4.16966d-05, 6.74744d-04, *-2.08736d-07,-3.42654d-03,-3.13957d-06,-6.31361d-06,-2.92981d-09, *-2.63883d-03,-1.32235d-07,-6.19406d-06, 3.54334d-09, 6.65986d-03, *-5.81949d-06,-1.88809d-04, 3.62055d-08,-4.64380d-04,-2.21159d-07, *-1.77496d-04, 4.95560d-08,-3.18867d-04,-3.17697d-07,-1.05815d-05, * 2.22220d-09/ DATA FF/-5.07092d+00, 4.71960d-03,-3.79851d-03, *-3.67309d-06,-6.02439d-01, 1.08490d-04, 5.09287d-04, 5.62210d-07, * 7.05718d-02, 5.13160d-06,-2.85571d+00,-4.31728d-05, 1.03185d-03, * 1.05332d-07, 1.04106d-02, 1.60749d-05, 4.18031d-05, 3.32759d-08, * 1.20113d-01, 1.40486d-05,-3.37993d-05, 5.48340d-09, 9.10815d-02, *-4.00608d-04, 3.75393d-03,-4.69939d-07,-2.48561d-02, 1.31836d-04, *-2.67755d-04,-7.60285d-08, 3.04443d-03,-3.28956d-06, 5.82367d-01, * 5.39496d-06,-6.15261d-04, 4.05316d-08, 1.13546d-02,-4.26493d-06, *-2.72007d-02, 5.72523d-08,-2.98576d+00, 3.07325d-05, 1.51645d-03, * 1.25098d-06, 4.07213d-02, 1.05964d-05, 1.04232d-04, 1.77381d-08, * 1.92781d-01, 2.15734d-05,-1.65741d-05,-1.88683d-09, 2.44803d-01, * 1.51316d-05,-3.01157d-04, 8.47006d-08, 1.86971d-02,-6.94074d-06, * 9.13198d-03,-2.38052d-07, 1.28552d-01, 6.92595d-06,-8.36792d-05, *-6.10021d-08/ DATA TILTL /99.0D0/ c write(6,*)tilt c stop C C SET UP SOME OF THE INITIAL POSITION VARIABLES. Y2 = Y**2 Z2 = Z**2 R2 = X**2 + Y2 + Z2 C C SET MAGNETIC FIELD VARIABLES TO ZERO. BX = 0.0D0 BY = 0.0D0 BZ = 0.0D0 C C CHECK TO SEE IF POSITION IS WITHIN REGION OF VALIDITY. CON = 1.0D0 C IF DISTANCE TOO LARGE TAKE ERROR EXIT. IF (R2 .GT. 225.0D0) GO TO 50 C IF DISTANCE TOO SMALL SET FIELD TO ZERO AND EXIT. IF (R2 .LT. 4.0D0) GO TO 40 IF (R2 .LT. 6.25D0) CON = CON * (R2-4.0D0) / 2.25D0 C C IF TILT HAS NOT CHANGED, GO DIRECTLY TO FIELD CALCULATION. IF (TILTL .EQ. TILT) GO TO 6 C SET UP POWERS OF TILT. TILTL = TILT TT(1) = 1.0D0 TT(2) = TILTL TT(3) = TILTL**2 TT(4) = TILT * TT(3) C C SET UP THE X AND Z COMPONENT TILT INDEPENDENT COEFFICIENTS. DO 1 I=1,32 J = (I-1) * 2 + 1 K = ITA(I) A(I) = AA(J) * TT(K) + AA(J+1) * TT(K+2) B(I) = BB(J) * TT(K) + BB(J+1) * TT(K+2) K = ITC(I) E(I) = EE(J) * TT(K) + EE(J+1) * TT(K+2) F(I) = FF(J) * TT(K) + FF(J+1) * TT(K+2) 1 CONTINUE C C SET UP THE Y COMPONENT TILT INDEPENDENT COEFFICIENTS. DO 2 I=1,22 J = (I-1) * 2 + 1 K = ITB(I) C(I) = CC(J) * TT(K) + CC(J+1) * TT(K+2) D(I) = DD(J) * TT(K) + DD(J+1) * TT(K+2) 2 CONTINUE 6 EXPR = DEXP(-0.06D0*R2) C C INITIALIZE THE POINTERS. II=1 JJ=1 KK=1 XB=1.0D0 C C BEGIN SUM OVER X. DO 30 I=1,5 YEXB = XB C C BEGIN SUM OVER Y. DO 20 J=1,3 IF (I+2*J .GT. 8) GO TO 25 ZEYEXB = YEXB IJK = I + 2 * J + 1 K = 1 C C Z LOOP STARTS HERE. 10 BZ = BZ + (E(KK)+F(KK)*EXPR) * ZEYEXB KK = KK + 1 BX = BX + (A(II)+B(II)*EXPR) * ZEYEXB II = II + 1 IF (IJK .GT. 8) GO TO 15 BY = BY + (C(JJ)+D(JJ)*EXPR) * ZEYEXB * Y JJ = JJ + 1 ZEYEXB = ZEYEXB * Z IJK = IJK + 1 K = K + 1 IF ((IJK .LE. 9) .AND. (K .LE. 5)) GO TO 10 15 YEXB = YEXB * Y2 20 CONTINUE 25 XB = XB * X 30 CONTINUE C C SET UP THE OUTPUT ARRAY, MULTIPLY BY CON. CON IS NORMALLY ONE C BUT INSIDE OF R=2.5 IT GOES TO ZERO. INSIDE R=2 IT IS ZERO. 40 ABX = BX * CON ABY = BY * CON ABZ = BZ * CON RETURN C ERROR EXIT IF OUTSIDE OF R = 15. 50 CONTINUE c WRITE(6,*)' OUTSIDE THE VALID REGION--POWER SERIES DIVERGES c &BFIELD IS SET TO ZERO' 60 FORMAT(4H X= ,E10.3,4H Y= ,E10.3,4H Z= ,E10.3,76H 'IS OUTSIDE THE' &'VALID REGION--POWER SERIES DIVERGES BFIELD IS SET TO ZERO') GO TO 40 END