%****************************************************************************
% 	calcul et trace les rayons allant du GPS au LEO durant une RO. model 2D
%****************************************************************************

clear all;

%**************************************************************
%	modelisation index de refraction constant ds chaque couche
%	modele standard
%**************************************************************

NN =input (' entrer nombre de couches : ');

[n,r,nombre_cercle] = f_refraction2(NN);

%**************************************************************
%	rayon de la terre en metre [m]
% 	rayon orbite LEO
%**************************************************************

ae=6378.137*10^3;

r_leo=ae+800*10^3;

r_gps=ae+20000*10^3;


%**************************************************************
% position du LEO sur son orbite
% on peut choisir n'importe quelle coor
% on obtiendra une Radio Occultation
%**************************************************************

xl_in=ae+500*10^3;
yl_in=+sqrt(r_leo^2-xl_in^2);

%**************************************************************
% Initiallisation
%**************************************************************

[phi_L_ini] = g_position(xl_in,yl_in,ae,r_gps)

flag_touch = 0;
k=0;
while (flag_touch == 0)
teta_inibis= phi_L_ini +k*0.00010;
[flag,flag_touch,X_ini,Y_ini,L_ini,delta_t_0,rt_0,delta_phase_0] = g_raytracing(teta_inibis,r_leo,r_gps,xl_in,yl_in,r,n)

k=k+1;
end

xg_ini=X_ini(L_ini);
yg_ini=Y_ini(L_ini);


%**************************************************************
%	 Raytracing signal a travers differentes couches
%	choix periode des mesures du recepteur LEO
%**************************************************************

X=ones(200,1000);
P=ones(200,1000);

km=7/4;
tp_leo=km/7; %vitesse leo 7km/s, tp est le tps parcouru par leo pour la valeur de km
% une mesure tte les tp_leo [Hz] 

t(1) = 0;
teta_out(1) = teta_inibis ;
delta_t(1) = delta_t_0


i=2

while (flag ~= 0) %tant qu'on ne sort pas de l'atmosphere

[flag,t(i),d(i),xg(i),yg(i),teta_out(i),Xa,Ya,delta_t_a(i),L(i),Pa] = 
f_detection2(teta_out(i-1),delta_t(i-1),xl_in,yl_in,i,xg_ini,yg_ini,tp_leo,r_gps,n,r);
X(i,1:L(i))=Xa;
Y(i,1:L(i))=Ya;

end