CDPP / AMDA_COTS

Blame view

maglib/src/k899k.f 3.79 KB
Edit Raw Normal View History
eb2d6c5c   Hacene SI HADJ MOHAND   adding maglib 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
 
      subroutine k899k (tilt,facrc,c20,rre,thet,phi,bre,bte,bpe,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. CNES - JC KOSIK - janvier 1998
c*AUT port. CISI
c*
c*ROL Theme : Modeles de champs magnetiques
c*ROL         Calcul du courant de l'anneau du champ magnetique.
c*
c*PAR tilt  (I) : angle de tilt (radians)
c*
c*PAR facrc (I) : coefficient du courant de l'anneau
c*PAR c20   (I) : coefficient du courant de l'anneau
c*
c*PAR rre   (I) : distance radiale (kilometres ou rayons terrestres) 
c*PAR thet  (I) : colatitude geocentrique GSM (radians)
c*PAR phi   (I) : colatitude geocentrique GSM  (radians)
c*
c*PAR bre   (O) : distante radiale geocentrique (idem input)
c*PAR bte   (O) : colatitude geocentrique (radians)
c*PAR bpe   (O) : colatitude geocentrique (radians)
c*
c*PAR ier   (O) : code de retour
c*
c*NOT ier       : sans objet
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.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 tilt, facrc, c20, rre, thet, phi
      double precision bre, bte, bpe
      integer ier 
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      double precision br,bt,bp
c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ
c
      double precision brtilt,bttilt,bptilt
c*LOC brtilt,bttilt,bptilt : composantes radiale, tangentielle et azimuthale
c
      double precision c10,c21,xk1,xk2,xk2t,r2,r3
      double precision fons1,fons2,fons2t,ct,st,s2t,c2t,cp,sp
      double precision dfons1, dfons2, dfons2t
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     Debut de programme
c     ******************
c
      ier = 0
c
      c10 = -1.5d0
      c21 =  0.11d0
c
      xk1  = 0.04d0
      xk2  = 0.01d0
      xk2t = 0.005d0
c
      r2 = rre * rre
      r3 = r2 * rre
c
      fons1  = r3 * dexp(-xk1 * r2)
      fons2  = r3 * dexp(-xk2 * r2)
      fons2t = r2 * dexp(-xk2t * r2)
c
      ct  = dcos(thet)
      st  = dsin(thet)
      s2t = 2.d0 * st * ct
      c2t = ct * ct - st * st
      cp  = dcos(phi)
      sp  = dsin(phi)
c
      dfons1  = (4.d0 - 2.d0 * xk1  * r2) * fons1  / rre
      dfons2  = (4.d0 - 2.d0 * xk2  * r2) * fons2  / rre
      dfons2t = (3.d0 - 2.d0 * xk2t * r2) * fons2t / rre
c
      br = 2.d0 * c10 * fons1*ct
     >     / rre + 6.d0 * c21 * ct * st * cp *  fons2 / rre
      bt = -c10 * dfons1 * st + 1.732d0 * c21 * c2t * cp * dfons2
      bp = -1.732d0 * c21 * ct * sp * dfons2
c
      br = br * facrc
      bt = bt * facrc
      bp = bp * facrc
c
      brtilt = +c20 * (6.d0 * c2t - 2.d0) * dsin(tilt)
     >         * facrc * fons2t / rre
      bttilt = -c20 * s2t * dsin(tilt) * dfons2t * facrc
      bptilt = 0.d0
c
      bre = br * dcos(tilt) + brtilt
      bte = bt * dcos(tilt) + bttilt
      bpe = bp * dcos(tilt) + bptilt      
c
c     ****************
c     Fin de programme
c     ****************
c
 999    continue
      return
      end