Skip to content

Commit

Permalink
Update gammapkt.cc
Browse files Browse the repository at this point in the history
  • Loading branch information
@lukeshingles
lukeshingles committed Oct 26, 2024
1 parent 9aecd73 commit faafca5
Showing 1 changed file with 27 additions and 27 deletions.
54 changes: 27 additions & 27 deletions gammapkt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,8 +57,8 @@ std::vector<NucGammaLine> allnuc_gamma_line_list;
void read_gamma_spectrum(const int nucindex, const char filename[50])
// reads in gamma_spectra and returns the average energy in gamma rays per nuclear decay
{
printout("reading gamma spectrum for Z=%d A=%d from %s...", decay::get_nuc_z(nucindex), decay::get_nuc_a(nucindex),
filename);
printoutf("reading gamma spectrum for Z=%d A=%d from %s...", decay::get_nuc_z(nucindex), decay::get_nuc_a(nucindex),
filename);

FILE *filein = fopen_required(filename, "r");
int nlines = 0;
Expand All @@ -79,11 +79,11 @@ void read_gamma_spectrum(const int nucindex, const char filename[50])

decay::set_nucdecayenergygamma(nucindex, E_gamma_avg);

printout("nlines %d avg_en_gamma %g MeV\n", nlines, E_gamma_avg / MEV);
printoutf("nlines %d avg_en_gamma %g MeV\n", nlines, E_gamma_avg / MEV);
}

