diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000000000000000000000000000000000000..27f1f01c527bf0f0cc11901c95f941d42aef86e8
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,50 @@
+Language: Cpp
+AccessModifierOffset: -4
+AlignEscapedNewlinesLeft: true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: true
+AllowShortFunctionsOnASingleLine: None
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AlwaysBreakBeforeMultilineStrings: true
+AlwaysBreakTemplateDeclarations: true
+BinPackParameters: false
+BreakBeforeBinaryOperators: false
+BreakBeforeBraces: Allman
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+ColumnLimit: 100
+CommentPragmas: ''
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+ConstructorInitializerIndentWidth: 4
+ContinuationIndentWidth: 4
+DerivePointerAlignment: false
+DisableFormat: false
+ExperimentalAutoDetectBinPacking: false
+IndentCaseLabels: false
+IndentWidth: 2
+IndentWrappedFunctionNames: false
+IndentFunctionDeclarationAfterType: false
+MaxEmptyLinesToKeep: 1
+KeepEmptyLinesAtTheStartOfBlocks: false
+NamespaceIndentation: None
+PenaltyBreakBeforeFirstCallParameter: 1
+PenaltyBreakComment: 300
+PenaltyBreakString: 1000
+PenaltyBreakFirstLessLess: 120
+PenaltyExcessCharacter: 1
+PenaltyReturnTypeOnItsOwnLine: 1000
+PointerAlignment: Left
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeParens: ControlStatements
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 1
+SpacesInAngles: false
+SpacesInCStyleCastParentheses: false
+SpacesInContainerLiterals: true
+SpacesInParentheses: false
+Cpp11BracedListStyle: true
+Standard: Cpp11
+TabWidth: 4
+UseTab: Never
diff --git a/compute_stokes.c b/compute_stokes.c
index 532df030815e655e5029b82b01ea54a682efb88d..4c468619581594a56e5b6e674fee071a9c1f8964 100644
--- a/compute_stokes.c
+++ b/compute_stokes.c
@@ -91,6 +91,12 @@ void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_para
double* Pdeg = (double*)malloc(nstepangles * sizeof(double));
dat_integral_polar = fopen(integral, "w");
+ fprintf(dat_integral_polar, "!Col1: cos(theta)\n");
+ fprintf(dat_integral_polar, "!Col2: Pdeg\n");
+ fprintf(dat_integral_polar, "!Col3: Q\n");
+ fprintf(dat_integral_polar, "!Col4: U\n");
+ fprintf(dat_integral_polar, "!Col5: Itot/I(u=1)*1/u\n\n");
+
for (jj = 0; jj < nstepangles; jj++)
{
u_center = (array_ubounds[jj] + array_ubounds[jj + 1]) / 2;
diff --git a/main_program.c b/main_program.c
index 375a1ffd69c42eba8db23b66403afa22ebe61f78..d29571be27cdad0f3d0949632d32bf43b01820ef 100644
--- a/main_program.c
+++ b/main_program.c
@@ -2,24 +2,29 @@
#include <mpi.h>
#include <stdio.h>
+char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method);
+
int seed;
double kte;
double ktbb;
double tau_c;
double disktau;
+double albedobase;
+double gmin;
+double gmax;
+double p1;
+
char* polarization_file;
char* diffusion_file;
char* polardisk;
-double gmin;
-double gmax;
-double p1;
int FlagEleDist;
int nph;
int rank;
-//int nproc;
+
int kiter;
int DiskPolarizationFlag;
+int Albedoflag;
char* outspec;
char* polarfile;
@@ -27,8 +32,6 @@ char* integral;
char* diffusion;
char* polarinput;
-char* set_filename(char* root, char* tau_char, char* seed_char);
-
FILE* dat_spectrum;
FILE* dat_diffusion;
FILE* dat_integral_polar;
@@ -57,7 +60,7 @@ int main(int argc, char* argv[])
printf("Shared parameters\n\n");
printf("disktau=vaue: slab optical depth\n\n");
- printf("\nSeed photons distributions: three options are possible\n");
+ printf("\nSeed photons distributions: three options are possible with parameter seed=N\n");
printf("1: photons at the center of the slab\n");
printf("2: photons uniform distribution\n");
printf("3: photons at the base (valid only for ode and mc methods)\n\n");
@@ -85,6 +88,7 @@ int main(int argc, char* argv[])
char* tau_char;
char* seed_char;
+ char* albedo_char;
for (t = 1; t < argc; t++)
{
@@ -155,6 +159,13 @@ int main(int argc, char* argv[])
{
ktbb = numpar(argv[t]);
}
+
+ if (strstr(argv[t], "albedobase=") != NULL)
+ {
+ albedobase = numpar(argv[t]);
+ albedo_char = stringpar(argv[t]);
+ Albedoflag = 1;
+ }
}
/*=========================================================*/
@@ -167,26 +178,6 @@ int main(int argc, char* argv[])
return 1;
}
- if (diffusion == 0x0)
- {
- diffusion = "mc_diffusion.qdp";
- }
-
- if (integral == 0x0)
- {
- integral = "integral_polar_mc.qdp";
- }
-
- if (polarfile == 0x0)
- {
- polarfile = "energy_polar_mc.qdp";
- }
-
- if (outspec == 0x0)
- {
- outspec = "espectrum.qdp";
- }
-
if (!kte)
{
kte = 1;
@@ -213,8 +204,8 @@ int main(int argc, char* argv[])
if (strcmp(method, "ode") == 0)
{
- diffusion_file = set_filename("diffusion_ode_tau", tau_char, seed_char);
- polarization_file = set_filename("intensity_ode_tau", tau_char, seed_char);
+ diffusion_file = set_filename("diffusion_ode_tau", tau_char, seed_char, albedo_char, "ode");
+ polarization_file = set_filename("intensity_ode_tau", tau_char, seed_char, albedo_char, "ode");
if (kiter == 0)
{
@@ -226,8 +217,8 @@ int main(int argc, char* argv[])
}
else if (strcmp(method, "st85") == 0)
{
- polarization_file = set_filename("intensity_st85_tau", tau_char, seed_char);
- diffusion_file = set_filename("diffusion_st85_tau", tau_char, seed_char);
+ polarization_file = set_filename("intensity_st85_tau", tau_char, seed_char, albedo_char, "st85");
+ diffusion_file = set_filename("diffusion_st85_tau", tau_char, seed_char, albedo_char, "st85");
if (kiter == 0)
{
@@ -239,8 +230,8 @@ int main(int argc, char* argv[])
}
else if (strcmp(method, "mc") == 0)
{
- MPI_Init(&argc, &argv);
- init_mpi();
+ MPI_Init(&argc, &argv);
+ init_mpi();
if (polardisk == 0x0 && seed == 3 && rank == 0)
{
@@ -253,6 +244,18 @@ int main(int argc, char* argv[])
polardisk = "n";
}
+ if (Albedoflag != 1)
+ {
+ printf("Please select the albedo at the base of the slab with albedobase=value\n");
+ exit(1);
+ }
+
+ if (albedobase < 0 || albedobase > 1)
+ {
+ printf("Error: albedo must be in the range 0-1\n");
+ exit(1);
+ }
+
if (nph == 0)
{
printf("Please select the number of photons for MC simulation with parametere nph=nph\n");
@@ -279,7 +282,13 @@ int main(int argc, char* argv[])
exit(1);
}
+ diffusion = set_filename("mc_diffusion_tau", tau_char, seed_char, albedo_char, "mc");
+ outspec = set_filename("mc_spectrum_tau", tau_char, seed_char, albedo_char, "mc");
+ integral = set_filename("mc_integralpolar_tau", tau_char, seed_char, albedo_char, "mc");
+ polarfile = set_filename("mc_energypolar_tau", tau_char, seed_char, albedo_char, "mc");
+
status = slab_mc(nph, seed);
+
// check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);
}
@@ -295,7 +304,7 @@ int main(int argc, char* argv[])
return (0);
}
-char* set_filename(char* root, char* tau_char, char* seed_char)
+char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method)
{
// root="intensity_ode";
@@ -303,6 +312,7 @@ char* set_filename(char* root, char* tau_char, char* seed_char)
char* tmp_a;
char* tmp_b;
char* tmp_c;
+ char* tmp_d;
tmp_a = stringconcat(root, tau_char);
@@ -316,9 +326,20 @@ char* set_filename(char* root, char* tau_char, char* seed_char)
// printf("tmp_c = %s\n", tmp_c);
- filename = stringconcat(tmp_c, ".dat");
+ if (strcmp(method, "mc") == 0)
+ {
+ tmp_d = stringconcat(tmp_c, stringconcat("_A", albedo));
+
+ filename = stringconcat(tmp_d, ".qdp");
+ }
+ else
+ {
+ filename = stringconcat(tmp_c, ".qdp");
+ }
// printf("Final name %s\n\n", filename);
return filename;
}
+
+/*==============================================================================*/
diff --git a/mc_slab.c b/mc_slab.c
index 5279bcb0e3c806bea2c9f6bb83f54ca0a716b65b..5cd29c754c438623984565a16d3dc47b6a0e7915 100644
--- a/mc_slab.c
+++ b/mc_slab.c
@@ -2,6 +2,8 @@
#define NSC_MAX 1000
+extern double albedobase;
+
/*==================================================================*/
/*Function prototypes*/
/*==================================================================*/
@@ -40,17 +42,14 @@ int nstepangles;
char* matrix_gb;
-
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_maxw;
-gsl_root_fsolver *s_maxw;
-
+gsl_root_fsolver* s_maxw;
/*==================================================================*/
@@ -92,13 +91,12 @@ int slab_mc(int nphot, int seed)
T_phi = gsl_root_fsolver_bisection;
s_phi = gsl_root_fsolver_alloc(T_phi);
-
+
/*This is for sampling beta from subrelativistics maxwellian*/
T_maxw = gsl_root_fsolver_bisection;
s_maxw = gsl_root_fsolver_alloc(T_maxw);
-
nstepangles = 40;
obsmindeg = 0.;
obsmaxdeg = 90.;
@@ -303,6 +301,8 @@ int slab_mc(int nphot, int seed)
Magfield_lab[2] = Magfield_lab_new[2];
}
else
+
+ /*Photon potentially escaping from the slab */
{
if (k_lab[2] > 0)
{
@@ -313,18 +313,46 @@ int slab_mc(int nphot, int seed)
collect_photons(E_new_lab, k_lab, array_ene, array_counts, 0);
nph_esc++;
+ FlagLocation = 0;
}
else
{
- // printf("Photon %d escaping\n", ii);
- }
+ csi = clock_random();
- FlagLocation = 0;
- }
+ /*Photon absorbed and newly re-emitted with probability A=albedobase*/
+ if (csi <= albedobase)
+ {
+ // printf("Photon %d reflected\n", ii);
+
+ cos_theta = clock_random();
+ sin_theta = sqrt(1 - cos_theta * cos_theta);
+ phi = 2 * PI * clock_random();
+
+ k_lab[0] = sin_theta * cos(phi);
+ k_lab[1] = sin_theta * sin(phi);
+ k_lab[2] = cos_theta;
+
+ pos[0] = 0;
+ pos[1] = 0;
+ pos[2] = 0;
+
+ unpolarized_seedphot(k_lab, Efield_lab, Magfield_lab);
+ PolarFlag = 0;
+ E_lab = bb_spec(ktbb);
+ }
+
+ /*Photon absorbed*/
+ else
+ {
+ FlagLocation = 0;
+ }
+ } // Close the else with condition k_lab[2] < 0
+
+ } // End of confition if (tau > eps)... else...
} // end of while (FlagLocation)
- // newphot:;
+
}
/*==================================================================*/
@@ -426,9 +454,9 @@ void comptonization(double E_lab,
versor_boost(k_lab, el_vel_vect, k_erf, 1);
- double beta_new=subrel_maxwellian (kte);
+ double beta_new = subrel_maxwellian(kte);
- //printf("pluto %lf %lf\n", beta_el, beta_new);
+ // printf("pluto %lf %lf\n", beta_el, beta_new);
elmag_lorentztransform(Efield_lab, Magfield_lab, el_vel_vect, ElectField_old_erf, MagField_Erf, 1);
diff --git a/run_mc.sh b/run_mc.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b55613d212a3d0679bc1abed45466675ecd048ac
--- /dev/null
+++ b/run_mc.sh
@@ -0,0 +1,43 @@
+array_tau=(0.2 0.5 1 2 5 8)
+array_seed=(2 3)
+array_albedo=(0 0.5 1)
+
+ii="0"
+jj="0"
+kk="0"
+
+nph=10000
+
+#Loop over tau
+while [ $ii -lt ${#array_tau[@]} ]
+do
+
+ #Loop over seed distribution
+ jj="0"
+ while [ $jj -lt ${#array_seed[@]} ]
+ do
+
+ #Loop over slabrmax
+ kk="0"
+ while [ $kk -lt ${#array_albedo[@]} ]
+ do
+
+
+ command="start_iteration method=mc disktau=${array_tau[ii]} seed=${array_seed[jj]} polardisk=n nph=$nph albedobase=${array_albedo[kk]}"
+
+ echo $command
+
+ $command
+
+ kk=$[$kk+1]
+ done
+
+ jj=$[$jj+1]
+ done
+
+ ii=$[$ii+1]
+ done
+
+
+
+