diff --git a/legendre_integration.c b/legendre_integration.c
index 0bd565659d3def7a2e70c1c787dff007bdd8695d..3183e5f5dff72f3eaac7bb6ebcbcf354df772dbe 100644
--- a/legendre_integration.c
+++ b/legendre_integration.c
@@ -33,6 +33,9 @@ double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, do
status = gsl_spline2d_eval_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)
{
printf("Warning: interpolated Ik_l < 0 for tau=%5.4f u_prime=%5.4f\n", tau, u_prime[ii]);
@@ -43,7 +46,7 @@ double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, do
{
printf("Error at line %d\n", __LINE__);
printf("File %s\n", __FILE__);
- printf("tau=%5.4f u=%5.4f\n", tau, u_prime[ii]);
+ printf("tau=%5.4f u=%5.4f I^{k}(l)=%5.4e\n", tau, u_prime[ii], Ik_l);
exit(1);
}
diff --git a/main_program.c b/main_program.c
index 75196361c2c95231d177292579182d31a43974b5..aa9bf093b71ca598ba357527bdb37ce4747bb9cb 100644
--- a/main_program.c
+++ b/main_program.c
@@ -213,7 +213,7 @@ int main(int argc, char* argv[])
exit(1);
}
- solve_rte(kiter, seed);
+ optimize_taustep(kiter, seed);
}
else if (strcmp(method, "st85") == 0)
{
diff --git a/solve_rte.c b/solve_rte.c
index e57a9fbd7a8b8334eac17b09620f342a3dcc8a31..a39d4e8f7cf79596702347f4303487a482d7ee96 100644
--- a/solve_rte.c
+++ b/solve_rte.c
@@ -1,4 +1,5 @@
#include "iteration_functions.h"
+#include <gsl/gsl_interp2d.h>
double I0_center_updown(double tau, double u);
double I0_upward(double tau, double u);
@@ -11,12 +12,28 @@ double tau0;
/*================================================================================*/
-void solve_rte(int k_iter, int seed_distribution)
+void optimize_taustep(int k_iter, int seed_distribution)
{
+ Nstep_tau = 100;
+
+ tau0 = disktau / 2.;
+
+ while (solve_rte(k_iter, seed_distribution, Nstep_tau))
+ {
+ Nstep_tau = (int)Nstep_tau * 1.2;
+ }
+}
+
+/*================================================================================*/
+
+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_tau = 100;
+ /* Nstep_tau = 100;*/
tau0 = disktau / 2.;
@@ -47,7 +64,6 @@ void solve_rte(int k_iter, int seed_distribution)
{
array_mu[ii - 1] = roots[ii];
array_weights[ii - 1] = weights[ii];
-
// printf("%lf 1 \n", array_mu[ii - 1]);
}
@@ -109,6 +125,9 @@ void solve_rte(int k_iter, int seed_distribution)
}
}
}
+
+
+
/*=================================================================*/
/*Setup for 2D interpolation*/
/*=================================================================*/
@@ -129,10 +148,8 @@ void solve_rte(int k_iter, int seed_distribution)
xacc_2d = gsl_interp_accel_alloc();
yacc_2d = gsl_interp_accel_alloc();
- step_tau = (2 * tau0) / (Nstep_tau - 1);
+ step_tau = (2 * tau0) / Nstep_tau;
- // printf("Step tau %5.4f\n", step_tau);
- /*double check_val=0;*/
Iklr_intensity* ptr_data = malloc(sizeof(Iklr_intensity));
Iklr_intensity* params;
@@ -173,13 +190,23 @@ void solve_rte(int k_iter, int seed_distribution)
ptr_data->Spline_I2d_r_upstream = Spline_I2d_r_upstream;
ptr_data->Spline_I2d_r_downstream = Spline_I2d_r_downstream;
+
+
+
+ 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*/
/*=======================================================================*/
status = GSL_SUCCESS + 1;
- gsl_set_error_handler_off();
+ //gsl_set_error_handler_off();
for (jj = 0; jj < Nstep_mu; jj++)
{
@@ -204,53 +231,26 @@ void solve_rte(int k_iter, int seed_distribution)
y[2] = 0;
y[3] = 0;
- /* h_step=1e-7;
- gsl_odeiv_evolve * e_rte = gsl_odeiv_evolve_alloc (4);
- gsl_odeiv_control * c_rte = gsl_odeiv_control_y_new (1e-8, 1e-8);
- gsl_odeiv_step * s_rte = gsl_odeiv_step_alloc (T, 4);
- */
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-4, 1e-4);
+ gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, step_tau, 1e-7, 1e-7);
while (tau < 2 * tau0)
{
- status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
-
- while (status != GSL_SUCCESS)
- {
- tau = 0;
- step_tau = step_tau * 0.75;
- printf("step_tau now is %lf\n", step_tau);
-
- // gsl_odeiv2_driver_free(d);
+ // printf("tau prima %lf\n", tau);
- gsl_odeiv2_driver* d =
- gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, step_tau, 1e-3, 1e-3);
-
- status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
- }
-
- idx = set_array_idx(tau_prev, tau, array_tau, Nstep_tau);
+ status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
if (status)
{
printf("Error during integration: %s\n", gsl_strerror(status));
- printf("Status %d\n", status);
- exit(1);
+ // printf("Status %d\n", status);
+ return (1);
}
- // printf("tau_prev=%lf tau=%lf idx %d\n", tau_prev, tau, idx);
-
- // status = gsl_odeiv_evolve_apply(e_rte, c_rte, s_rte, &sys, &tau, 2*tau0, &h_step, y);
-
- // printf("tau=%lf y[0]=%4.3e y[1]=%4.3e y[2]=%4.3e y[3]=%4.3e idx %d\n", tau, y[0], y[1],
- // y[2], y[3], idx);
+ idx = set_array_idx(tau_prev, tau, array_tau, Nstep_tau);
if (idx > 0)
{
- // printf("tau=%lf y[0]=%4.3e y[1]=%4.3e y[2]=%4.3e y[3]=%4.3e idx %d\n", tau, y[0], y[1],
- // y[2], y[3], idx);
-
ptr_Iklr[kk].Il_matrix_upstream[idx][jj] = y[0];
ptr_Iklr[kk].Ir_matrix_upstream[idx][jj] = y[1];
@@ -263,10 +263,6 @@ void solve_rte(int k_iter, int seed_distribution)
gsl_odeiv2_driver_free(d);
- // gsl_odeiv_evolve_free (e_rte);
- // gsl_odeiv_control_free (c_rte);
- // gsl_odeiv_step_free (s_rte);
-
} // End of loop over array_mu
} // End of loop over k_iter
@@ -276,7 +272,20 @@ void solve_rte(int k_iter, int seed_distribution)
compute_results(1, k_iter, Nstep_tau, Nstep_mu, array_mu, array_weights, ptr_Iklr, ptr_a, ptr_b);
- exit(1);
+ free(Il_upstream_za);
+ free(Il_downstream_za);
+ free(Ir_upstream_za);
+ free(Ir_downstream_za);
+ free(array_tau);
+ free(ptr_Iklr);
+
+ gsl_spline2d_free(Spline_I2d_l_upstream);
+ gsl_spline2d_free(Spline_I2d_l_downstream);
+ gsl_spline2d_free(Spline_I2d_r_upstream);
+ gsl_spline2d_free(Spline_I2d_r_downstream);
+
+
+ return 0;
}
/*================================================================================*/
@@ -295,15 +304,14 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
gsl_spline2d* Spline_I2d_r_upstream = ((Iklr_intensity*)params)->Spline_I2d_r_upstream;
gsl_spline2d* Spline_I2d_r_downstream = ((Iklr_intensity*)params)->Spline_I2d_r_downstream;
- if (s > 2 * tau0 - 1e-5)
- {
- s = 2 * tau0;
- }
/*==============================================*/
/*Equation for I_l upstream*/
/*==============================================*/
+ //printf("UNO s=%10.7f\n", 2*tau0);
+
+
f[0] = -1 / u * y[0] +
3 / (8 * u) *
(legendre_integration_A(Spline_I2d_l_upstream, s, u, array_u, weights_u) +
@@ -312,10 +320,16 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
legendre_integration_B(Spline_I2d_r_upstream, s, u, array_u, weights_u) +
legendre_integration_B(Spline_I2d_r_downstream, 2 * tau0 - s, u, array_u, weights_u));
+
+
+
+
/*==============================================*/
/*Equation for I_r upstream*/
/*==============================================*/
+
+
f[1] = -1 / u * y[1] +
3 / (8 * u) *
(legendre_integration_C(Spline_I2d_l_upstream, s, u, array_u, weights_u) +
@@ -327,6 +341,8 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
/*Equation for I_l downstream*/
/*==============================================*/
+
+
f[2] = -1 / u * y[2] +
3 / (8 * u) *
(legendre_integration_A(Spline_I2d_l_downstream, s, u, array_u, weights_u) +
@@ -339,6 +355,8 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
/*Equation for I_r downstream*/
/*==============================================*/
+
+
f[3] = -1 / u * y[3] +
3 / (8 * u) *
(legendre_integration_C(Spline_I2d_l_downstream, s, u, array_u, weights_u) +