subroutine ex89kp (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 1989 04 Leningrad USSR c*AUT port. CISI c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ interne d'origine magnetospherique en c*ROL fonction de l'angle de tilt et de l'indice geomagnetique Kp. 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 indval : 1 = Kp = 0 , 0+ c*NOT indval : 2 = Kp = 1- , 1 , 1+ c*NOT indval : 3 = Kp = 2- , 2 , 2+ c*NOT indval : 4 = Kp = 3- , 3 , 3+ c*NOT indval : 5 = Kp = 4- , 4 , 4+ c*NOT indval : 6 = Kp > 5- c* c*NOT ier : sans objet c* c*NOT Traite les composants GSM du champ magnetique produits par les c*NOT systemes extraterrestres courantsdans 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* 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 i c*LOC i : indice de boucles c integer ip c*LOC ip : valeur initiale de indval (= 100) c double precision ga1(28),ga2(28),ga3(28),ga4(28),ga5(28) double precision ga6(28),pa(28) double precision a,a02,aa4sps,adr,adrt,adrt2,adt2r2,at,brzr1 double precision brzr2,bxc,bxcf,bxdr,bxt,byc,bycf,byt double precision bzc,bzcf,bzt,cps,d,d0,dadd,dadd05,dadd18 double precision dd,ddel,ddop,ddopdx,ddr,ddx,ddy,del,delx double precision dfdrx,drdy2m,drdycm,dwcx,dwcy,dwx,dwy double precision dy,dyc,dzsx,dzsy,ex,f1,f3,f5,f7,f9,fcy,fdr double precision fk1,fk2,fr,fs,fxm,fxp,fy,fym,fyp,fzm,fzp double precision g,gam,gspm,h,ha02,hldxm,hlwc2m,hsx,htps double precision hxlw2m,hy,hys,psi,rc,rdy2,rdyc2,ro2,rq,rqc2 double precision rsmrt,rsprt,rsqxdl,rt,rx2a2,s1,sm,sp,sps double precision srq,srqc2,srx2a2,sx,sxc,sxrc16,sxsix double precision t,tdr,w,wc,wcsm,wcsp,wt,wtfs,x2sm,xd double precision xld2,xlw2,xlwc2,xrc,xsixt,xsixtd,xsm double precision xsx,xsxc,xwcywc,xwyw,xxd, xxd2l2 double precision y2,y4,y410,yz,z2,zm,zp,zr,zs,zsm 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/ -98.72d0, -10014.d0, 15.03d0, 76.62d0, > -10237.d0, 1.813d0, 31.1d0, -0.07464d0, > -0.07764d0, 0.003303d0, -1.129d0, 0.001663d0, > 0.000988d0, 18.21d0, -0.03018d0, 0.d0, > 0.d0, 0.d0, 0.d0, 24.74d0, > 8.16d0, 2.08d0, -0.88d0, 9.08d0, > 3.84d0, 13.55d0, 26.94d0, 5.75d0/ c data ga2/-35.65d0, -12800.d0, 14.37d0, 124.5d0, > -13543.d0, 2.316d0, 35.64d0, -0.0741d0, > -0.1081d0, 0.003924d0, -1.451d0, 0.00202d0, > 0.00111d0, 21.37d0, -0.04567d0, 0.d0, > 0.d0, 0.d0, 0.d0, 22.33d0, > 8.12d0, 1.664d0, 0.932d0, 9.24d0, > 2.43d0, 13.81d0, 28.83d0, 6.05d0/ c data ga3/ -77.45d0, -14588.d0, 64.85d0, 123.9d0, > -16229.d0, 2.641d0, 42.46d0, -0.07611d0, > -0.1579d0, 0.004078d0, -1.391d0, 0.00153d0, > 0.000727d0, 21.86d0, -0.04199d0, 0.d0, > 0.d0, 0.d0, 0.d0, 20.9d0, > 6.28d0, 1.54d0, 4.18d0, 9.61d0, > 6.59d0, 15.08d0, 30.57d0, 7.43d0/ c data ga4/ -70.12d0, -16125.d0, 90.71d0, 38.08d0, > -19630.d0, 3.181d0, 47.5d0, -0.1327d0, > -0.1864d0, 0.01382d0, -1.488d0, 0.002962d0, > 0.000897d0, 22.74d0, -0.04095d0, 0.d0, > 0.d0, 0.d0, 0.d0, 18.64d0, > 6.27d0, 0.935d0, 5.39d0, 8.57d0, > 5.94d0, 15.63d0, 31.47d0, 8.10d0/ c data ga5/ -162.5d0, -15806.d0, 160.6d0, 5.888d0, > -27534.d0, 3.607d0, 51.1d0, -0.1006d0, > -0.1927d0, 0.03353d0, -1.392d0, 0.001594d0, > 0.002439d0, 22.41d0, -0.04925d0, 0.d0, > 0.d0, 0.d0, 0.d0, 18.31d0, > 6.2d0, 0.768d0, 5.07d0, 10.06d0, > 6.67d0, 16.1d0, 30.04d0, 8.26d0/ c data ga6/ -128.4d0, -16184.d0, 149.1d0, 215.5d0, > -36435.d0, 4.09d0, 49.09d0, -0.0231d0, > -0.1359d0, 0.01989d0, -2.298d0, 0.004911d0, > 0.003421d0, 21.79d0, -0.05447d0, 0.d0, > 0.d0, 0.d0, 0.d0, 19.48d0, > 5.83d0, 0.332d0, 6.47d0, 10.47d0, > 9.08d0, 15.85d0, 25.27d0, 7.98d0/ 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 Coefficients pa(16)-pa(19) sont trouves dans pa(6)-pa(13) c*FON afin 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 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 la moyenne de la c*FON 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