#include "prototype.h"

/****** Function f_refraction2 */

void f_refraction2(double *n, double *r, int nombre_cercle);
 {
   double Rdry, Cdry, dens0, Ha, ae, delta_z, a, b, hc, N;
   double *ns, *nc, *rc, *re, *rv, *nv, *r, *n, *z;
   int i, j, p, pp;

   Rdry=287.16;
   Cdry=77.6*1.e-2;
   dens0=1.25;
   Ha=7.5;
   ae=6378.137*1.e+3;

   ns=(double *) malloc (nombre_cercle * sizeof(double));
   nc=(double *) malloc ((nombre_cercle-1) * sizeof(double));
   z=(double *) malloc (nombre_cercle * sizeof(double));
   rc=(double *) malloc ((nombre_cercle-1) * sizeof(double));
   re=(double *) malloc ((nombre_cercle-1) * sizeof(double));
   rv=(double *) malloc (nombre_cercle * sizeof(double));
   nv=(double *) malloc ((nombre_cercle+1) * sizeof(double));
   r=(double *) malloc (nombre_cercle * sizeof(double));
   n=(double *) malloc ((nombre_cercle+1) * sizeof(double));

   delta_z=100/nombre_cercle;

   for (i=1 ; i <= nombre_cercle ; i++)
     {
       z[i-1]=i*delta_z;
       N=Rdry*Crdy*dens0*exp(-z/Ha);
       ns[i-1]=(N*1.e-6)+1;
     }

   for (j=0 ; j < nombre_cercle-1 ; j++)
     {
       a=(ns[j+1]-ns[j])/(z[j+1]-z[j]);
       b=ns[j]-a*z[j];
       hc=(z[j+1]-z[j])/2;
       nc[j]=a*hc+b;
       rc[j]=z[j+1];
     }

   for (i=0 ; i<nombre_cercle-1 ; i++)
       re[i]=ae+rc[i]*1.e+3;

   rv[0]=ae;
   nv[0]=1;
   nv[nombre_cercle]=1;

   for (j=1 ; j<nombre_cercle ; j++)
     {
       rv[j]=re[j-1]; 
       nv[j]=rc[j-1];
     }

   p=nombre_cercle-1;
   
   for (k=0 ; k<nombre_cercle ; k++)
     {
       r[k]=rv[p];
       p=p-1;
     }

   pp=nombre_cercle;

   for (k=0 ; k<nombre_cercle+1 ; k++)
     {
       n[k]=nv[pp];
       pp=pp-1;
     }

 }