function [flag,flag_touch,X,Y,L,delta_t,rt,delta_phase] = g_raytracing(phi_L_rad,r_leo,r_gps,xl,yl,r,n)

N=length(n);
c=2.998*10^8;
%c=3*10^8;

[flag,x(1),y(1),angle_rad,phi(1)] = g_ini(phi_L_rad,r_leo,r(1),xl,yl);
[tps(1)] = f_delay(xl,yl,x(1),y(1),n(1));

i=2;
flag_R = 1;
flag_touch=1;

if (flag == 1)
   
 while ( (flag_R == 1) & (i < N) )
    
   
    [flag_R,x(i),y(i),angle_rad,phi(i),phi_prim] = g_snell(angle_rad,phi(i-1),n(i-1),n(i),r(i-1),r(i));
    [tps(i)] = f_delay(x(i),y(i),x(i-1),y(i-1),n(i));
 
    i=i+1;

   
   if ( (flag_R == 1) & (i == N) )
      flag_touch=0;
      rt=0;
   break;
   end

end

if ( flag_touch == 1)
[flag_R,x(i-1),y(i-1),angle_rad,phi(i-1),phi_prim] = g_snell(angle_rad,phi(i-2),n(i-2),n(i-1),r(i-2),r(i-2));
[tps(i-1)] = f_delay(x(i-1),y(i-1),x(i-2),y(i-2),n(i-1));
rt=r(i-2)*sin(phi_prim);

t=i;
j=i;

while ( t > 3 )
[flag_R,x(j),y(j),angle_rad,phi(j)] = g_snell2(angle_rad,phi(j-1),n(t-1),n(t-2),r(t-2),r(t-3));
[tps(j)] = f_delay(x(j),y(j),x(j-1),y(j-1),n(t-2));

j=j+1;
t=t-1;
end

[flag_R,x(j),y(j),angle_rad,phi(j)] = g_snell2(angle_rad,phi(j-1),n(2),n(1),r(1),r_gps);
[tps(j)] = f_delay(x(j),y(j),x(j-1),y(j-1),n(1));

end

else
   
   x=xl;
   y=yl;
   rt=0;
end

X=[xl x];
Y=[yl y];

L=length(X);
delta_t=sum(tps);

RLG=sqrt((X(L)-xl)^2+(Y(L)-yl)^2);
delta_phase=delta_t-(RLG/c);