diff --git a/compute_results.c b/compute_results.c
index 2c899a541acfdd82941a6ae64e748b98c990b0a5..bee1f7739ffd05ea2a97600ccdc49c3e6535d8a9 100644
--- a/compute_results.c
+++ b/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],
diff --git a/legendre_integration.c b/legendre_integration.c
index 3183e5f5dff72f3eaac7bb6ebcbcf354df772dbe..81a76191f4576a4df0b0dc7fcf99f7a6c00d5f39 100644
--- a/legendre_integration.c
+++ b/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,10 +20,12 @@ 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++)
{
angular_term = scattering_matrix(1, 1, u, u_prime[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)
{
diff --git a/main_program.c b/main_program.c
index aa9bf093b71ca598ba357527bdb37ce4747bb9cb..f585dc189d584ccc4fea6c6e3ea16f33a8c9a188 100644
--- a/main_program.c
+++ b/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);
}
diff --git a/mc_slab.c b/mc_slab.c
index 70204eab9083e2edc2b2abb5c8e11e8812ac8327..a127c0c21d804c33c316702fbe60a2513c779e83 100644
--- a/mc_slab.c
+++ b/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;
diff --git a/solve_rte.c b/solve_rte.c
index a39d4e8f7cf79596702347f4303487a482d7ee96..44139b8009e452f0c6bec189f3194dda16f05f92 100644
--- a/solve_rte.c
+++ b/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*/