sphere()

void sphere(double s_gp[SDIV+1], 
	    double log_e_tab[201], 
	    double log_p_tab[201], 
	    double log_h_tab[201],
	    double log_n0_tab[201], 
	    int n_tab,                 
	    char eos_type[],
	    double Gamma_P, 
	    double e_center,
	    double p_center, 
	    double h_center,
	    double p_surface,
	    double e_surface,
	    double **rho,
	    double **gama,
	    double **alpha,
	    double **omega,
	    double *r_e)

{
 int s,
     m,
     n_nearest;

 double r_is_s,
        r_is_final,
        r_final, 
        m_final,
        lambda_s,
        nu_s,
        r_is_gp[RDIV+1],
        lambda_gp[RDIV+1],
        nu_gp[RDIV+1],
        gama_mu_0[SDIV+1],
        rho_mu_0[SDIV+1],
        gama_eq,
        rho_eq,
        s_e=0.5;

 /* The function TOV integrates the TOV equations. The function
	can be found in the file equil.c */

 TOV(1, eos_type, e_center, p_center, p_surface, e_surface, Gamma_P,
              log_e_tab, log_p_tab, log_h_tab, n_tab, r_is_gp, lambda_gp, 
              nu_gp, &r_is_final, &r_final, &m_final);

 TOV(2, eos_type, e_center, p_center, p_surface, e_surface, Gamma_P,
              log_e_tab, log_p_tab, log_h_tab, n_tab, r_is_gp, lambda_gp, 
              nu_gp, &r_is_final, &r_final, &m_final);

 TOV(3, eos_type, e_center, p_center, p_surface, e_surface, Gamma_P,
              log_e_tab, log_p_tab, log_h_tab, n_tab, r_is_gp, lambda_gp, 
              nu_gp, &r_is_final, &r_final, &m_final);



 n_nearest=RDIV/2;
 for(s=1;s<=SDIV;s++) {
    r_is_s=r_is_final*(s_gp[s]/(1.0-s_gp[s]));

    if(r_is_s<r_is_final) {
      lambda_s=interp(r_is_gp,lambda_gp,RDIV,r_is_s,&n_nearest);
      nu_s=interp(r_is_gp,nu_gp,RDIV,r_is_s,&n_nearest);
    }
    else {
      lambda_s=2.0*log(1.0+m_final/(2.0*r_is_s));
      nu_s=log((1.0-m_final/(2.0*r_is_s))/(1.0+m_final/(2*r_is_s)));
    }

    gama[s][1]=nu_s+lambda_s;
    rho[s][1]=nu_s-lambda_s;

    for(m=1;m<=MDIV;m++) {
        gama[s][m]=gama[s][1];        
        rho[s][m]=rho[s][1];
        alpha[s][m]=(gama[s][1]-rho[s][1])/2.0;
        omega[s][m]=0.0; 
    }
 
    gama_mu_0[s]=gama[s][1];                   /* gama at \mu=0 */
    rho_mu_0[s]=rho[s][1];                     /* rho at \mu=0 */

 }

   n_nearest=SDIV/2;
   gama_eq = interp(s_gp,gama_mu_0,SDIV,s_e,&n_nearest); /* gama at equator */
   rho_eq = interp(s_gp,rho_mu_0,SDIV,s_e,&n_nearest);   /* rho at equator */
 
   (*r_e)= r_final*exp(0.5*(rho_eq-gama_eq)); 

}