subroutine bdbspst2k (tilt,rb,gg,rr,thetr,phir,bmxe,bmye,bmze,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 00.03.29 - V 1.1 c*VER 01.06.01 - V 2.0 c*VER 03.01.06 - V 2.1 c* c*AUT spec. CNES - JC KOSIK - janvier 1998 c*AUT port. CISI c* c*ROL Theme : Modeles de champs magnetiques c*ROL Calcul du champ lointain. c* c*PAR tilt (I) : angle de tilt (radians) c* c*PAR rb (I) : distance subsolaire normalisee de la magnetopause c*PAR : (rayons terrestres) c* c*PAR gg (I) : coefficients de l'expansion spherique harmonique (6,6) c* c*PAR rr (I) : distance radiale geocentrique (rayons terrestres) c*PAR thetr (I) : colatitude geocentrique (radians) c*PAR phir (I) : longitude geocentrique (radians) c* c*PAR bmxe (O) : composante xgsm du champ lointain c*PAR bmye (O) : composante ygsm du champ lointain c*PAR bmze (O) : composante zgsm du champ lointain c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c* c*NOT common : util c* c*INF utilise : vspvcar c* c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP c*HST version 1.1 - 00.03.29 - amelioration c*HST version 2.0 - 01.06.01 - 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 tilt,rb double precision gg(6,6) double precision rr,thetr,phir double precision bmxe,bmye,bmze 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 kmax c*LOC kmax : nombre de pas de calcul (= 6) c integer mtot c*LOC mtot : nombre de pas calcules c integer k,l,jj,m c*LOC k,l,jj,m : indices de boucles c double precision g(6,6) c double precision br,bt,bp c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ c double precision dbr(9),dbt(9),dbp(9) c*LOC dbr,dbt,dbp : composantes radiales, tangentielles et azimuthales du champ c double precision thetlim c*LOC thetlim : limite max de la colatitude c double precision sh(6,6),fn(6),fm(6) double precision cp(6),sp(6),p(6,6),dp(6,6),const(6,6) double precision cph,sph,ct,st double precision ar,par double precision temp 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 kmax /6/ c data thetlim /1.d-15/ 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 sh(1,1) /0.d0/ c c ****************** c Debut de programme 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 10 k = 1, kmax fn(k) = k - 1 do 20 l = 1, k fm(l) = l - 1 const(k,l) = dble((k - 2)**2 - (l - 1)**2) > / dble((2 * k - 3) * (2 * k - 5)) 20 continue 10 continue c do 30 k = 2, kmax sh(k,1) = dble(2 * k - 3) * sh(k-1,1) / dble(k - 1) jj = 2 do 40 l = 2, k sh(k,l) = sh(k,l-1) * dsqrt(dble((k - l + 1) * jj) > / (k + l - 2)) jj = 1 40 continue 30 continue c do 50 k = 2, kmax g(k,1) = gg(k,1) * sh(k,1) do 60 l = 2, k g(k,l) = gg(k,l) * sh(k,l) 60 continue 50 continue endif c cph = dcos(phir) sph = dsin(phir) c ct = cos(thetr) st = sin(thetr) c if (thetr .le. thetlim) then c ct = 1.0d0 st = 1.d-15 cph = 1.d0 sph = 0.d0 c else if (thetr .ge. pi) then c ct = -1.d0 st = 1.d-15 cph = 1.d0 sph = 0.d0 c endif c sp(2) = sph cp(2) = cph c do 70 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 70 continue c p(2,1) = ct dp(2,1) = -st p(2,2) = st dp(2,2) = ct c ar = 1.d0 c br = ar * g(2,1) * p(2,1) * cos(tilt) bt = ar * g(2,1) * dp(2,1) * cos(tilt) bp = 0.d0 c m = 1 c dbr(m) = ar * p(2,1) * cos(tilt) dbt(m) = ar * dp(2,1) * cos(tilt) dbp(m) = 0.d0 c do 90 k = 3, kmax c ar = ar * (rr / rb) c do 80 l = 1, k c mtot = k + l 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 (mod(mtot,2) .eq. 1) then c m = m + 1 dbp(m) = 0.d0 c if (l .eq. 1) then temp = 1.d0 else temp = cp(l) dbp(m) = -sp(l) * fm(l) * par * cos(tilt) / st bp = bp + g(k,l) * dbp(m) endif dbr(m) = temp * fn(k) * par * cos(tilt) br = br + g(k,l) * dbr(m) dbt(m) = temp * dp(k,l) * ar * cos(tilt) bt = bt + g(k,l) * dbt(m) c endif c 80 continue c 90 continue c call vspvcar(thetr,phir,br,bt,bp,bmxe,bmye,bmze,ier) c c **************** c Fin de programme c **************** c return end