%****************************************************************************
% 	calcul et trace les rayon en provenancedu GPS immobile avec un multitrajet
%****************************************************************************

clear all;

NN =input (' entrer nombre de couches : ');
layer =input (' numero de la couche ou vous avez realiser un multitrajet: ');
Nv =input (' nouvelle valeur de la couche : ');


[n,r,nombre_cercle] = f_refraction2(50);
n(NN-layer+1)=Nv;

figure(1)
plot(r,n(1:nombre_cercle),'+-');

% rayon Terre
ae=6378.137*10^3;

% rayon orbite LEO
r_leo=ae+800*10^3;

% rayon orbite GPS
r_gps=ae+20000*10^3;

% echantillonnage angle_rad des cercles
NN=500;
nn=1:NN;
delta_teta_cercle = (2*pi)/NN;
teta_cecle_rad = nn * delta_teta_cercle;

% cercle de rayon r(1),r(2)...
cx=ones(nombre_cercle,NN);
cy=ones(nombre_cercle,NN);

for u=1:nombre_cercle,
	cx(u,:)=r(u)*cos(teta_cecle_rad);
	cy(u,:)=r(u)*sin(teta_cecle_rad);
end

LEO_x=ones(NN);
LEO_y=ones(NN);

LEO_x=r_leo*cos(teta_cecle_rad);
LEO_y=r_leo*sin(teta_cecle_rad);


% position du LEO sur son orbite
flag=1;

xl=ae+630*10^3
yl=sqrt(r_leo^2-xl^2)

%raytracing using bouger's rule

[phi_L_ini] = g_position(xl,yl,ae,r_gps)

figure(2)

for k=1:nombre_cercle,
  plot(cx(k,:),cy(k,:),'r:');
  hold on;
end


plot(LEO_x,LEO_y,'g--');
hold on;

flag_touch = 0;
k=0;

Xv=ones(500,1000);
Yv=ones(500,1000);


while (flag_touch == 0)
   phi_L_inibis=phi_L_ini+k*0.0000010;
   [flag,flag_touch,X,Y,L,delta_t,rt,delta_phase] = g_raytracing(phi_L_inibis,r_leo,r_gps,xl,yl,r,n);
   k=k+1;
end

Xv(1,1:L)=X;
Yv(1,1:L)=Y;

DT(1)=delta_phase;
T(1)=delta_t;
Rt(1)=rt;

plot(X,Y,'m-');
hold on;

k=1;
j=2;
phi_L = phi_L_inibis;
i=nombre_cercle;

while (flag == 1)
   
   while ( (r(i)<=rt) & (rt<r(i-1)) )
      
       phi_L_all=phi_L+k*0.0001;
       [flag,flag_touch,X,Y,L,delta_t,rt,delta_phase] = g_raytracing(phi_L_all,r_leo,r_gps,xl,yl,r,n);
       k=k+1
       
   end
    
   Xv(j,1:L)=X;
   Yv(j,1:L)=Y;
   DT(j)=delta_phase;
   T(j)=delta_t;
   Rt(j)=rt;
   plot(X,Y,'-');
   hold on;
   j=j+1;
   i=i-1;
   if (i < 2)
      break;
   end
   

end

G=length(DT);

figure (3)
plot(DT(1:G-1),'+-');

figure (4)
plot(Rt(1:G-1)-ae,'+-');

figure (5)

for k=1:nombre_cercle,
  plot(DT(1:G-1),(r(k)-ae)*ones(1,G-1),'r');
  hold on;
end

plot(DT(1:G-1),Rt(1:G-1)-ae,'+-');