subroutine chp05 (year,rre,thetr,phir,br,bt,bp,bb,ier) c* c*********************************************************************** c* c*PRO MAGLIB c* c*VER 01.05.30 - V 2.0 c*VER 03.01.06 - V 2.1 c*VER 10.09.27 - V 3.0 c*VER 17.02.22 - V 4.0 c* c*AUT spec. CNES - JC KOSIK - fevrier 2001 c*AUT port. CS-SI c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul des composantes du champ magnetique apres 2000. c* c*PAR year (I) : annee fractionnaire (>= 2005.) c* c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) c*PAR thetr (I) : colatitude geocentrique (radians) c*PAR phir (I) : longitude geocentrique (radians) c* c*PAR br (O) : composante radiale du champ magnetique (gauss) c*PAR bt (O) : composante tangentielle du champ magnetique le long c*PAR : du meridien positive vers le bas (gauss) c*PAR bp (O) : composante azimuthale du champ magnetique, positive c*PAR : vers l'est (gauss) c* c*PAR bb (O) : module du champ magnetique (gauss) c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c* c*NOT common : util, util2 c* c*INF utilise : sans objet c* c*HST version 2.0 - 01.05.30 - Enrichissement de la maglib au CDPP c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 c*HST version 3.0 - 10.09.27 - Passage aux coefficiens de Schmidt c*HST définitifs et aux variations seculaires c*HST interpolees c*HST Modif de fn1 (10->100) car précision c*HST des coeff 2005 sur 2 décimales c*HST Modif de fn2 (100->1000) car précision c*HST des coeff 2005 sur 2 décimales c*HST version 5.0 - 22.02.17 - Passage aux variations séculaires c* definitives 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,rre,thetr,phir,br,bt,bp,bb 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 acos(-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(195),lgt(195) c*LOC lg, lgt : coefficients apres l'an 2005 c integer kmax c*LOC kmax : nombre de pas de calcul (= 14) c integer l,k,jj,nm,ntot,inc c*LOC l,k,jj,nm,ntot,inc : indices de boucles et de tableaux c double precision gg(14,14),ggt(14,14),hh(14,14),hht(14,14) double precision g(14,14),h(14,14),sh(14,14),fn(14),fm(14) double precision cp(14),sp(14),p(14,14),dp(14,14),const(14,14) double precision f1, f2, tzero, tmold,t,cph,sph,ct,st double precision rrmag,ar,aor,temp,par 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 Revision septembre 2010 - coefficients lg c Revision fevrier 2017 - coefficients lgt c data lg/ > -2955463,-166905,507799,-233724,304769,-259450,165776,-51543, > 133630,-230583,-19886,124639,26972,67251, > -52472,92055,79796,28207,21065,-22523,-37986,14515,10000, > -30536,-22700,35441,4272,20895, > 18025,-13654,-12345,-16805,-1957,-1355,10385,7360,6956,-2033, > 7674,5475,-15134,6363, > -1458,-6353,1458,24,-8636,5094,7988,-7446,-6114,-165,-2257, > 3873,682,1230, > 2535,937,1093,542,-2632,194,-464,2480,762,1120,-1173,-2088, > -688,983, > -1811,-1971,1017,1622,936,761,-1125,-1276,-487,-6,558,976, > -2011,358, > 1269,-694,1267,501,-672,-1076,-816,-125,810,876,292,-666,-773, > -922, > 601,-217,-612,219,142,10,-235,446,-15,476,306,-658,29,-101, > 206,-347,377,-86,-21,-231,-209,-793,295,-160,26,-188,144,144, > -77,-31,-227,29,90,-79,-58,53,-269,180,-108,16,-158,96, > -190,399,-139,-215,-29,-55,21,23,89,238,-38,-263,96,61, > -30,40,46,1,-35,2,-36,28,8,-87,-49,-34,-8,88, > -16,-88,-76,30,33,28,172,-43,-54,118,-107,-37,-4,75, > 63,-26,21,35,53,-5,38,41,-22,-10,-57,-18,-82 > / c data lgt/ > 11612,16526,-26746,-11764,-4270,-22808,2082,-12060,710,-4142, > 7692,-2858,-3594,-7756, > -2462,-1578,2202,882,-8814,2840,4606,3862,-2120,-872,-774,576, > 372,-1738, > 1752,-902,1078,976,3912,1104,-562,-164,-174,-114,-164,-2114, > 1988,-418, > -1650,-546,-296,556,1654,892,112,-108,668,-580,274,1302,-56,340, > -78,218,-780,-756,-258,596,272,-78,118,-72,-554,170,258,400, > -246,460,288,98,298,-130,-560,404,266,340,-16,-62,-86,-26, > -236,334,16,-376,-84,-324,148,98,-26,-66,-156,-352,330,-172, > 200,46,-24,108,-106,-40,256,50,-2,-64,-122,-128,-124,10, > 14,-96,-136,-226,-164,68,-142,-76,20,24,-26,-30,46,42, > 22,-40,102,50,-10,0,38,-32,36,-2,-38,-8,-106,-42, > -8,-48,-94,6,16,-64,18,8,30,-50,-50,28,-2,-24, > 38,38,12,-2,-8,22,-2,-2,26,2,-56,22,24,-2, > 14,-2,-22,2,-6,28,-12,-4,-10,-20,-14,12,-6,6, > -18,16,-22,6,-8,14,12,2,-6,-32,8,-16,6 > / c data kmax /14/ c data f1 /100.d0/, f2 /1000.d0/ 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 tzero /2005.d0/ c data tmold /0.d0/ c data sh(1,1) /0.d0/ c c ****************** c Debut de programme c ****************** c c --------------- c*FON Initialisations c --------------- c ier = 0 c if (sh(1,1) .ne. -1.d0) then sh(1,1) = -1.d0 c c --------------------------------------------------- c*FON Calcul des coefficients de Shmidt et des constantes c --------------------------------------------------- c do 20 k = 1, kmax fn(k) = k do 10 l = 1, k fm(l) = l - 1 const(k,l) = dble((k - 2)**2 - (l - 1)**2) / > dble((2 * k - 3) * (2 * k - 5)) 10 continue 20 continue c do 40 k = 2, kmax sh(k,1) = dble(2 * k - 3) * sh(k-1,1) / dble(k - 1) jj = 2 do 30 l = 2,k sh(k,l) = sh(k,l-1) * dsqrt(dble((k - l + 1) * jj) > / (k + l - 2)) jj = 1 30 continue 40 continue c ntot = 0 do 60 k = 2, kmax c gg(k,1) = dble(lg(ntot+1)) ggt(k,1) = dble(lgt(ntot+1)) c nm = 2 * k - 1 do 50 l = 2, k inc = 2 * l - 2 gg(k,l) = dble(lg(ntot+inc)) ggt(k,l) = dble(lgt(ntot+inc)) hh(k,l) = dble(lg(ntot+inc+1)) hht(k,l) = dble(lgt(ntot+inc+1)) 50 continue c ntot = ntot + nm 60 continue c do 80 k = 2, kmax gg(k,1) = gg(k,1) * sh(k,1) / f1 ggt(k,1) = ggt(k,1) * sh(k,1) / f2 do 70 l = 2, k gg(k,l) = gg(k,l) * sh(k,l) / f1 ggt(k,l) = ggt(k,l) * sh(k,l) / f2 hh(k,l) = hh(k,l) * sh(k,l) / f1 hht(k,l) = hht(k,l) * sh(k,l) / f2 70 continue 80 continue endif c if (year .ne. tmold) then tmold = year t = year - tzero do 100 k = 2, kmax g(k,1) = gg(k,1) + t * ggt(k,1) do 90 l = 2,k g(k,l) = gg(k,l) + t * ggt(k,l) h(k,l) = hh(k,l) + t * hht(k,l) 90 continue c 100 continue c endif c cph = dcos(phir) sph = dsin(phir) c if (thetr .ne. 0.d0) then ct = dcos(thetr) st = dsin(thetr) else ct = 1.0d0 st = 1.d-15 endif sp(2) = sph cp(2) = cph c do 110 k = 3, kmax c sp(k) = sp(2) * cp(k-1) + cp(2) * sp(k-1) cp(k) = cp(2) * cp(k-1) - sp(2) * sp(k-1) c 110 continue 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 c br = -(2.d0 * ar) * (g(2,1) * p(2,1) + p(2,2) > * (g(2,2) * cp(2) + h(2,2) * sp(2))) bt = ar * (g(2,1) * dp(2,1) + dp(2,2) * (g(2,2) > * cp(2) + h(2,2) * sp(2))) bp = ar * (h(2,2) * cp(2) - g(2,2) * sp(2)) * p(2,2) c do 130 k = 3, kmax c ar = aor * ar c do 120 l = 1, k c if (k .eq. l) then p(k,k) = st * p(k-1,k-1) dp(k,k) = st * dp(k-1,k-1) + ct * p(k-1,k-1) else p(k,l) = ct * p(k-1,l) - const(k,l) * p(k-2,l) dp(k,l) = ct * dp(k-1,l) - st * p(k-1,l) - const(k,l) > * dp(k-2,l) endif c par = p(k,l) * ar c if (l .eq. 1) then temp = g(k,l) else temp = g(k,l) * cp(l) + h(k,l) * sp(l) bp = bp - (g(k,l) * sp(l) - h(k,l) * > cp(l)) * fm(l) * par endif br = br - temp * fn(k) * par bt = bt + temp * dp(k,l) * ar c 120 continue c 130 continue c bp = bp / st c c --------------------------------------- c*FON Expression du champ magnetique en gauss c --------------------------------------- c br = br / 100000.d0 bt = bt / 100000.d0 bp = bp / 100000.d0 c bb = dsqrt(br * br + bt * bt + bp * bp) c c **************** c Fin de programme c **************** c return end