GCC Code Coverage Report

Directory: ./
File: phys/iniorbit.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 31 31 100.0%
Branches: 4 6 66.7%
Line Branch Exec Source
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
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 DO iter = 1, 100
88 3 zdx = -(zx0-r_ecc*sin(zx0)-zxref)/(1.-r_ecc*cos(zx0))
89
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 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/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 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
104