subroutine ex89ae (indval,tilt,x,y,z,bx,by,bz,ier) c* c*********************************************************************** c* c* "Copyright [c] CNES 98 - tous droits reserves" c* ********************************************** c* c*PRO MAGLIB c* c*VER 99.03.31 - V 1.0 c*VER 01.05.30 - V 2.0 c*VER 03.01.06 - V 2.1 c* c*AUT spec. Nikolai A. Tsyganenko c*AUT Institute of Physics, Leningrad University c*AUT Stary Petergof 198904 Leningrad USSR c*AUT port. CISI c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ externe d'origine magnetospherique en c*ROL fonction de l'angle de tilt et de l'indice geomagnetique Ae. c* c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ c* c*PAR tilt (I) : angle de tilt (radians) c* c*PAR x (I) : coordonnee solaire magnetique en x (rayons terrestres) c*PAR y (I) : coordonnee solaire magnetique en y (rayons terrestres) c*PAR z (I) : coordonnee solaire magnetique en z (rayons terrestres) c* c*PAR bx (O) : composante en x dans le systeme solaire magnetospherique c*PAR : (gauss) c*PAR by (O) : composante en y dans le systeme solaire magnetospherique c*PAR : (gauss) c*PAR bz (O) : composante en z dans le systeme solaire magnetospherique c*PAR : (gauss) c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c* c*NOT indval : 1 = Ae = 0 - 50 c*NOT indval : 2 = Ae = 50 - 100 c*NOT indval : 3 = Ae = 100 - 150 c*NOT indval : 4 = Ae = 150 - 250 c*NOT indval : 5 = Ae = 250 - 400 c*NOT indval : 6 = Ae >= 400 c* c*NOT Traite les composants GSM du champ magnetique produits par les c*NOT systemes extraterrestres courants dans la geomagnetosphere. c*NOT Le modele est valide jusqu'a des distances geocentriques de c*NOT 70 rayons terrestres et bases sur l'union des jeux de donnees c*NOT satellite imp-a,c,d,e,f,g,h,i,j (1966-1974) and heos-1,-2 c*NOT (1969-1974). c*NOT Reference: n.a. Tsyganenko, A magnetospheric c*NOT magnetic field model with a warped tail current sheet: planet. c*NOT space. sci., v.37, pp.5-20, 1989. c* c*INF utilise : sans objet c* c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP c*HST version 2.0 - 01.05.30 - correction de commentaires de code c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 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 integer indval double precision tilt double precision x, y, z double precision bx, by, bz integer ier c c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c integer ip c*LOC ip : valeur initiale de indval (=100) c integer i c*LOC i : indice de boucles c double precision ga1(28),ga2(28),ga3(28),ga4(28),ga5(28) double precision ga6(28),pa(28) c double precision a,a02,aa4sps,adr,adrt,adrt2,adt2r2,at double precision brzr1,brzr2,bxc,bxcf,bxdr,bxt,byc,bycf double precision byt,bzc,bzcf,bzt,cps,d,d0,dadd,dadd05,dd,ddel double precision ddop,ddopdx,ddr,ddx,dadd18,ddy,del,delx double precision dfdrx,drdycm,dwcx,dwcy,dwx,dwy,dy,dyc double precision dzsx,dzsy,ex,f3,f5,f7,f9,fcy,fdr,fk1,fk2 double precision fr,fs,fxm,fxp,fy,fym,fyp,fzm,f1,fzp double precision g,gam,gspm,h,ha02,hldxm,hlwc2m double precision hsx,htps,hxlw2m,hy,hys,psi,rc,rdy2,rdyc2 double precision ro2,rq,rqc2,rsmrt,rsprt,rsqxdl,rt double precision rx2a2,s1,sm,sp,sps,srq,srqc2,srx2a2 double precision sx,sxc,sxrc16,sxsix,t,tdr double precision w,wc,wcsm,wcsp,wt,wtfs,x2sm,xd double precision xld2,xlw2,xlwc2,xrc,xsixt,xsixtd double precision xsm,xsx,xsxc,xwcywc,xwyw,xxd,xxd2l2,y2,y4 double precision y410,yz,z2,zm,zp,zr,zs,zsm,drdy2m 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 data ga1/ -77.38d0, -10921.d0, -2.89d0, 193.9d0, > -10954.d0, 2.048d0, 28.83d0, -0.0367d0, > -0.0625d0, 0.00104d0, -1.246d0, 0.00166d0, > 0.000792d0, 20.81d0, -0.03608d0, 0.d0, > 0.d0, 0.d0, 0.d0, 25.82d0, > 7.656d0, 2.106d0, -0.30d0, 8.753d0, > 2.733d0, 13.17d0, 27.09d0, 5.184d0/ c data ga2/ -59.2d0, -13647.d0, 11.6d0, 237.2d0, > -14465.d0, 2.326d0, 36.54d0, -0.07084d0, > -0.1182d0, 0.000146d0, -1.626d0, 0.002466d0, > 0.001355d0, 23.49d0, -0.04602d0, 0.d0, > 0.d0, 0.d0, 0.d0, 23.15d0, > 7.6d0, 1.6d0, 1.457d0, 8.373d0, > 3.883d0, 14.18d0, 28.79d0, 5.97d0/ c data ga3/ -65.91d0, -15267.d0, 64.48d0, 230.6d0, > -14870.d0, 2.712d0, 43.28d0, -0.09687d0, > -0.1671d0, 0.009814d0, -1.835d0, 0.0026d0, > 0.0017d0, 22.73d0, -0.05049d0, 0.d0, > 0.d0, 0.d0, 0.d0, 20.76d0, > 6.314d0, 1.42d0, 4.141d0, 9.436d0, > 6.266d0, 15.07d0, 30.35d0, 6.852d0/ c data ga4/ -57.97d0, -16106.d0, 86.28d0, 199.7d0, > -16796.d0, 3.033d0, 46.80d0, -0.1057d0, > -0.1992d0, 0.01509d0, -1.534d0, 0.001783d0, > -0.000542d0, 22.29d0, -0.04664d0, 0.d0, > 0.d0, 0.d0, 0.d0, 19.79d0, > 5.788d0, 1.177d0, 5.007d0, 9.60d0, > 7.32d0, 15.78d0, 33.53d0, 7.401d0/ c data ga5/ -118.1d0, -16473.d0, 136.3d0, 158.7d0, > -21134.d0, 3.135d0, 49.78d0, -0.1088d0, > -0.2245d0, 0.01759d0, -1.874d0, 0.003243d0, > -0.000268d0, 21.97d0, -0.04913d0, 0.d0, > 0.d0, 0.d0, 0.d0, 18.48d0, > 6.25d0, 0.8643d0, 5.821d0, 9.552d0, > 6.051d0, 16.14d0, 31.27d0, 9.062d0/ c data ga6/ -223.7d0, -15054.d0, 219.1d0, 83.84d0, > -31140.d0, 3.777d0, 51.08d0, -0.1261d0, > -0.2393d0, 0.02507d0, -1.426d0, 0.001678d0, > 0.002039d0, 19.10d0, -0.05344d0, 0.d0, > 0.d0, 0.d0, 0.d0, 16.73d0, > 6.532d0, 0.382d0, 5.807d0, 8.099d0, > 6.726d0, 16.21d0, 25.63d0, 9.431d0/ c data del /0.01d0/ c data gam /4.d0/ c data dyc /20.d0/ c data xd /0.d0/ c data xld2 /40.d0/ c data xlw2 /170.d0/ c data a02 /25.d0/ c data rt /30.d0/ c data sxc /4.d0/ c data xlwc2 /50.d0/ c data dadd /1.d0/ c data ip /100/ c data psi /10.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c c -------------------------------------------------------------- c*FON Calcul initial des constantes a partir d'un jeu donne de c*FON parametres d'un modele specifie par le numero en option (indval) c -------------------------------------------------------------- c c Cas ou indval est egal a 100 if (indval.ne.ip) then c ip = indval c do 10 i = 1, 28 if (ip .eq. 1) pa(i) = ga1(i) if (ip .eq. 2) pa(i) = ga2(i) if (ip .eq. 3) pa(i) = ga3(i) if (ip .eq. 4) pa(i) = ga4(i) if (ip .eq. 5) pa(i) = ga5(i) if (ip .eq. 6) pa(i) = ga6(i) 10 continue c delx = pa(20) adr = pa(21) d0 = pa(22) dd = pa(23) rc = pa(24) g = pa(25) a = pa(26) dy = pa(27) sx = pa(28) ha02 = 0.5d0 * a02 rdyc2 = 1.d0 / dyc**2 hlwc2m = -0.5d0 * xlwc2 drdycm = -2.d0 * rdyc2 hldxm = -0.5d0 * xld2 ddel = 2.d0 * del rdy2 = 1.d0 / dy**2 drdy2m = -2.d0 * rdy2 hxlw2m = -0.5d0 * xlw2 dadd05 = dadd * 0.5d0 dadd18 = -18.d0 * dadd c c ------------------------------------------------------------ c*FON Les coefficients pa(16)-pa(19) sont deduits de pa(6)-pa(13) c*FON de facon a ce que le champ magnetique soit moins divergent c ----------------------------------------------------------- c pa(16) = -0.5d0 * (pa(6) / delx + pa(10)) pa(17) = -(pa(7) / delx + pa(11)) pa(18) = -(pa(8) / delx + 3.d0 * pa(12)) pa(19) = -(pa(9) / delx + pa(13)) / 3.d0 c psi = tilt sps = sin(tilt) cps = cos(tilt) htps = sps / (2.d0 * cps) gspm = -g * sps c c cas ou indval prend les autres valeurs else c if (abs(tilt - psi) .lt. 1.d-6) then c aucun traitement else psi = tilt sps = sin(tilt) cps = cos(tilt) htps = sps / (2.d0 * cps) gspm = -g * sps endif c endif c c ------------------------------------------------------------ c*FON Le traitement commence ici si indval et tilt sont identiques c ------------------------------------------------------------ c xsm = x * cps - z * sps zsm = z * cps + x * sps x2sm = xsm**2 y2 = y * y ro2 = x2sm + y2 xxd = xsm - xd xxd2l2 = 1.d0 / (xxd**2 + xld2) rsqxdl = sqrt(xxd2l2) h = 0.5d0 * (1.d0 + xxd * rsqxdl) hsx = -hldxm * xxd2l2 * rsqxdl xsixt = xsm + 16.d0 xsixtd = 1.d0 / (xsixt**2 + 36.d0) sxsix = sqrt(xsixtd) ddop = dadd05 * (1.d0 - xsixt * sxsix) ddopdx = dadd18 * xsixtd * sxsix d = d0 + del * y2 + gam * h + ddop ddx = gam * hsx + ddopdx ddy = ddel * y xrc = xsm + rc sxrc16 = sqrt(xrc**2 + 16.d0) y4 = y2 * y2 y410 = 1.d0 / (y4 + 1.d4) hy = y2 * y410 * y hys = hy * y410 * 4.d4 hy = hy * y zs = htps * (xrc - sxrc16) + gspm * hy dzsx = htps * (1.d0 - xrc / sxrc16) dzsy = gspm * hys c c ------------------------------------------------------- c*FON zs = zs(xsm,ysm) definissent le format du champ courant c ------------------------------------------------------- c xsx = xsm - sx rq = 1.d0 /(xsx**2 + xlw2) srq = sqrt(rq) fy = 1.d0 / (1.d0 + y2 * rdy2) w = 0.5d0 * (1.d0 - xsx * srq) * fy dwx = hxlw2m * rq * srq * fy dwy = drdy2m * w * y * fy zr = zsm - zs t = sqrt(zr**2 + d**2) at = a + t s1 = sqrt(at**2 + ro2) f5 = 1.d0 / s1 f7 = 1.d0 / (s1 + at) f1 = f5 * f7 f3 = f5**3 f9 = at * f3 xwyw = xsm * dwx + y * dwy fr = zr *(xsm * dzsx + y * dzsy) fs = fr - d * (xsm * ddx + y * ddy) wt = w / t wtfs = wt * fs brzr1 = wt * f1 brzr2 = wt * f3 bxt = (pa(1) * brzr1 + pa(2) * brzr2) * zr byt = bxt * y bxt = bxt * xsm bzt = pa(1) * (w * f5 + xwyw * f7 + wtfs * f1) + > pa(2) * (w * f9 + xwyw * f1 + wtfs * f3) c c --------------------------------------------------------- c*FON La contribution de la page courante de la queue centrale c*FON (bxt,byt,bzt) est trouvee. Maintenant traitons les champs c*FON peripheriques c ---------------------------------------------------------- c rx2a2 = 1.d0 / (x2sm + a02) srx2a2 = sqrt(rx2a2) fdr = 0.5d0 * (1.d0 + xsm * srx2a2) dfdrx = ha02 * rx2a2 * srx2a2 ddr = d0 + dd * fdr + ddop tdr = sqrt(zr**2 + ddr**2) adrt = adr + tdr adrt2 = adrt**2 adt2r2 = 1.d0 / (adrt2 + ro2) fk1 = adt2r2**2 * sqrt(adt2r2) fk2 = 3.d0 * adrt * fk1 / tdr bxdr = pa(5) * zr * fk2 byt = byt + bxdr * y bxt = bxt + bxdr * xsm bzt = bzt + pa(5) * ((2.d0 * adrt2 - ro2) * fk1 + fk2 > * (fr - ddr * (dd * dfdrx + ddopdx) * xsm)) c c --------------------------------------------------- c*FON Calcul du champ Chapman-Ferraro et de la moyenne de c*FON la contribution des champs alignes courants c --------------------------------------------------- c ex = exp(x / delx) z2 = z * z yz = y * z bxcf = ex * (cps * pa(6) * z + sps * > (pa(7) + pa(8) * y2 + pa(9) * z2)) bycf = ex * (cps * pa(10) * yz + sps * y * > (pa(11) + pa(12) * y2 + pa(13) * z2)) bzcf = ex * (cps * (pa(14) + pa(15) * y2 + pa(16) * z2) > + sps * z * (pa(17) + pa(18) * y2 + pa(19) * z2)) c c ----------------------------------------------------------- c*FON La queue magnetique renvoie les composants du champ courant c*FON (bxc,byc,bzc) c ----------------------------------------------------------- c fcy = 1.d0 / (1.d0 + y2 * rdyc2) xsxc = x - sxc rqc2 = 1.d0 / (xsxc**2 + xlwc2) srqc2 = sqrt(rqc2) wc = 0.5d0 * (1.d0 - xsxc * srqc2) * fcy dwcx = hlwc2m * rqc2 * srqc2 * fcy dwcy = drdycm * wc * y * fcy xwcywc = x * dwcx + y * dwcy ro2 = y2 + x**2 zp = z + rt zm = z - rt sp = sqrt(zp**2 + ro2) sm = sqrt(zm**2 + ro2) wcsp = wc / sp wcsm = wc / sm rsprt = 1.d0 / (sp + zp) rsmrt = 1.d0 / (sm - zm) fxp = wcsp * rsprt fxm = - wcsm * rsmrt fyp = fxp * y fym = fxm * y fxp = fxp * x fxm = fxm * x fzp = wcsp + xwcywc *rsprt fzm = wcsm + xwcywc *rsmrt aa4sps = pa(4) * sps bxc = pa(3) * (fxp + fxm) + (fxp - fxm) * aa4sps byc = pa(3) * (fyp + fym) + (fyp - fym) * aa4sps bzc = pa(3) * (fzp + fzm) + (fzp - fzm) * aa4sps c c --------------------------------------------------- c*FON Somme des champs. Les composants centraux du champ c*FON courant sont transformes en coordonnees GSM c --------------------------------------------------- c bx = bxc + bxt * cps + bzt * sps + bxcf by = byc + byt + bycf bz = bzc + bzt * cps - bxt * sps + bzcf c c **************** c Fin de programme c **************** c return end