diff --git a/Makefile b/Makefile index 94d23c8fc12cf56a938add645120559e64ca2bd4..9f0802a0417c4e9ae3a5131f431fef6ad177571d 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ LDFLAGS = -L ${LIB_CFITSIO} -L ${LIB_GSL} -L ${LIB_MPI} -lcfitsio -lm -lgsl -lgslcblas -lmpi INCLUDE = -I ${PWD} -I ${HEADERS_CFITSIO} -I ${HEADERS_GSL} -I ${HEADERS_MPI} -I ${MC_INSTDIR}/../../include/ -CFLAGS = -g -Ofast -Wall ${INCLUDE} +CFLAGS = -g -O3 -Wall ${INCLUDE} LOCAL_DIR = ${PWD} @@ -20,7 +20,7 @@ MC_SYNC = ${MC_INSTDIR}/../../gr_code/synchrotron LOCAL_DIR_OBJ= ${LOCAL_DIR}/main_program.o ${LOCAL_DIR}/check_results.o ${LOCAL_DIR}/compute_results.o ${LOCAL_DIR}/function_integrals.o ${LOCAL_DIR}/legendre_integration.o ${LOCAL_DIR}/set_deltatau.o ${LOCAL_DIR}/slab_polarization.o ${LOCAL_DIR}/solve_rte.o ${LOCAL_DIR}/mc_slab.o ${LOCAL_DIR}/compute_stokes.o ${LOCAL_DIR}/collect_photons.o ${LOCAL_DIR}/simul_sync.o -MC_POLARIZATION_OBJ=${MC_POLARIZATION}/ortonormal_axes.o ${MC_POLARIZATION}/unpolarized_seedphot.o ${MC_POLARIZATION}/polarization_probability.o ${MC_POLARIZATION}/depolarization_probability.o ${MC_POLARIZATION}/electfield_as_new.o ${MC_POLARIZATION}/sample_azim_angle.o ${MC_POLARIZATION}/allocate_stokes_pointer.o ${MC_POLARIZATION}/disk_polarizationelectricfield.o ${MC_POLARIZATION}/pdeg_atmosphere.o ${MC_POLARIZATION}/write_polardata.o ${MC_POLARIZATION}/select_quadrant.o $ ${MC_POLARIZATION}/detector_axis.o +MC_POLARIZATION_OBJ=${MC_POLARIZATION}/ortonormal_axes.o ${MC_POLARIZATION}/unpolarized_seedphot.o ${MC_POLARIZATION}/polarization_probability.o ${MC_POLARIZATION}/depolarization_probability.o ${MC_POLARIZATION}/electfield_as_new.o ${MC_POLARIZATION}/sample_azim_angle.o ${MC_POLARIZATION}/allocate_stokes_pointer.o ${MC_POLARIZATION}/disk_polarizationelectricfield.o ${MC_POLARIZATION}/pdeg_atmosphere.o ${MC_POLARIZATION}/write_polardata.o ${MC_POLARIZATION}/select_quadrant.o ${MC_POLARIZATION}/detector_axis.o MC_ELE_DISTRIBUTION_OBJ=${MC_ELE_DISTRIBUTION}/hybrid/eledist_hybrid.o ${MC_ELE_DISTRIBUTION}/hybrid/integrate_maxwgamma.o ${MC_ELE_DISTRIBUTION}/sample_electron.o ${MC_ELE_DISTRIBUTION}/powerlaw/eledist_nth.o ${MC_ELE_DISTRIBUTION}/maxwellian/subrel_maxwellian.o diff --git a/main_program.c b/main_program.c index 00308f304e6e42b20364d4d28549499a47104e26..f10aea9c8d6d4e42a48daf9b18bea4f26802bb7e 100644 --- a/main_program.c +++ b/main_program.c @@ -1,5 +1,5 @@ -#include "iteration_functions.h" #include "globals.h" +#include "iteration_functions.h" #include <mpi.h> #include <stdio.h> @@ -9,22 +9,20 @@ int seed; double tau_c; double albedobase; - char* polarization_file; char* diffusion_file; - -//int nph; -//int rank; +// int nph; +// int rank; int kiter; int DiskPolarizationFlag; int Albedoflag; -//char* outspec; -//char* polarfile; +// char* outspec; +// char* polarfile; char* integral; -//char* diffusion; +// char* diffusion; char* polarinput; FILE* dat_spectrum; @@ -152,6 +150,16 @@ int main(int argc, char* argv[]) ktbb = numpar(argv[t]); } + if (strstr(argv[t], "emin=") != NULL) + { + emin = numpar(argv[t]); + } + + if (strstr(argv[t], "emax=") != NULL) + { + emax = numpar(argv[t]); + } + if (strstr(argv[t], "albedobase=") != NULL) { albedobase = numpar(argv[t]); @@ -180,6 +188,16 @@ int main(int argc, char* argv[]) ktbb = 1; } + if (!emin) + { + emin = 0.1; + } + + if (!emax) + { + emax = 50.; + } + /*=========================================================*/ if (!disktau || disktau == 0) @@ -274,23 +292,37 @@ int main(int argc, char* argv[]) exit(1); } - diffusion = set_filename("diffusion_mc_tau", tau_char, seed_char, albedo_char, "mc"); - outspec = set_filename("spectrum_mc_tau", tau_char, seed_char, albedo_char, "mc"); - integral = set_filename("integralpolar_mc_tau", tau_char, seed_char, albedo_char, "mc"); - polarfile = set_filename("energypolar_mc_tau", tau_char, seed_char, albedo_char, "mc"); + if (diffusion == 0x0) + { + diffusion = set_filename("diffusion_mc_tau", tau_char, seed_char, albedo_char, "mc"); + } + + if (outspec == 0x0) + { + outspec = set_filename("spectrum_mc_tau", tau_char, seed_char, albedo_char, "mc"); + } + + if (integral == 0x0) + { + integral = set_filename("integralpolar_mc_tau", tau_char, seed_char, albedo_char, "mc"); + } + + if (polarfile == 0x0) + { + polarfile = set_filename("energypolar_mc_tau", tau_char, seed_char, albedo_char, "mc"); + } status = slab_mc(nph, seed); - } else { - printf( - "\nPlease select method=ode (solve ODE) or method=mc (Montecarlo)\n"); + printf("\nPlease select method=ode (solve ODE) or method=mc (Montecarlo)\n"); return (1); } - if (strcmp(method, "mc") == 0) MPI_Finalize(); + if (strcmp(method, "mc") == 0) + MPI_Finalize(); return (0); } diff --git a/mc_slab.c b/mc_slab.c index 8cdf0c6c22c518e1e5e64ade296a78dd9baf5cec..f609a6de5e17f0a277e7c0825379952e628c8b6a 100644 --- a/mc_slab.c +++ b/mc_slab.c @@ -17,7 +17,6 @@ void comptonization(double E_lab, double* Efield_lab_new, double* Magfield_lab_new); - double distance_bottom(double z0, double kz); double distance_top(double z0, double kz, double H); @@ -25,22 +24,20 @@ double distance_top(double z0, double kz, double H); /*Datatypes defined in the header globals.h*/ /*==================================================================*/ -//double obsmindeg; -//double obsmaxdeg; - +// double obsmindeg; +// double obsmaxdeg; -//double* array_angles; -//double* array_ubounds; - -//double** array_gb; -//double** ptr_2colfile; +// double* array_angles; +// double* array_ubounds; +// double** array_gb; +// double** ptr_2colfile; const gsl_root_fsolver_type* T_phi; gsl_root_fsolver* s_phi; -//const gsl_root_fsolver_type* T_disk; -//gsl_root_fsolver* s_disk; +// const gsl_root_fsolver_type* T_disk; +// gsl_root_fsolver* s_disk; const gsl_root_fsolver_type* T_maxw; gsl_root_fsolver* s_maxw; @@ -59,7 +56,6 @@ int slab_mc(int nphot, int seed) int photon_diffusion[NSC_MAX]; int status; - double E_lab; double E_new_lab; double k_lab[3]; @@ -77,7 +73,7 @@ int slab_mc(int nphot, int seed) double Hphot; double L_top; double tau; - //double eps; + // double eps; double l_int; double csi; double Pdeg; @@ -105,7 +101,7 @@ int slab_mc(int nphot, int seed) array_gb = read_matrix(matrix_gb); double array_ene[NSTEP_ENE]; - make_array_energy(array_ene, 0.01, 20, NSTEP_ENE); + make_array_energy(array_ene, emin, emax, NSTEP_ENE); double** array_counts; array_counts = dmatrix(0, NSTEP_ENE - 1, 0, nstepangles - 1); @@ -139,8 +135,7 @@ int slab_mc(int nphot, int seed) setup_gsl_objects(radiation_field, MC_SLAB); } - - Pdeg=0; + Pdeg = 0; nph_esc = 0; cos_theta = 0; @@ -160,7 +155,7 @@ int slab_mc(int nphot, int seed) stokes_parameters struct_stokes_average[POLAR_NBOUNDS - 1]; allocate_stokes_pointer(struct_stokes_average); - n_zero_sc=0; + n_zero_sc = 0; /*=============================================================*/ /*Here starts loop over photon */ /*=============================================================*/ @@ -190,7 +185,10 @@ int slab_mc(int nphot, int seed) cos_theta = 2 * clock_random() - 1; } else if (seed == 2) - { + {if (!emin) + { + emin=0.1; + } pos[2] = clock_random() * Hphot; cos_theta = 2 * clock_random() - 1; } @@ -224,7 +222,6 @@ int slab_mc(int nphot, int seed) k_lab[1] = sin_theta * sin(phi); k_lab[2] = cos_theta; - if (DiskPolarizationFlag == FROMFILE) { electricfield_fromfile(k_lab, Efield_lab, Magfield_lab, &Pdeg); @@ -264,23 +261,22 @@ int slab_mc(int nphot, int seed) pos[1] = pos[1] + k_lab[1] * l_int; pos[2] = pos[2] + k_lab[2] * l_int; - //printf("Photon %d before scattering %d Pdeg = %5.4f\n", ii, n_sc, Pdeg); + // printf("Photon %d before scattering %d Pdeg = %5.4f\n", ii, n_sc, Pdeg); - comptonization(E_lab, k_lab, Efield_lab, Magfield_lab, &E_new_lab, k_new_lab, + comptonization(E_lab, k_lab, Efield_lab, Magfield_lab, &E_new_lab, k_new_lab, Efield_lab_new, Magfield_lab_new); - //printf("After scattering Pdeg = %5.4f\n\n", Pdeg); + // printf("After scattering Pdeg = %5.4f\n\n", Pdeg); n_sc++; E_lab = E_new_lab; - for (jj=0; jj<3; jj++) + for (jj = 0; jj < 3; jj++) { - k_lab[jj] = k_new_lab[jj]; - Efield_lab[jj] = Efield_lab_new[jj]; - Magfield_lab[jj] = Magfield_lab_new[jj]; + k_lab[jj] = k_new_lab[jj]; + Efield_lab[jj] = Efield_lab_new[jj]; + Magfield_lab[jj] = Magfield_lab_new[jj]; } - } else @@ -290,10 +286,16 @@ int slab_mc(int nphot, int seed) { // printf("Photon %d escaping after %d scattering\n", ii, n_sc); - if (n_sc==0) n_zero_sc++; + if (n_sc == 0) + n_zero_sc++; compute_stokes(E_lab, k_lab, Efield_lab, struct_stokes, array_ubounds, 0); - photon_diffusion[n_sc]++; + + if (n_sc < NSC_MAX) + { + photon_diffusion[n_sc]++; + } + collect_photons(E_new_lab, k_lab, array_ene, array_counts, 0); nph_esc++; @@ -364,23 +366,24 @@ int slab_mc(int nphot, int seed) // Get the number of processes - MPI_Comm_size(MPI_COMM_WORLD, &nproc); - + MPI_Comm_size(MPI_COMM_WORLD, &nproc); /*==================================================================*/ + // printf("Rank %d si trova qui %d\n", rank, __LINE__); + if (rank == 0) { compute_stokes(1, k_lab, Efield_lab, struct_stokes_average, array_ubounds, 1); printf("Percentage of escaping photons from the top layer %4.3f\n", 1. * nph_esc / nphot); - printf("Percentage of escaping non scattered photons from the top layer %4.3f\n", (double) n_zero_sc/ nphot); + printf("Percentage of escaping non scattered photons from the top layer %4.3f\n", (double)n_zero_sc / nphot); /*============================================================*/ dat_diffusion = fopen(diffusion, "w"); fprintf(dat_diffusion, "!Percentage of escaping photons %4.3f\n", 1. * nph_esc / nphot); - Fval_max=-10; + Fval_max = -10; for (ii = 0; ii < 300; ii++) { @@ -392,28 +395,35 @@ int slab_mc(int nphot, int seed) for (ii = 0; ii < 300; ii++) { - fprintf(dat_diffusion, " %d %3.2f\n", ii, (double)photon_diffusion[ii]/Fval_max); + fprintf(dat_diffusion, " %d %3.2f\n", ii, (double)photon_diffusion[ii] / Fval_max); } fclose(dat_diffusion); - } - - collect_photons(0, k_lab, array_ene, reduce_array_counts, 1); - dat_energy_polar = fopen(polarfile, "w"); + dat_energy_polar = fopen(polarfile, "w"); + if (dat_energy_polar == NULL) + { + fprintf(stderr, "Errore: fopen fallita su %s\n", polarfile); + MPI_Abort(MPI_COMM_WORLD, 1); + } - write_polardata(struct_stokes_average, ene_bound, dat_energy_polar); + write_polardata(struct_stokes_average, ene_bound, dat_energy_polar); + fclose(dat_energy_polar); + + collect_photons(0, k_lab, array_ene, reduce_array_counts, 1); + } - fclose(dat_energy_polar); + - gsl_root_fsolver_free(s_phi); - gsl_root_fsolver_free(s_disk); + if (s_phi != NULL) + gsl_root_fsolver_free(s_phi); + if (s_disk != NULL) + gsl_root_fsolver_free(s_disk); return 0; } // End of function - void comptonization(double E_lab, double* k_lab, double* Efield_lab, @@ -424,7 +434,6 @@ void comptonization(double E_lab, double* Magfield_lab_new) { - int jj; int status; int PolarFlag; @@ -444,9 +453,8 @@ void comptonization(double E_lab, double ElectField_old_erf[3]; double MagField_Erf[3]; - - PolarFlag=0; - Pdeg=0; + PolarFlag = 0; + Pdeg = 0; //================================================================== */ // Find the K' components in the ERF */ @@ -462,7 +470,6 @@ void comptonization(double E_lab, E_erf = energy_lorentz_boost(E_lab, el_vel_vect, k_lab, 1); - csi = clock_random(); costheta_erf = (-1 + pow(-2 + 4 * csi + sqrt(5 - 16 * csi + 16 * pow(csi, 2)), 2 / 3.)) / @@ -476,7 +483,6 @@ void comptonization(double E_lab, E_new_erf = new_energy_erf(E_erf, costheta_erf); - phik = sample_azim_angle(E_erf, E_new_erf, costheta_erf, 1); //================================================================== */ @@ -496,9 +502,8 @@ void comptonization(double E_lab, // and compute the new electric field vector //================================================================== */ - PolarFlag=1; + PolarFlag = 1; status = depolarization_probability(E_erf, E_new_erf, costheta_erf, phik, &Pdeg, &PolarFlag); - if (status == 1) {