1 |
|
5 |
SUBROUTINE iniorbit(paphelie, pperiheli, pyear_day, pperi_day, pobliq) |
2 |
|
|
IMPLICIT NONE |
3 |
|
|
|
4 |
|
|
! ======================================================================= |
5 |
|
|
|
6 |
|
|
! Auteur: |
7 |
|
|
! ------- |
8 |
|
|
! Frederic Hourdin 22 Fevrier 1991 |
9 |
|
|
|
10 |
|
|
! Objet: |
11 |
|
|
! ------ |
12 |
|
|
! Initialisation du sous programme orbite qui calcule |
13 |
|
|
! a une date donnee de l'annee de duree year_day commencant |
14 |
|
|
! a l'equinoxe de printemps et dont le perihelie se situe |
15 |
|
|
! a la date peri_day, la distance au soleil et la declinaison. |
16 |
|
|
|
17 |
|
|
! Interface: |
18 |
|
|
! ---------- |
19 |
|
|
! - Doit etre appele avant d'utiliser orbite. |
20 |
|
|
! - initialise une partie du common planete.h |
21 |
|
|
|
22 |
|
|
! Arguments: |
23 |
|
|
! ---------- |
24 |
|
|
|
25 |
|
|
! Input: |
26 |
|
|
! ------ |
27 |
|
|
! aphelie \ aphelie et perihelie de l'orbite |
28 |
|
|
! periheli / en millions de kilometres. |
29 |
|
|
|
30 |
|
|
! ======================================================================= |
31 |
|
|
|
32 |
|
|
! ----------------------------------------------------------------------- |
33 |
|
|
! Declarations: |
34 |
|
|
! ------------- |
35 |
|
|
|
36 |
|
|
include "planete.h" |
37 |
|
|
include "YOMCST.h" |
38 |
|
|
|
39 |
|
|
! Arguments: |
40 |
|
|
! ---------- |
41 |
|
|
|
42 |
|
|
REAL paphelie, pperiheli, pyear_day, pperi_day, pobliq |
43 |
|
|
|
44 |
|
|
! Local: |
45 |
|
|
! ------ |
46 |
|
|
|
47 |
|
|
REAL zxref, zanom, zz, zx0, zdx, pi |
48 |
|
|
INTEGER iter |
49 |
|
|
|
50 |
|
|
! ----------------------------------------------------------------------- |
51 |
|
|
|
52 |
|
|
pi = 2.*asin(1.) |
53 |
|
|
|
54 |
|
1 |
aphelie = paphelie |
55 |
|
1 |
periheli = pperiheli |
56 |
|
1 |
year_day = pyear_day |
57 |
|
1 |
obliquit = pobliq |
58 |
|
1 |
peri_day = pperi_day |
59 |
|
|
|
60 |
|
1 |
PRINT *, 'Perihelie en Mkm ', periheli |
61 |
|
1 |
PRINT *, 'Aphelie en Mkm ', aphelie |
62 |
|
1 |
PRINT *, 'obliquite en degres :', obliquit |
63 |
|
1 |
PRINT *, 'Jours dans l annee : ', year_day |
64 |
|
1 |
PRINT *, 'Date perihelie : ', peri_day |
65 |
|
1 |
unitastr = 149.597870 |
66 |
|
1 |
e_elips = (aphelie-periheli)/(periheli+aphelie) |
67 |
|
1 |
p_elips = 0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr |
68 |
|
|
|
69 |
|
1 |
PRINT *, 'e_elips', e_elips |
70 |
|
1 |
PRINT *, 'p_elips', p_elips |
71 |
|
|
|
72 |
|
|
! ----------------------------------------------------------------------- |
73 |
|
|
! calcul de l'angle polaire et de la distance au soleil : |
74 |
|
|
! ------------------------------------------------------- |
75 |
|
|
|
76 |
|
|
! calcul de l'zanomalie moyenne |
77 |
|
|
|
78 |
|
1 |
zz = (year_day-pperi_day)/year_day |
79 |
|
1 |
zanom = 2.*pi*(zz-nint(zz)) |
80 |
|
1 |
zxref = abs(zanom) |
81 |
|
1 |
PRINT *, 'zanom ', zanom |
82 |
|
|
|
83 |
|
|
! resolution de l'equation horaire zx0 - e * sin (zx0) = zxref |
84 |
|
|
! methode de Newton |
85 |
|
|
|
86 |
|
1 |
zx0 = zxref + r_ecc*sin(zxref) |
87 |
✓✗ |
3 |
DO iter = 1, 100 |
88 |
|
3 |
zdx = -(zx0-r_ecc*sin(zx0)-zxref)/(1.-r_ecc*cos(zx0)) |
89 |
✓✓ |
3 |
IF (abs(zdx)<=(1.E-12)) GO TO 120 |
90 |
|
3 |
zx0 = zx0 + zdx |
91 |
|
|
END DO |
92 |
|
|
120 CONTINUE |
93 |
|
1 |
zx0 = zx0 + zdx |
94 |
✗✓ |
1 |
IF (zanom<0.) zx0 = -zx0 |
95 |
|
1 |
PRINT *, 'zx0 ', zx0 |
96 |
|
|
|
97 |
|
|
! zteta est la longitude solaire |
98 |
|
|
|
99 |
|
1 |
timeperi = 2.*atan(sqrt((1.+r_ecc)/(1.-r_ecc))*tan(zx0/2.)) |
100 |
|
1 |
PRINT *, 'longitude solaire du perihelie timeperi = ', timeperi |
101 |
|
|
|
102 |
|
1 |
RETURN |
103 |
|
|
END SUBROUTINE iniorbit |