void set_trivial_gamma_spectrum(const int nucindex) {
// printout("Setting trivial gamma spectrum for z %d a %d engamma %g\n", z, a, decay::nucdecayenergygamma(z, a));
// printoutf("Setting trivial gamma spectrum for z %d a %d engamma %g\n", z, a, decay::nucdecayenergygamma(z, a));
const int nlines = 1;
gamma_spectra[nucindex].resize(nlines, {});
gamma_spectra[nucindex][0].energy = decay::nucdecayenergygamma(nucindex);
Expand All @@ -93,13 +93,13 @@ void set_trivial_gamma_spectrum(const int nucindex) {
void read_decaydata() {
// migrate from old filename
if (!std::filesystem::exists("ni56_lines.txt") && std::filesystem::exists("ni_lines.txt")) {
printout("Moving ni_lines.txt to ni56_lines.txt\n");
printoutf("Moving ni_lines.txt to ni56_lines.txt\n");
std::rename("ni_lines.txt", "ni56_lines.txt");
}

// migrate from old filename
if (!std::filesystem::exists("co56_lines.txt") && std::filesystem::exists("co_lines.txt")) {
printout("Moving co_lines.txt to co56_lines.txt\n");
printoutf("Moving co_lines.txt to co56_lines.txt\n");
std::rename("co_lines.txt", "co56_lines.txt");
}

Expand Down Expand Up @@ -129,7 +129,7 @@ void read_decaydata() {
} else if (std::ifstream(filename2)) {
read_gamma_spectrum(nucindex, filename2);
} else if (decay::nucdecayenergygamma(nucindex) > 0.) {
// printout("%s does not exist. Setting 100%% chance of single gamma-line with energy %g MeV\n",
// printoutf("%s does not exist. Setting 100%% chance of single gamma-line with energy %g MeV\n",
// filename, decay::nucdecayenergygamma(z, a) / EV / 1e6);
set_trivial_gamma_spectrum(nucindex);

Expand All @@ -140,7 +140,7 @@ void read_decaydata() {
assert_always(z != 28 || a != 57); // Ni-57 must have a gamma spectrum if present in list of nuclides
assert_always(z != 28 || a != 57); // Co-57 must have a gamma spectrum if present in list of nuclides
} else {
// printout("%s does not exist. No gamma decay from this nuclide.\n", filename);
// printoutf("%s does not exist. No gamma decay from this nuclide.\n", filename);
}
}

Expand All @@ -162,7 +162,7 @@ void init_gamma_linelist() {
for (int nucindex = 0; nucindex < decay::get_num_nuclides(); nucindex++) {
total_lines += std::ssize(gamma_spectra[nucindex]);
}
printout("total gamma-ray lines %td\n", total_lines);
printoutf("total gamma-ray lines %td\n", total_lines);

allnuc_gamma_line_list = std::vector<NucGammaLine>();
allnuc_gamma_line_list.reserve(total_lines);
Expand Down Expand Up @@ -205,7 +205,7 @@ void init_gamma_linelist() {

void init_xcom_photoion_data() {
// read the file
printout("reading XCOM photoionization data...\n");
printoutf("reading XCOM photoionization data...\n");
// reserve memory
for (int Z = 0; Z < numb_xcom_elements; Z++) {
photoion_data[Z].reserve(100);
Expand Down Expand Up @@ -247,7 +247,7 @@ __host__ __device__ auto choose_gamma_ray(const int nucindex) -> double {
}
}

printout("Failure to choose line (pellet_nucindex %d). Abort. zrand %g runtot %g\n", nucindex, zrand, runtot);
printoutf("Failure to choose line (pellet_nucindex %d). Abort. zrand %g runtot %g\n", nucindex, zrand, runtot);
assert_always(false);
}

Expand Down Expand Up @@ -296,24 +296,24 @@ auto choose_f(const double xx, const double zrand) -> double
int count = 0;
double err = 1e20;

// printout("new\n");
// printoutf("new\n");

double ftry = (f_max + f_min) / 2;
while ((err > 1.e-4) && (count < 1000)) {
ftry = (f_max + f_min) / 2;
const double sigma_try = sigma_compton_partial(xx, ftry);
// printout("ftry %g %g %g %g %g\n",ftry, f_min, f_max, try, norm);
// printoutf("ftry %g %g %g %g %g\n",ftry, f_min, f_max, try, norm);
if (sigma_try > norm) {
f_max = ftry;
err = (sigma_try - norm) / norm;
} else {
f_min = ftry;
err = (norm - sigma_try) / norm;
}
// printout("error %g\n",err);
// printoutf("error %g\n",err);
count++;
if (count == 1000) {
printout("Compton hit 1000 tries. %g %g %g %g %g\n", f_max, f_min, ftry, sigma_try, norm);
printoutf("Compton hit 1000 tries. %g %g %g %g %g\n", f_max, f_min, ftry, sigma_try, norm);
}
}

Expand Down Expand Up @@ -373,7 +373,7 @@ auto thomson_angle() -> double {

// handle physical Compton scattering event
void compton_scatter(Packet &pkt) {
// printout("Compton scattering.\n");
// printoutf("Compton scattering.\n");

const double xx = H * pkt.nu_cmf / ME / CLIGHT / CLIGHT;

Expand Down Expand Up @@ -605,7 +605,7 @@ auto get_chi_pair_prod_rf(const Packet &pkt) -> double {
double chi_rf = chi_cmf * calculate_doppler_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time);

if (chi_rf < 0) {
printout("Negative pair production sigma. Setting to zero. Abort? %g\n", chi_rf);
printoutf("Negative pair production sigma. Setting to zero. Abort? %g\n", chi_rf);
chi_rf = 0.;
}

Expand Down Expand Up @@ -730,19 +730,19 @@ void transport_gamma(Packet &pkt, const double t2) {
? globals::rmax * pkt.prop_time / globals::tmin
: 2 * globals::rmax * (pkt.prop_time + sdist / CLIGHT_PROP) / globals::tmin;
if (sdist > maxsdist) {
printout("Unreasonably large sdist (gamma). Abort. %g %g %g\n", globals::rmax, pkt.prop_time / globals::tmin,
sdist);
printoutf("Unreasonably large sdist (gamma). Abort. %g %g %g\n", globals::rmax, pkt.prop_time / globals::tmin,
sdist);
assert_always(false);
}

if (sdist < 0) {
printout("Negative distance (sdist). Abort?\n");
printoutf("Negative distance (sdist). Abort?\n");
sdist = 0;
}

if (((snext < 0) && (snext != -99)) || (snext >= grid::ngrid)) {
printout("Heading for inappropriate grid cell. Abort.\n");
printout("Current cell %d, target cell %d.\n", pkt.where, snext);
printoutf("Heading for inappropriate grid cell. Abort.\n");
printoutf("Current cell %d, target cell %d.\n", pkt.where, snext);
assert_always(false);
}

Expand Down Expand Up @@ -780,7 +780,7 @@ void transport_gamma(Packet &pkt, const double t2) {

assert_always(tdist >= 0);

// printout("sdist, tdist, edist %g %g %g\n",sdist, tdist, edist);
// printoutf("sdist, tdist, edist %g %g %g\n",sdist, tdist, edist);

if ((sdist < tdist) && (sdist < edist)) {
move_pkt_withtime(pkt, sdist / 2.);
Expand Down Expand Up @@ -976,8 +976,8 @@ void guttman_thermalisation(Packet &pkt) {
for (int i = 0; i < numb_rnd_dirs; i++) {
const double summand =
width * (1 - std::exp(-std::pow(t_gamma, 2.) / std::pow(t, 2.) * column_densities[i] / avg_column_density));
printout("width: %f t_gamma: %f t: %f column_densities[i]: %f avg_column_density: %f summand: %f", width, t_gamma,
t, column_densities[i], avg_column_density, summand);
printoutf("width: %f t_gamma: %f t: %f column_densities[i]: %f avg_column_density: %f summand: %f", width, t_gamma,
t, column_densities[i], avg_column_density, summand);
f_gamma += summand;
}
f_gamma /= (4 * PI);
Expand Down Expand Up @@ -1056,9 +1056,9 @@ __host__ __device__ void pellet_gamma_decay(Packet &pkt) {
}

pkt.pol_dir = vec_norm(pkt.pol_dir);
// printout("initialise pol state of packet %g, %g, %g, %g,
// printoutf("initialise pol state of packet %g, %g, %g, %g,
// %g\n",pkt.stokes_qu[0],pkt.stokes_qu[1],pkt.pol_dir[0],pkt.pol_dir[1],pkt.pol_dir[2]);
// printout("pkt direction %g, %g, %g\n",pkt.dir[0],pkt.dir[1],pkt.dir[2]);
// printoutf("pkt direction %g, %g, %g\n",pkt.dir[0],pkt.dir[1],pkt.dir[2]);
}

__host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) {
Expand Down

0 comments on commit faafca5

Please sign in to comment.