subroutine tsyg87 (indval,ps,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 externe d'origine magnetospherique en c*ROL fonction de l'angle de tilt (ps) et de l'indice c*ROL geomagnetique (indval). c* c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ c* c*PAR ps (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 solaire magnetospherique x du champ c*PAR : externe (nanoteslas) c*PAR by (O) : composante solaire magnetospherique y du champ c*PAR : externee (nanoteslas) c*PAR bz (O) : composant solaire magnetospherique z du champ c*PAR : externe (nanoteslas) 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 common : util 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 ps double precision x, y, z double precision bx, by, bz 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 c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c integer ip c*LOC ip : valeur initiale de indval (= 100) c integer j c*LOC j : indice de boucles c double precision tga(6,32),pa(32) c double precision c1,rrc2,xn,rh,dy,delx,b0,b1,xn2,xn22,xn21,psi double precision sps,cps,zs,rps,zp,rt,fy,s2m,g2m,fc,xnx2,xnx,xc1 double precision x1,xc2,x2,xc22,xr2,xc12,b,b20,bp,b2p,bm,b2m double precision xa1,xap1,xam1,xa2,xap2,xam2,xna,xnap,xnam,f,fp double precision fm,s0,s0p,dstr,b2,xnr,adln,zm,xln1,xlnp1,xlnm1 double precision s1p,s2p,g1p,g2,g2p,ex1,ex2,y2,z2,s0m,xln2,xlnp2 double precision yz,zsm,rr2,xsm,zn,bzsm,by1,brsm,bxsm,xlnm2 double precision aln,s1,s2,s1m,g1,g1m 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 (tga(1,j),j = 1,32)/ > -0.09673d0, -10.63d0, 1.210d0, 34.57d0, > -0.04502d0, -0.06553, -0.02952d0, 0.3852d0, > -0.03665d0, -2.084d0, 0.001795d0, 0.00638d0, > -23.49d0, 0.06082d0, 0.01642d0, -0.02137d0, > 32.21d0, -0.04373d0, -0.02311d0, -0.2832d0, > -0.002303d0, -0.000631d0, -6.397d0, -967.d0, > -8650.d0, -20.55d0, 5.18d0, -2.796d0, > 2.715d0, 13.58d0, 8.038d0, 29.21d0/ c data (tga(2,j),j = 1,32)/ > -0.4850d0, -12.84d0, 1.856d0, 40.06d0, > -0.0294d0, -0.09071d0, -0.02993d0, 0.5465d0, > -0.04928d0, -2.453d0, 0.001587d0, 0.007402d0, > -29.41d0, 0.08101d0, 0.02322d0, -0.1091d0, > 40.75d0, -0.07995d0, -0.03859d0, -0.2755d0, > -0.002759d0, -0.000408d0, -6.189d0, -957.8d0, > -7246.d0, -25.51d0, 5.207d0, -4.184d0, > 2.641d0, 16.56d0, 7.795d0, 29.36d0/ c data (tga(3,j),j = 1,32)/ > -1.132d0, -18.05d0, 2.625d0, 48.55d0, > -0.004868d0, -0.1087d0, -0.03824d0, 0.8514d0, > -0.0522d0, -2.881d0, -0.000295d0, 0.009055d0, > -29.48d0, 0.06394d0, 0.03864d0, -0.2288d0, > 41.77d0, -0.05849d0, -0.06443d0, -0.4683d0, > 0.001222d0, -0.000519d0, -3.696d0, -991.1d0, > -6955.d0, -31.43d0, 4.878d0, -3.151d0, > 3.277d0, 19.19d0, 7.248d0, 28.99d0/ c data (tga(4,j),j = 1,32)/ > -1.003d0, -16.98d0, 3.140d0, 52.81d0, > -0.08625d0, -0.1478d0, -0.03501d0, 0.5500d0, > -0.07778d0, -2.970d0, 0.002086d0, 0.01275d0, > -26.79d0, 0.06328d0, 0.03622d0, 0.08345d0, > 39.72d0, -0.06009d0, -0.07825d0, -0.9698d0, > 0.000178d0, -0.000573d0, -0.9328d0, -872.5d0, > -5851.d0, -39.68d0, 4.902d0, -3.848d0, > 2.790d0, 20.91d0, 6.193d0, 26.81d0/ c data (tga(5,j),j = 1,32)/ > -1.539d0, -14.29d0, 3.479d0, 53.36d0, > -0.004201d0, -0.2043d0, -0.03932d0, 0.6409d0, > -0.1058d0, -3.221d0, -0.00114d0, 0.02166d0, > -30.43d0, 0.04049d0, 0.05464d0, 0.008884d0, > 42.00d0, -0.01035d0, -0.1053d0, -1.630d0, > 0.003802d0, -0.001029d0, 4.204d0, -665.6d0, > -1011.0d0, -43.49d0, 4.514d0, -2.948d0, > 2.99d0, 21.59d0, 6.005d0, 22.00d0/ c data (tga(6,j),j = 1,32)/ > -2.581d0, -7.726d0, 5.045d0, 53.31d0, > 0.02262d0, -0.1972d0, -0.01981d0, 0.4280d0, > -0.1055d0, -5.075d0, 0.002762d0, 0.03277d0, > -27.35d0, 0.04986d0, 0.06119d0, -0.1211d0, > 47.48d0, -0.0502d0, -0.1477d0, 0.838d0, > -0.01008d0, -0.0057d0, 9.231d0, -674.3d0, > -900.0d0, -74.43d0, 4.658d0, -3.245d0, > 3.39d0, 21.80d0, 5.620d0, 25.17d0/ c data ip /100/ c data fc /0.3183099031d0/ c data rt /30.d0/ c data x1 /4.d0/ data x2 /5.d0/ c data psi /10.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c if (indval .eq. ip) go to 2 c ip = indval c do 10 j = 1, 32 pa(j) = tga(ip,j) 10 continue c c1 = pa(29)**2 rrc2 = pa(27)**2 dstr = pa(26) / rrc2 * 4.d0 xn = pa(28) rh = pa(31) dy = pa(30) delx = pa(32) b0 = pa(23) b1 = pa(24) b2 = pa(25) xn21 = (xn-x1)**2 xn2 = xn - x2 xnr = 1.d0 / xn2 xn22 = xn2**2 adln = log(xn22 / xn21) c 2 continue c if (abs(ps - psi) .lt. 1.d-6) go to 3 c psi = ps sps = sin(ps) cps = cos(ps) rps = rh * sps c c ---------------------------------------------------------- c*FON Le traitement commence ici si les parametress indval et ps c*FON restent inchanges apres l'appel precedant c ---------------------------------------------------------- c 3 continue c zs = z - rps zp = z - rt zm = z + rt fy = fc / (1.d0 + (y / dy)**2) xnx = xn - x xnx2 = xnx**2 xc1 = x - x1 xc2 = x - x2 xc22 = xc2**2 xr2 = xc2 * xnr xc12 = xc1**2 b20 = zs**2 + c1 b2p = zp**2 + c1 b2m = zm**2 + c1 b = sqrt(b20) bp = sqrt(b2p) bm = sqrt(b2m) xa1 = xc12 + b20 xap1 = xc12 + b2p xam1 = xc12 + b2m xa2 = 1.d0 / (xc22 + b20) xap2 = 1.d0 / (xc22 + b2p) xam2 = 1.d0 / (xc22 + b2m) xna = xnx2 + b20 xnap = xnx2 + b2p xnam = xnx2 + b2m f = b20 - xc22 fp = b2p - xc22 fm = b2m - xc22 xln1 = log(xn21 / xna) xlnp1 = log(xn21 / xnap) xlnm1 = log(xn21 / xnam) xln2 = xln1 + adln xlnp2 = xlnp1 + adln xlnm2 = xlnm1 + adln aln = 0.25d0 * (xlnp1 + xlnm1 -2.d0 * xln1) s0 = (atan(xnx / b) + pid)/ b s0p = (atan(xnx / bp) + pid) / bp s0m = (atan(xnx / bm) + pid) / bm s1 = (xln1 * .5d0 + xc1 * s0) / xa1 s1p = (xlnp1 * .5d0 + xc1 * s0p) / xap1 s1m = (xlnm1 * .5d0 + xc1 * s0m) / xam1 s2 = (xc2 * xa2 * xln2 - xnr - f * xa2 * s0) * xa2 s2p = (xc2 * xap2 * xlnp2 - xnr - fp * xap2 * s0p) * xap2 s2m = (xc2 * xam2 * xlnm2 - xnr - fm * xam2 * s0m) * xam2 g1 = (b20 * s0 - 0.5d0 * xc1 * xln1) / xa1 g1p = (b2p * s0p - 0.5d0 * xc1 * xlnp1) / xap1 g1m = (b2m * s0m - 0.5d0 * xc1 * xlnm1) / xam1 g2 = ((0.5d0 * f * xln2+ 2.d0 * s0 * b20 * xc2) > * xa2 + xr2) * xa2 g2p = ((0.5d0 * fp * xlnp2 + 2.d0 * s0p * b2p * xc2) > * xap2 + xr2) * xap2 g2m = ((0.5d0 * fm * xlnm2 + 2.d0 * s0m *b2m * xc2) > * xam2 + xr2) * xam2 bx = fy * (b0 * (zs * s0 - 0.5d0 * (zp * s0p + zm * s0m)) > + b1 * (zs * s1 - 0.5d0 * (zp * s1p + zm * s1m)) > + b2 * (zs * s2 - 0.5d0 * (zp * s2p + zm * s2m))) by = 0.d0 bz = fy * (b0 * aln + b1 *(g1 - 0.5d0 * (g1p + g1m)) > + b2 * (g2 - 0.5d0 * (g2p + g2m))) c c --------------------------------------------------------- c*FON La contribution courante au calcul de la queue magnetique c*FON est terminee c --------------------------------------------------------- c ex1 = exp(x / delx) ex2 = ex1**2 y2 = y**2 z2 = z**2 yz = y * z c bx = bx + (ex1 * pa(1) + ex2 * pa(3)) * z * cps > + (ex1 * pa(2) + ex2 * (pa(4) + pa(5) * y2 > + pa(6) * z2)) * sps by = (ex1 * pa(7) + ex2 * pa(9)) * yz * cps > + (ex1 * pa(8) + ex2 * (pa(10) + pa(11) * y2 > + pa(12) * z2)) * y * sps bz = bz+(ex1 *(pa(13) + pa(14) * y2 + pa(15) * z2) > + ex2 * (pa(17) + pa(18) * y2 + pa(19) * z2)) > * cps + (ex1 * pa(16) + ex2 * (pa(20) + pa(21) > * y2 + pa(22) * z2)) * z * sps c c ------------------------------------------------------------ c*FON Les contributions dcf et fac ont ete ajoutees a bx, by et bz c ------------------------------------------------------------ c xsm = x * cps - z * sps zsm = x * sps + z * cps z2 = zsm**2 rr2 = xsm**2 + y**2 zn = sqrt((rr2 + z2) / rrc2 + 4.d0)**5 brsm = dstr * 3.d0*zsm/zn bzsm = dstr * (2.d0 * z2-rr2+8.d0*rrc2)/zn by1 = brsm * y bxsm = brsm * xsm c c -------------------------------------------------- c*FON Reconstitution des composantes du champ magnetique c -------------------------------------------------- c bx = bx + bxsm * cps + bzsm * sps bz = bz - bxsm * sps + bzsm * cps by = by + by1 c c **************** c Fin de programme c **************** c return end