soltervo.f
6.16 KB
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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
subroutine soltervo (iyear,imonth,iday,ihour,imin,isec,year,
> alfag,alfas,deltas,obliq,tetdip,phidip,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 : Astronomie et calculs d'orbite
c*ROL Calcul de l'ascension droite de Greenwich.
c*ROL Calcul de l'ascension droite et de la declinaison
c*ROL du soleil.
c*ROL Calcul de la position geographique du point ou le dipole
c*ROL coupe l'hemisphere nord du 1er janvier 45 au 31 decembre 69.
c*
c*PAR iyear (I) : annee de la date (annee >= 1945 et < 1970)
c*PAR imonth (I) : mois de la date
c*PAR iday (I) : jour de la date
c*PAR ihour (I) : heures de la date
c*PAR imin (I) : minutes de la date
c*PAR isec (I) : secondes de la date
c*
c*PAR year (O) : annee fractionnaire (>= 1945. et < 1970.)
c*
c*PAR alfag (O) : ascension droite de Greenwich (radians)
c*
c*PAR alfas (O) : ascension droite du soleil (radians)
c*PAR deltas (O) : declinaison du soleil (radians)
c*
c*PAR obliq (O) : obliquite de l'ecliptique (radians)
c*
c*PAR tetdip (O) : colatitude geocentrique du dipole (radians)
c*PAR phidip (O) : longitude geocentrique du dipole (radians)
c*
c*PAR ier (O) : code de retour
c*
c*NOT ier : = 1 si iyear < 1945 ou iyear >= 1970
c*
c*NOT L'annee 1990 represente l'origine des temps pour les
c*NOT calculs du champ magnetique terrestre, il faut donc
c*NOT que le parametre iyear soit superieur ou egal a 1990.
c*
c*INF utilise : julg, sun
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
integer iyear, imonth, iday, ihour, imin, isec
double precision year
double precision alfag, alfas, deltas
double precision obliq, tetdip, phidip
integer ier
c
c ----------------------------------
c*FON Declaration des fonctions externes
c ----------------------------------
c
external julg
integer julg
c
c ---------------------------------
c*FON Declaration des variables locales
c ---------------------------------
c
integer lg(3,6)
integer julc, jul1, nujour, imin1, imax
integer i,k
c*LOC i,k : indices de boucles et de tableaux
c
double precision gg(3,6)
double precision delt
c*LOC delt : intervalle entre l'annee courante et l'annee de reference
c
double precision tzero(6)
c*LOC tzero : tableau des annees de reference de 1945 a 1970
c
double precision fjour,g10,g11,c11,h11,gdt
double precision tm, tmold
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 (lg(i,1),i=1,3) /-30594, -2285, 5810/
data (lg(i,2),i=1,3) /-30554, -2250, 5815/
data (lg(i,3),i=1,3) /-30500, -2215, 5820/
data (lg(i,4),i=1,3) /-30421, -2169, 5791/
data (lg(i,5),i=1,3) /-30334, -2119, 5776/
data (lg(i,6),i=1,3) /-30220, -2068, 5737/
c
data tzero /1945.d0, 1950.d0, 1955.d0, 1960.d0,
> 1965.d0, 1970.d0/
data tmold /0.d0/
c
c ******************
c Debut de programme
c ******************
c
ier = 0
c
if (iyear .lt. 1970 .and. iyear .ge. 1945) then
c
c -----------------------------
c*FON Calcul de la date en secondes
c -----------------------------
c
fjour = dble ((ihour * 60 + imin) * 60 + isec) / 86400.d0
c
c -------------------------------
c*FON Calcul du jour julien CNES 1950
c -------------------------------
c
julc = julg (iday,imonth,iyear,ier)
c
jul1 = julg(1,1,iyear,ier)
c
c -----------------------------------------------
c*FON Calcul du numero du jour compte a partir du 1er
c*FON janvier de l'annee courante
c -----------------------------------------------
c
nujour = julc - jul1 + 1
c ----------------------------------------------
c*FON Calcul de la position du soleil dans le repere
c*FON equatorial a gamma 50
c ----------------------------------------------
c
call sun (iyear,nujour,fjour,alfag,alfas,deltas,
+ obliq,ier)
year = dble(iyear) + dble(nujour-1)/365.d0
c
c ---------------------------------------------------
c*FON Calculs intermediaires pour obtenir les coordonnees
c*FON du dipole
c ---------------------------------------------------
c
tm = dble(iyear) + dble(imonth) / 12.d0
if (tm .ne. tmold) then
c
tmold = tm
imin1 = idint((tm - 1945.d0) / 5.d0) + 1
imax = imin1 + 1
delt = (tm - tzero(imin1)) / 5.d0
c
do 30 k = 1, 3
gg(k,imin1) = dble(lg(k,imin1))
> + dble(lg(k,imax)
> - lg(k,imin1)) * delt
30 continue
c
g10 = - gg(1,imin1)
g11 = - gg(2,imin1)
h11 = - gg(3,imin1)
c
c11 = sqrt(g11**2 + h11**2)
gdt = sqrt(g11**2 + h11**2 + g10**2)
c
c ------------------------------------------
c*FON Calcul de la longitude et de la colatitude
c*FON du dipole en radians
c ------------------------------------------
c
phidip = atan(h11 / g11)
tetdip = asin(c11 / gdt)
c
endif
else
ier = 1
endif
c
c ****************
c Fin de programme
c ****************
c
return
end