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