diff --git a/iteration_functions.h b/iteration_functions.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6ea2c8aedf50e1cd8dfcc607756699e057bfcec
--- /dev/null
+++ b/iteration_functions.h
@@ -0,0 +1,242 @@
+#include "local_header.h"
+#include "globals.h"
+#include "common_header.h"
+#include "gsl/gsl_odeiv2.h"
+#include "define_values.h"
+
+
+#ifndef ITERATION_FUNCIONS_H_
+#define ITERATION_FUNCIONS_H_
+
+#define SEED_CENTER 1
+#define SEED_UNIFORM 2
+#define SEED_BASE 3
+#define SEED_DELTA 4
+#define SEED_EXP 5
+
+#define RTE 1
+#define MC_SLAB 2
+
+#define U_POSITIVE 1
+#define U_NEGATIVE -1
+
+extern FILE* dat_spectrum;
+extern FILE* dat_diffusion;
+extern FILE* dat_integral_polar;
+extern FILE* dat_energy_polar;
+
+extern int seed_location;
+extern int Nstep_mu;
+extern int Nstep_tau;
+
+extern double tau0;
+extern double tau_c;
+
+extern gsl_spline2d *Spline_I2d_l;
+extern gsl_spline2d *Spline_I2d_r;
+
+
+extern gsl_interp_accel* xacc_2d;
+extern gsl_interp_accel* yacc_2d;
+
+extern gsl_spline* Spline_sample_Iu;
+extern gsl_spline* Spline_eval_limbdark;
+extern gsl_spline* Spline_eval_pdeg;
+extern gsl_spline *Spline_tau_Ak;
+extern gsl_spline *Spline_tau_Bk;
+extern gsl_spline *Spline_tau_Ck;
+
+extern gsl_interp_accel* xacc_1d;
+extern gsl_interp_accel* yacc_1d;
+
+extern char* polarization_file;
+extern char* diffusion_file;
+extern char* integral;
+extern char* polarinput;
+
+typedef struct {
+
+ double* array_mu;
+ double* array_tau;
+ double** matrix;
+
+
+
+} Ikl_intensity;
+
+typedef struct {
+
+ double* array_mu;
+ double* array_tau;
+ double** matrix;
+
+} Ikr_intensity;
+
+
+typedef struct {
+
+ int idx_tau;
+ double * weights_u;
+ double u;
+ double* array_u;
+ double* array_tau;
+ double** Il_matrix_upstream;
+ double** Ir_matrix_upstream;
+
+ double** Il_matrix_downstream;
+ double** Ir_matrix_downstream;
+
+ gsl_spline2d *Spline_I2d_l_upstream;
+ gsl_spline2d *Spline_I2d_l_downstream;
+ gsl_spline2d *Spline_I2d_r_upstream;
+ gsl_spline2d *Spline_I2d_r_downstream;
+
+
+} Iklr_intensity;
+
+
+typedef struct {
+
+ int k;
+ int seed_location;
+ double tau_s;
+ double tau0;
+ double tau;
+ double* array_mu;
+ double* array_tau;
+ double** matrix;
+
+} AngularFunctionParams;
+
+
+typedef struct {
+
+ int nlines;
+ double u_min;
+ double u_max;
+ double* u;
+ double* pdeg;
+ double* Q;
+ double* U;
+ double* limbdark;
+ gsl_spline* Spline;
+ gsl_interp_accel* acc;
+
+
+} radiation_field_t;
+
+typedef struct {
+ double *array_rand;
+ double *array_u;
+ double *array_integ;
+ double total_integral;
+ int nlines;
+ gsl_spline *spline;
+ gsl_interp_accel *acc1;
+ } cumulative_angdist_t;
+
+
+
+/*========================================================================*/
+/*Functions prototypes*/
+/*========================================================================*/
+
+
+double f_interp(double u, double tau);
+double I0_center(double tau0, double tau, double u);
+double I0_delta(double tau, double tau_s, double mu);
+double I0_uniform(double tau0, double tau, double mu);
+double I0_base(double tau0, double tau, double mu, int sign_mu);
+double I0_exp(double tau, double tau0, double mu);
+
+double theta_funct(double tau, double tau_c);
+void make_linarray(double a, double b, int nstep, double* my_array);
+double** dmatrix(long nrl, long nrh, long ncl, long nch);
+
+double Ak_integral (double x, void * params);
+double Bk_integral (double x, void * params);
+double Ck_integral (double x, void * params);
+
+double Il_integral (double x, void * params);
+double Ir_integral (double x, void * params);
+
+void Compute_Ak_array(AngularFunctionParams* ptr_data,
+ double** Il_funct, double* za, double*array_tau, double*array_mu, int k,
+ gsl_integration_workspace * w, double* Ak_array);
+
+void Compute_Bk_array(AngularFunctionParams* ptr_data,
+ double** Il_funct, double* za, double*array_tau, double*array_mu, int k,
+ gsl_integration_workspace * w, double* Bk_array);
+
+void Compute_Ck_array(AngularFunctionParams* ptr_data,
+ double** Ir_funct, double* za, double*array_tau, double*array_mu, int k,
+ gsl_integration_workspace * w, double* Ck_array);
+
+double round_to_digits(double value, int digits);
+
+void slab_polarization(int kiter, int seed_location);
+
+void read_resultsfile(char* filename, radiation_field_t* radiation_field, int Flag);
+void compute_cumulative_angdist(radiation_field_t* radiation_field, cumulative_angdist_t* cumulative_angdist);
+
+void check_results(gsl_spline* Spline_sample_Iu, gsl_spline* Spline_eval_pdeg, gsl_spline* Spline_eval_limbdark);
+
+/* ========================================================================= */
+/* Functions used to integrate the system of two ODE */
+/* ========================================================================= */
+
+
+extern gsl_spline* Spline_I2d_l_upstream;
+extern gsl_spline* Spline_I2d_l_downstream;
+extern gsl_spline* Spline_I2d_r_upstream;
+extern gsl_spline* Spline_I2d_l_downstream;
+
+double legendre_integration_A(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
+double legendre_integration_B(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
+double legendre_integration_C(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
+double legendre_integration_D(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
+
+
+double scattering_matrix (int row, int col, double u, double u_prime);
+
+int solve_rte(int k_iter, int seed_distribution, int Nstep_tau);
+int RTE_Equations(double s, const double y[], double f[], void* params);
+
+
+void compute_results(int method,
+ int k_iter,
+ int Nstep_tau,
+ int Nstep_mu,
+ double* array_mu,
+ double* weights_u,
+ Iklr_intensity* ptr_Iklr,
+ Ikl_intensity* ptr_Ikl,
+ Ikr_intensity* ptr_Ikr);
+
+
+
+
+//==============================================================
+//Routines related to local MC
+//==============================================================
+
+int slab_mc(int nphot, int seed);
+void compute_stokes(double ene, double* k_lab, double* polvect_cart,
+ stokes_parameters* ptr_stokes, double* array_ubounds, int flag);
+
+
+void collect_photons(double ene,
+ double* k,
+ double* array_ene,
+ double** array_counts,
+ int flag);
+
+void setup_gsl_objects(radiation_field_t* radiation_field, int Flag);
+void electricfield_fromfile (double* k_phot, double* elect_field_lab, double* mag_field_lab, double *Pdeg);
+
+double subrel_maxwellian (double kte);
+
+void optimize_taustep(int k_iter, int seed_distribution);
+
+
+#endif