subroutine lines (r1,r2,r3,b,arc,err,j,vp,vn,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.06.05 - V 2.0 c*VER 03.01.06 - V 2.1 c* c*AUT spec. Mac Ilwain c*AUT adap. CNES - JC KOSIK - juin 1995 c*AUT port. CISI c* c*ROL Theme : Calculs de geophysique c*ROL Calcul de la ligne de champ passant par les points. c*ROL Logiciel modifie pour tenir compte du champ total. c* c*PAR r1 (I/O) : valeur du champ magnetique B pour le 1 er point c*PAR r2 (I/O) : valeur du champ magnetique B pour le 2 eme point c*PAR r3 (I/O) : coordonnees du 3eme point c* c*PAR b (O) : champ magnetique c* c*PAR arc (I/O) : tableau des arcs c* c*PAR err (I) : parametre d'erreurs c* c*PAR j (O) : nombre de points de la ligne de force c* c*PAR vp (O) : coordonnees du 2 eme point c*PAR vn (O) : coordonnees du 3 eme point c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c* c*NOT common : util c* c*INF utilise : gsfc65 c* c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP c*HST version 2.0 - 01.06.05 - 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 double precision year double precision r1(3), r2(3), r3(3) double precision b(200) double precision arc(200) double precision err integer j double precision vn(3), vp(3) 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 ilp,is,i c*LOC ilp,is,i : indices de boucles c integer ier1 c*LOC ier1 : code retour des modules appeles c double precision br,bt,bp c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ c double precision ra(3) double precision a1,a2,a3,aab,aa,ad,am,ao6,arcj,asum double precision bb,bd,cc,cd,cre,dd,dn double precision pre1,pre2,pre3,qrt,r,rbar,rt,sit,sna,x 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 year /1965.d0/ c c ****************** c Debut de programme c ****************** c ier = 0 c cre = 0.25d0 c if (err - 0.15625d0) 10,20,20 c 10 continue c cre = (err**0.333333333d0) c 20 continue c a3 = arc(3) aab = abs (a3) sna = a3 / aab a1 = arc(1) a2 = arc(2) ao6 = a3 * a3 / 6.0d0 j = 3 ilp = 1 is = 1 c go to 60 c 30 continue c is = 1 j = j + 1 ao6 = a3 * a3 / 6.0d0 arcj = a1 + a2 + a3 ad = (asum + a1) / aa bd = asum / bb cd = a1 / cc c 40 continue c do 50 i = 1 ,3 c dd = r1(i) / aa - r2(i) / bb + r3(i) / cc if (is .ne. 2) then c rt = r1(i) - (ad * r1(i) - bd * r2(i) + cd * r3(i) - > dd * arcj) * arcj ra(i) = r1(i) r1(i) = r2(i) r2(i) = r3(i) r3(i) = rt vp(i) = vn(i) c endif c rbar = (r2(i) + r3(i)) / 2.d0 - dd * ao6 vn(i) = vp(i) + a3 * rbar c 50 continue c 60 continue c if (vn(2)) 70,80,80 c 70 continue c vn(2) = -vn(2) c 80 continue c if (vn(2) - pi) 100,100,90 c 90 continue c vn(2) = dpi - vn(2) c go to 80 c 100 continue c if (vn(3)) 110,120,120 c 110 continue c vn(3) = vn(3) + dpi c go to 100 c 120 continue c if (vn(3) - dpi) 140,140,130 c 130 continue c vn(3) = vn(3) - dpi c go to 120 c 140 continue c if (is .eq. 2) go to 150 c sit = abs(sin(vn(2))) pre1 = vn(1) pre2 = pre1 * vn(2) pre3 = pre1 * sit * vn(3) r = vn(1) c c -------------------------------------------- c*FON Calcul du champ magnetique d'origine interne c -------------------------------------------- c call gsfc65(year,r,vn(2),vn(3),br,bt,bp,b(j),ier1) c r3(1) = br / b(j) dn = b(j) * vn(1) r3(2) = bt / dn r3(3) = bp / (dn * sit) asum = a3 + a2 aa = asum * a2 bb = a3 * a2 cc = asum * a3 is = 2 c go to 40 c 150 continue c sit = abs(sin(vn(2))) b(j) = b(j) * ((pre1 / vn(1))**3) c qrt = .5d0 * abs(r3(1)) / (.1d0 + abs(r3(2) * vn(1))) x = (abs(vn(1) - pre1) + qrt * abs(vn(1) * vn(2) - pre2) + > abs(vn(1) * sit * vn(3) - pre3)) / > (aab * err * sqrt(1.d0 + qrt * qrt)) c if (ilp .ne. 2) go to 180 c if (x - 3.3d0) 180,160,160 c 160 continue c a3 = a3 * 0.2d0 * (8.0d0 + x) / (0.8d0 + x) j = j - 1 ilp = 3 asum = a2 + a1 aa = asum * a1 bb = a2 * a1 cc = asum * a2 c do 170 i = 1, 3 c vn(i) = vp(i) r3(i) = r2(i) r2(i) = r1(i) r1(i) = ra(i) c 170 continue c go to 250 c 180 continue c if (j - 200) 190,260,260 c 190 continue c a1 = a2 if (b(j) - b(2)) 200,200,260 c 200 continue c ilp = 2 a2 = a3 a3 = a3 * .2d0 * (8.d0 + x) / (.8d0 + x) am = (2.d0 - r3(2) * vn(1)) * vn(1) * cre if (abs(a3) - am) 220,220,210 c 210 continue c a3 = sna * am c 220 continue c if (sna * r3(1) + .5d0) 230,230,250 c 230 continue c am = -.5d0 * sna * vn(1) / r3(1) if (abs(a3) - am) 250,250,240 c 240 continue c a3 = sna * am c 250 continue c arc(j+1) = a3 aab = abs(a3) c go to 30 c 260 continue c c **************** c Fin de programme c **************** c return end