Latest working version of the code

compute_results.c

@@ -94,17 +94,13 @@ void compute_results(int method,
     }
   }
   else
-	  
+  {  
 	  fprintf(dat, "!Pdeg band 1\n");
 	    fprintf(dat, "!Pdeg band 2\n");
 	    fprintf(dat, "!Pdeg band 3\n");
 	    fprintf(dat, "!Pdeg band 4\n");
 	      fprintf(dat, "!Pdeg band 5\n");
   
-	     
-	   
-	  
-  {
     for (jj = 0; jj < Nstep_mu; jj++)
     {
       fprintf(dat, "%lf  %5.4e  %5.4e %5.4e  %5.4e  %5.4e  %5.4e  %5.4e \n", array_mu[jj], Pdeg[jj],

legendre_integration.c

@@ -12,6 +12,7 @@ Gauss-Legendre quadrature
 /* [2*(1-u^2)*(1-u'^2) + u^2 u'^2] I_l(u)  du' */
 /*========================================================*/
 
+
 double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u)
 {
   int status;
@@ -19,9 +20,11 @@ double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, do
   double integ_value;
   double angular_term;
   double Ik_l;
+  double tolerance;
 
   integ_value = 0;
   Ik_l = 0;
+  tolerance=1e-7;
   
   for (ii = 0; ii < Nstep_mu; ii++)
   {
@@ -31,15 +34,15 @@ double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, do
     /*Evaluate by 2D-interpolation the function I^{k-1}_l (tau, u')*/
     /*====================================================================*/
 
-    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l);
+    status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l);
 
     //printf("u %lf tau %lf Ik_l %5.4e\n", u_prime[ii], tau, Ik_l);
 
 
-    if (Ik_l < 0)
+    if (Ik_l < -tolerance)
     {
-      printf("Warning: interpolated Ik_l < 0 for tau=%5.4f u_prime=%5.4f\n", tau, u_prime[ii]);
-      exit(1);
+      //("Warning: interpolated Ik_l < 0 for tau=%5.4f u_prime=%5.4f\n", tau, u_prime[ii]);
+      //exit(1);
     }
 
     if (status)
@@ -92,7 +95,7 @@ double legendre_integration_B(gsl_spline2d* Spline_I2d, double tau, double u, do
     /*Evaluate by 2D-interpolation the function I^{k-1}_r (tau, u')*/
     /*====================================================================*/
 
-    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r);
+    status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r);
 
     if (status)
     {
@@ -139,7 +142,7 @@ double legendre_integration_C(gsl_spline2d* Spline_I2d, double tau, double u, do
     /*Evaluate by 2D-interpolation the function I^{k-1}_l (tau, u')*/
     /*====================================================================*/
 
-    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l);
+    status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l);
 
     if (status)
     {
@@ -185,7 +188,7 @@ double legendre_integration_D(gsl_spline2d* Spline_I2d, double tau, double u, do
     /*Evaluate by 2D-interpolation the function I^{k-1}_r (tau, u')*/
     /*====================================================================*/
 
-    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r);
+    status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r);
 
     if (status)
     {

main_program.c

@@ -299,7 +299,7 @@ int main(int argc, char* argv[])
     return (1);
   }
 
-  MPI_Finalize();
+  if (strcmp(method, "mc") == 0) MPI_Finalize();
 
   return (0);
 }

mc_slab.c

@@ -268,7 +268,7 @@ int slab_mc(int nphot, int seed)
       }
 
       /*====================================================================*/
-      /*Compute the optical depth from top of bottom of the atmosphere*/
+      /*Compute the optical depth from top or base of the atmosphere*/
       /*=====================================================================*/
 
       tau = NeDisk * SIGMA_T * r_units * L_top;

solve_rte.c

@@ -32,7 +32,7 @@ int solve_rte(int k_iter, int seed_distribution, int Nstep_tau)
   printf("\n\nReceveing Nstep_tau %d\n", Nstep_tau);
   int ii, jj, kk, status;
 
-  Nstep_mu = 30;
+  Nstep_mu = 20;
   /* Nstep_tau = 100;*/
 
   tau0 = disktau / 2.;
@@ -193,12 +193,12 @@ int solve_rte(int k_iter, int seed_distribution, int Nstep_tau)
 
 
 
-    printf("CHECK HERE !\n");
+   /* printf("CHECK HERE !\n");
     double zio=0;
         int status = gsl_spline2d_eval_extrap_e(Spline_I2d_l_downstream, 2*tau0, 0.5, xacc_2d, yacc_2d, &zio);
 
        printf("================> status %d val %lf\n", status, zio);
-
+*/
 
     /*=======================================================================*/
     /*Solve the RTE at the quadrature points*/
@@ -232,13 +232,13 @@ int solve_rte(int k_iter, int seed_distribution, int Nstep_tau)
       y[3] = 0;
 
       gsl_odeiv2_system sys = {RTE_Equations, jac, 4, params};
-      gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, step_tau, 1e-7, 1e-7);
+      gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, step_tau, 1e-4, 1e-4);
 
       while (tau < 2 * tau0)
       {
         // printf("tau prima %lf\n", tau);
 
-        status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
+        status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 1, y);
 
         if (status)
         {
@@ -309,7 +309,7 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
   /*Equation for I_l  upstream*/
   /*==============================================*/
 
-  //printf("UNO  s=%10.7f\n", 2*tau0);
+  //printf("UNO  s=%lf\n");
 
 
   f[0] = -1 / u * y[0] +
@@ -322,7 +322,7 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
 
 
 
-
+  //printf("UNO DOPO  s=%10.7f\n");
 
   /*==============================================*/
   /*Equation for I_r upstream*/