subroutine igrf95_15 (year,rre,thet,phi,br,bt,bp,b,ier) c* c*********************************************************************** c* c*PRO MAGLIB c* c*VER 2017.02.23 - V 1.0 c* c*AUT spec. N A S A c*AUT spec. CNES - JC KOSIK - janvier 1996 c*AUT port. AKKA c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ magnetique pour la periode c*ROL comprise entre 1995 et 2015. c* c*PAR year (I) : annee fractionnaire (entre 1995. et 2015.) c* c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) c*PAR thet (I) : colatitude geocentrique (radians), la colatitude c*PAR : est calculee du pole nord au pole sud c*PAR phi (I) : longitude geocentrique (radians), calculee c*PAR : positive vers l'est de 0 a 2*PI c* c*PAR br (O) : composante radiale du champ magnetique le long du c* : meridien positive vers l'exterieur (gauss) c*PAR bt (O) : composante tangentielle du champ magnetique le long c*PAR : du meridien positive vers le sud (gauss) c*PAR bp (O) : composante azimuthale du champ magnetique, positive c*PAR : vers l'est (gauss) c* c*PAR b (O) : module du champ magnetique (gauss) c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c c*NOT Un test est effectue pour les colatitudes inferieurs a 1.d-15 c*NOT ou superieures a PI - 1.d-15 et le cosinus est remplace par 1.d0 c*NOT et le sinus par 1.d-15 ou bien le cosinus est remplace par -1.d0 c*NOT et le sinus par 1.d-15. c* c*NOT common : util, util2 c* c*INF utilise : sans objet c* c*HST version 1.0 - 2017.02.23 - Modele 2010 definitif c*HST et ajout du modele provisoire 2015 c* c*********************************************************************** c* implicit none c c --------------------------------- c*FON Declaration identificateur rcs_id c --------------------------------- c character rcs_id*100 c c -------------------------- c*FON Declaration des parametres c -------------------------- c double precision year double precision rre, thet, phi double precision br, bt, bp, b integer ier c c ---------------------------------- c*FON Declaration des variables communes c ---------------------------------- c double precision pi,dpi,rad,deg,pid,xmu,rayt c c*COM pi : constante pi (obtenue a partir de cos(-1.)) c*COM dpi : constante 2 * pi c*COM pid : constante pi / 2 c*COM rad : facteur de conversion degres ----> radians c*COM deg : facteur de conversion radians ----> degres c*COM xmu : constante de gravitation terrestre (km**3/sec**2) c*COM rayt : rayon equatorial terrestre (km) c common/util/pi,dpi,rad,deg,pid,xmu,rayt c double precision rgmt c c*COM rgmt : rayon geocentrique terrestre (km) c common/util2/rgmt c c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c integer lg(14,14,5) c*LOC lg : coefficients entre 1995 et 2015 c integer kmax c*LOC kmax : nombre de pas de calcul (= 14) c integer n,m,jj,km c*LOC n,m,jj,km : indices de boucles et de tableaux c double precision cof1,cof2,delt double precision f1 double precision shmit(14,14) double precision g(14,14) double precision const(14,14),fn(14),fm(14) double precision p(14,14),dp(14,14),sp(14),cp(14) double precision tzero(5) double precision tmold,sph,st,temp,par double precision cph,ct,aor,ar,rrmag c*LOC Variables de travail intermediaires c SAVE c c --------------------------------- c*FON Affectation identificateur rcs_id c --------------------------------- c data rcs_id /" >$Id$"/ c c -------------------------- c*FON Affectation des constantes c -------------------------- c c* coefficients pour 1995,2000,2005,2010,2015 c c Nouveaux coefficients pour les champs c dgrf2010 et igrf2015 c data lg/ c 1995 > 10,-296920,-22000,13350,9400,-2140,680,770,250,40,-30,3*0, s 53060,-17840,30700,-22670,7800,3520,670,-720,60,90,-60,3*0, s -23660,-4130,16810,12490,2900,2350,680,10,-60,30,20,3*0, s -2620,3020,-4270,7590,-4180,-1180,-1700,280,-90,-100,-40,3*0, s 2620,-2360,970,-3060,1220,-1660,-10,50,-140,80,-10,3*0, s 460,1650,-1430,-550,1070,-170,190,40,90,-80,40,3*0, s -170,720,670,-580,10,360,-930,80,60,-10,20,3*0, s -690,-250,40,240,170,-240,-60,-20,-50,100,20,3*0, s 110,-210,80,-230,150,110,-160,-40,-70,-20,50,3*0, s -200,150,120,-60,-80,80,50,-80,30,-80,10,3*0, s 10,0,40,50,-50,-10,-20,10,-20,-70,0,45*0, c 2000 > 10,-296194,-22677,13396,9323,-2188,723,790,244,50,-26,27,-22,-2, s 51861,-17282,30684,-22880,7868,3514,682,-740,66,94,-60,-17,-3, s -9, s -24816,-4580,16709,12521,2500,2223,742,0,-92,30,17,-19,2,3, s -2276,2934,-4911,7145,-4030,-1304,-1609,333,-79,-84,-31,15,9,1, s 2726,-2319,1198,-3038,1113,-1686,-59,91,-166,63,-5,-1,-2,-4, s 438,1719,-1331,-393,1063,-129,169,69,91,-89,37,1,9,13, s -174,637,651,-612,7,438,-904,73,70,-15,10,-7,-5,-4, s -646,-242,62,240,148,-254,-58,-12,-79,93,20,7,3,7, s 119,-215,85,-215,155,89,-149,-21,-70,-43,42,17,-3,-4, s -197,134,125,-62,-84,84,38,-82,48,-82,3,1,-4,3, s 17,0,40,49,-59,-12,-29,2,-22,-74,-11,12,-1,-1, s 1,13,-9,-26,9,-7,-28,-9,-12,-19,-9,40,-2,4, s -4,3,25,-26,7,3,0,0,3,-9,-4,8,-4,0, s -9,2,18,-4,-10,-1,7,3,6,3,-2,-5,-9,1, c 2005 > 100,-2955463,-233724,133630,92055,-22700,7360,7988,2480,558, > -217,295,-215,-16, > 507799,-166905,304769,-230583,79796,35441,6956,-7446,762,976, > -612,-160,-29,-88, > -259450,-51543,165776,124639,21065,20895,7674,-165,-1173,358, > 142,-188,21,30, > -19886,26972,-52472,67251,-37986,-13654,-15134,3873,-688,-694, > -235,144,89,28, > 28207,-22523,14515,-30536,10000,-16805,-1458,1230,-1811,501, > -15,-31,-38,-43, > 4272,18025,-12345,-1957,10385,-1355,1458,937,1017,-1076,306, > 29,96,118, > -2033,5475,6363,-6353,24,5094,-8636,542,936,-125,29,-79,-30, > -37, > -6114,-2257,682,2535,1093,-2632,-464,194,-1125,876,206,53,46, > 75, > 1120,-2088,983,-1971,1622,761,-1276,-6,-487,-666,377,180,-35, > -26, > -2011,1269,1267,-672,-816,810,292,-773,601,-922,-21,16,-36,35, > 219,10,446,476,-658,-101,-347,-86,-231,-793,-209,96,8,-5, > 26,144,-77,-227,90,-58,-269,-108,-158,-190,-139,399,-49,41, > -55,23,238,-263,61,40,1,2,28,-87,-34,88,-8,-10, > -76,33,172,-54,-107,-4,63,21,53,38,-22,-57,-82,-18, c 2010 (passage en coeff définitifs le 23.02.2017) > 100,-2949657,-239606,133985,91266,-23087,7278,8044,2441,550,-194, > 305,-212,-9, > 494426,-158642,302634,-232654,80897,35729,6869,-7500,821,945, > -624,-148,-21,-89, > -270854,-57573,166817,123210,16658,20026,7592,-455,-1450,345, > 89,-203,30,31, > -16040,25175,-53703,63373,-35683,-14105,-14140,4524,-559,-527, > -107,165,104,42, > 28648,-21103,16446,-30972,8940,-16317,-2283,1400,-1934,313, > -16,-51,-63,-45, > 4458,18901,-11806,-1,10104,-803,1310,1046,1161,-1238,245,54, > 95,108, > -2090,4418,6154,-6626,302,5540,-7809,164,1085,-76,-33,-79,-11, > -31, > -5780,-2120,654,2496,703,-2761,-328,492,-1405,843,213,37,52,78, > 1084,-2003,1183,-1741,1671,696,-1074,164,-354,-842,309,179, > -39,-18, > -2054,1151,1275,-714,-742,797,214,-608,701,-1008,-103,12,-37, > 38, > 273,-10,471,444,-722,-96,-395,-199,-197,-831,-280,75,21,2, > 13,167,-66,-176,85,-39,-251,-127,-211,-194,-186,375,-77,42, > -87,27,213,-249,49,59,0,13,27,-86,-23,87,4,-26, > -87,30,166,-59,-114,-7,54,10,49,44,-25,-53,-79,-26, c 2015 (création des coeff provisoires le 23.02.2017) > 10,-294420,-24451,13507,9076,-2326,700,816,242,54,-19,31,-19,0, > 47971,-15010,30129,-23523,8137,3601,677,-761,88,88,-63,-15,-2, > -9, > -28456,-6419,16767,12256,1204,1924,727,-68,-169,31,1,-23,4,4, > -1153,2449,-5384,5820,-3349,-1409,-1299,518,-32,-33,5,20,12,5, > 2833,-1887,1809,-3295,704,-1575,-289,150,-206,7,-5,-8,-8,-5, > 473,1970,-1193,160,1002,41,132,94,134,-133,18,6,9,10, > -208,332,589,-667,73,626,-709,-28,117,-1,-7,-7,1,-2, > -541,-195,57,244,34,-274,-22,68,-159,87,21,2,5,8, > 101,-183,133,-146,162,57,-91,21,-20,-91,24,17,-3,-1, > -216,108,118,-68,-69,78,10,-40,84,-105,-18,-2,-4,3, > 32,-4,46,44,-79,-6,-42,-28,-12,-87,-36,4,2,1, > -1,20,-7,-11,8,-2,-22,-14,-25,-20,-24,35,-9,5, > -11,4,19,-22,3,7,-1,3,2,-9,-1,7,0,-4, > -9,4,16,-5,-12,-1,4,-1,4,5,-3,-4,-8,-3 s/ c data kmax /14/ c data sp(1) /0.d0/ data p(1,1) /1.d0/ data cp(1) /1.d0/ data dp(1,1) /0.d0/ c data shmit(1,1) /0.d0/ c data tzero /1995.d0,2000.d0,2005.d0,2010.d0,2015.d0/ c data tmold /0.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c c ------------------------------------------- c*FON Initialisation a 0 des composantes du champ c*FON ainsi que du module c ------------------------------------------- c br = 0.0d0 bt = 0.0d0 bp = 0.0d0 b = 0.0d0 c c ------------------------------------------------- c*FON Initialisations faites au premier appel du module c ------------------------------------------------- c if (shmit(1,1) .ne. -1.d0) then c c ----------------------------------- c*FON Calcul des tableaux fn,fm et const c ----------------------------------- c shmit(1,1) = -1.d0 do 10 n = 1, 14 fn(n) = n do 20 m = 1, 14 fm(m) = m - 1 const(n,m) = dble((n - 2)**2 - (m - 1)**2) / > dble((2 * n - 3) * (2 * n - 5)) 20 continue 10 continue c c ---------------------------- c*FON Calcul du tableau shmit(n,m) c ---------------------------- c do 30 n = 2, 14 shmit(n,1) = dble(2 * n - 3) * shmit(n-1,1) > / dble(n - 1) jj = 2 do 40 m = 2, n shmit(n,m) = shmit(n,m-1) > * sqrt(dble((n - m + 1) * jj) / > (n + m - 2)) shmit(m-1,n) = shmit(n,m) jj = 1 40 continue 30 continue c c ------------------------------------------- c*FON Traitements effectues si c'est le 1er appel c ------------------------------------------- c endif c c -------------------------------------------------- c*FON Calcul du g(n,m), interpolation entre deux jeux de c*FON coefficients avec un delai de 5 ans entre eux c -------------------------------------------------- c if (year .ne. tmold) then c tmold = year c km = idint((year - 1995.d0) / 5.d0) + 1 c delt = (year - tzero(km)) / 5.d0 c do 70 n = 1, kmax do 80 m = 1, kmax cof1 = dble(lg(n,m,km))/dble(lg(1,1,km)) cof2 = dble(lg(n,m,km+1))/dble(lg(1,1,km+1)) g(n,m) = cof1 + (cof2 - cof1) * delt 80 continue 70 continue c do 77 n = 1, kmax do 78 m= 1, kmax g(n,m) = g(n,m) * shmit(n,m) 78 continue 77 continue c endif c c ---------------------------------------------------- c*FON Calcul des cosinus directeurs de la colatitude et de c*FON la longitude. On test si la colatitude est comprise c*FON entre 0 et PI c ---------------------------------------------------- c cph = cos(phi) sph = sin(phi) c if (thet .ne. 0.0d0 .and. thet .ne. acos(-1.0d0)) then c ct = cos(thet) st = sin(thet) c else c if (thet .ne. acos(-1.0d0)) then ct = 1.0d0 st = 1.d-15 elseif (thet .ne. 0.0d0) then ct = -1.0d0 st = 1.d-15 endif c endif c sp(2) = sph cp(2) = cph c c ---------------------------- c*FON Calcul des tableaux cp et sp c ---------------------------- c do 90 m = 3, kmax sp(m) = sp(2) * cp(m-1) + cp(2) * sp(m-1) cp(m) = cp(2) * cp(m-1) - sp(2) * sp(m-1) 90 continue c c -------------------------------- c*FON Calcul des trois premiers termes c -------------------------------- c rrmag = rre * rayt / rgmt aor = 1.d0 / rrmag ar = aor * aor * aor p(2,1) = ct dp(2,1) = -st p(2,2) = st dp(2,2) = ct br = -(ar + ar) * (g(2,1) * p(2,1) + p(2,2) > * (g(2,2) * cp(2) + g(1,2) * sp(2))) bt = ar * (g(2,1) * dp(2,1) + dp(2,2) > * (g(2,2) * cp(2) + g(1,2) * sp(2))) bp = ar * (g(1,2) * cp(2) - g(2,2) * sp(2)) * p(2,2) c c ------------------------------------------------ c*FON Calcul final des composantes du champ magnetique c ------------------------------------------------ c do 100 n = 3, kmax ar = aor * ar do 110 m = 1, n c if (m .eq. n) then p(n,n) = st * p(n-1,n-1) dp(n,n) = st * dp(n-1,n-1) + ct * p(n-1,n-1) else if (const(n,m) .eq. 0) then p(n,m) = ct * p(n-1,m) dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m) else p(n,m) = ct * p(n-1,m) - const(n,m) * p(n-2,m) dp(n,m) = ct * dp(n-1,m) - st * p(n-1,m) > - const(n,m) * dp(n-2,m) endif endif c par = p(n,m) * ar c if (m .eq. 1) then temp = g(n,m) else temp = g(n,m) * cp(m) + g(m-1,n) * sp(m) bp = bp - (g(n,m) * sp(m) - g(m-1,n) * cp(m)) > * fm(m) * par endif c br = br - temp * fn(n) * par bt = bt + temp * dp(n,m) * ar c 110 continue c 100 continue c bp = bp / st c c --------------------------------------- c*FON Expression du champ magnetique en gauss c --------------------------------------- c b = sqrt(bp * bp + bt * bt + br * br) c c **************** c Fin de programme c **************** c return end