#include "prototype.h"

/****** Function g_raytracing */

void g_raytracing(int *flag, int *flag_touch, double *X, double *Y, double *delta_t, double *rt, 
             double *delta_phase, double phi_L_rad, double r_leo, double r_gps, 
             double xl, double yl, double r, double n,int nombre_cercle)
 {
   
  
   int N,i, t, j, flag_R, cpt;
   double phi_prim, angle_rad;
   double *phi, *tps


   phi=(double *) malloc (MAX * sizeof(double));
   tps=(double *) malloc (MAX * sizeof(double));


     N=nombre_cercle+1;
    

     g_ini(flag, &(X[0]), &(Y[0]), &angle_rad, &(phi[0]), phi_L_rad, r_leo, r[0], xl, yl); 
     f_delay(&(tps[0]), xl, yl, x[0], y[0], n[0]);

     i=1;
     cpt=1;

     flag_R=1;
     *flag_touch=1;]
   
     *delta_t=0;

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

	     i=i+1;
	     cpt=cpt+1;

	     if( (flag_R==1) && (i==nombre_cercle) )
	       {
		 *flag_touch=0;
		 *rt=0;
		 break;
	       }
	   }

	 if(*flag_touch == 1)
	   {
	     g_snell(&flag_R, x[i-1], y[i-1], &angle_rad, &(phi[i-1]), &phi_prim,
	     angle_rad, phi[i-2], n[i-2], n[i-1], r[i-2], r[i-2]);
	     f_delay(&(tps[i-1]), 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 > 2)
	       {
		 g_snell2(&flag_R, x[j], y[j], &angle_rad, &(phi[j]),
		      angle_rad, phi[j-1], n[t-1], n[t-2], r[t-2], r[t-3]);
		 f_delay(&(tps[j]), x[j], y[j], x[j-i], y[j-1], n[t-2]);
		
		 j=j+1;
		 cpt=cpt+1;
		 t=t-1;
	       }
 
	     g_snell2(&flag_R, x[j], y[j], &angle_rad, &(phi[j]),
		      angle_rad, phi[j-1], n[1], n[0], r[1], r_gps);
	     f_delay(&(tps[j]), x[j], y[j], x[j-i], y[j-1], n[0]);
	   }

	 else
	   {
	     X[0]=xl;
	     Y[0]=yl;
	     *rt=0;
	   }

	 for(i=0; i<cpt ; i++)

	   *delta_t=*delta_t+tps[i];

	 RLG=sqrt((X[cpt-1]-xl)*(X[cpt-1]-xl)+(Y[cpt-1]-yl)*(Y[cpt-1]-yl));
	 *delta_phase=*delta_t-(RLG/C);

       }