Skip to content

Commit

Permalink
Merge pull request #33 from karllark/fix_rad_field_unc
Browse files Browse the repository at this point in the history 
Fix for missing rad field unc and adding nphotons in rad field calculation
  • Loading branch information
@karllark
karllark committed Jan 29, 2024
2 parents b2624f1 + b162e32 commit fe96ee1
Show file tree
Hide file tree
Showing 4 changed files with 139 additions and 68 deletions.
3 changes: 2 additions & 1 deletion DIRTY/get_sed_parameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@ void get_sed_parameters (ConfigFile& param_data,

unsigned int i = 0;
for (i = 0; i < runinfo.wavelength.size(); i++)
runinfo.sed_lum.push_back(1.);
// reasonable for a star, needed to get radiation field correct for single wavelength case
runinfo.sed_lum.push_back(1.e45);

} else {
#ifdef DEBUG_GSP
Expand Down
39 changes: 33 additions & 6 deletions DIRTY/output_model_grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,9 @@ void output_model_grid (geometry_struct& geometry,
string filename_unc = "!" + output.file_base;
filename_unc += "_rad_field_unc.fits";

string filename_npts = "!" + output.file_base;
filename_npts += "_rad_field_nphotons.fits";

// filename of the wavelength grid output file
string filename_wave = "!" + output.file_base;
filename_wave += "_wave_grid.fits";
Expand Down Expand Up @@ -46,6 +49,10 @@ void output_model_grid (geometry_struct& geometry,
fits_create_file(&out_unc_ptr,filename_unc.c_str(), &status);
check_fits_io(status, "fits_create_file : output_model_grid (rad_field unc)");

fitsfile *out_npts_ptr; // pointer to the output fits file
fits_create_file(&out_npts_ptr,filename_npts.c_str(), &status);
check_fits_io(status, "fits_create_file : output_model_grid (rad_field npts)");

fitsfile *out_tau_ptr; // pointer to the output fits file
fits_create_file(&out_tau_ptr,filename_tau.c_str(), &status);
check_fits_io(status, "fits_create_file : output_model_grid (tau)");
Expand Down Expand Up @@ -80,6 +87,10 @@ void output_model_grid (geometry_struct& geometry,
NumUtils::FourVector<float> tmp_rad_field_unc;
tmp_rad_field_unc.FVSize(geometry.grids[m].index_dim[0],geometry.grids[m].index_dim[1],geometry.grids[m].index_dim[2],n_waves);

// create a 4d matrix to copy the grid info into for output
NumUtils::FourVector<int> tmp_rad_field_npts;
tmp_rad_field_npts.FVSize(geometry.grids[m].index_dim[0],geometry.grids[m].index_dim[1],geometry.grids[m].index_dim[2],n_waves);

// create a 3d matrix to copy the grid info into for output
NumUtils::Cube<float> tmp_tau;
tmp_tau.CSize(geometry.grids[m].index_dim[0],geometry.grids[m].index_dim[1],geometry.grids[m].index_dim[2]);
Expand Down Expand Up @@ -119,7 +130,8 @@ void output_model_grid (geometry_struct& geometry,
tmp_num_H(i,j,k) = geometry.grids[m].grid(i,j,k).num_H;
for (n = 0; n < n_waves; n++) {
tmp_rad_field(i,j,k,n) = geometry.grids[m].grid(i,j,k).absorbed_energy[n];
if (geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n] >= 5) {
tmp_rad_field_npts(i,j,k,n) = geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n];
if (geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n] >= 1) {
rad_unc = geometry.grids[m].grid(i,j,k).absorbed_energy_x2[n]/geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n] -
pow(double(geometry.grids[m].grid(i,j,k).absorbed_energy[n]/geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n]),double(2.0));
if (rad_unc > 0.0)
Expand Down Expand Up @@ -165,6 +177,17 @@ void output_model_grid (geometry_struct& geometry,
cout << "done writing rad field unc" << endl;
#endif

// create and output each grid (rad_field npts)
fits_create_img(out_npts_ptr, 16, 4, four_vector_size, &status);
check_fits_io(status,"fits_create_image : output_model_grid (rad_field npts)");

fits_write_img(out_npts_ptr, TSHORT, 1, geometry.grids[m].index_dim[0]*geometry.grids[m].index_dim[1]*geometry.grids[m].index_dim[2]*n_waves,
&tmp_rad_field_npts[0], &status);

#ifdef DEBUG_OMG
cout << "done writing rad field npts" << endl;
#endif

// create and output each grid (tau)
fits_create_img(out_tau_ptr, -32, 3, geometry.grids[m].index_dim, &status);
check_fits_io(status,"fits_create_image : output_model_grid (tau)");
Expand Down Expand Up @@ -208,15 +231,15 @@ void output_model_grid (geometry_struct& geometry,
fits_write_comment(out_ptr, "**---------------------------------**",&status);
fits_write_comment(out_ptr, "Output of the DIRTY model",&status);
fits_write_comment(out_ptr, "Karl D. Gordon & Karl A. Misselt", &status);
fits_write_comment(out_ptr, "version v2.0prealpha (Oct 2009)", &status);
fits_write_comment(out_ptr, "version v2.2dev (Jan 2024)", &status);
fits_write_comment(out_ptr, "**---------------------------------**",&status);
check_fits_io(status,"fits_write_comment : output_model_grid");

// final stuff for primary header
fits_write_comment(out_unc_ptr, "**---------------------------------**",&status);
fits_write_comment(out_unc_ptr, "Output of the DIRTY model",&status);
fits_write_comment(out_unc_ptr, "Karl D. Gordon & Karl A. Misselt", &status);
fits_write_comment(out_unc_ptr, "version v2.0prealpha (Oct 2009)", &status);
fits_write_comment(out_unc_ptr, "vversion v2.2dev (Jan 2024)", &status);
fits_write_comment(out_unc_ptr, "**---------------------------------**",&status);
check_fits_io(status,"fits_write_comment : output_model_grid (unc)");

Expand All @@ -228,23 +251,23 @@ void output_model_grid (geometry_struct& geometry,
fits_write_comment(out_tau_ptr, "**---------------------------------**",&status);
fits_write_comment(out_tau_ptr, "Output of the DIRTY model",&status);
fits_write_comment(out_tau_ptr, "Karl D. Gordon & Karl A. Misselt", &status);
fits_write_comment(out_tau_ptr, "version v2.0prealpha (Oct 2009)", &status);
fits_write_comment(out_tau_ptr, "version v2.2dev (Jan 2024)", &status);
fits_write_comment(out_tau_ptr, "**---------------------------------**",&status);
check_fits_io(status,"fits_write_comment : output_model_grid (tau)");

// final stuff for primary header
fits_write_comment(out_num_H_ptr, "**---------------------------------**",&status);
fits_write_comment(out_num_H_ptr, "Output of the DIRTY model",&status);
fits_write_comment(out_num_H_ptr, "Karl D. Gordon & Karl A. Misselt", &status);
fits_write_comment(out_num_H_ptr, "version v2.0prealpha (Oct 2009)", &status);
fits_write_comment(out_num_H_ptr, "version v2.2dev (Jan 2024)", &status);
fits_write_comment(out_num_H_ptr, "**---------------------------------**",&status);
check_fits_io(status,"fits_write_comment : output_model_grid (num_H)");

// final stuff for primary header
fits_write_comment(out_pos_ptr, "**---------------------------------**",&status);
fits_write_comment(out_pos_ptr, "Output of the DIRTY model",&status);
fits_write_comment(out_pos_ptr, "Karl D. Gordon & Karl A. Misselt", &status);
fits_write_comment(out_pos_ptr, "version v2.0prealpha (Oct 2009)", &status);
fits_write_comment(out_pos_ptr, "version v2.2dev (Jan 2024)", &status);
fits_write_comment(out_pos_ptr, "**---------------------------------**",&status);
check_fits_io(status,"fits_write_comment : output_model_grid (pos)");
}
Expand All @@ -267,6 +290,10 @@ void output_model_grid (geometry_struct& geometry,
fits_close_file(out_unc_ptr, &status);
check_fits_io(status,"fits_close_file : output_model_grid (unc)");

// close FITS File
fits_close_file(out_npts_ptr, &status);
check_fits_io(status,"fits_close_file : output_model_grid (npts)");

// close FITS File
fits_close_file(out_wave_ptr, &status);
check_fits_io(status,"fits_close_file : output_model_grid (wave)");
Expand Down
156 changes: 95 additions & 61 deletions DIRTY/setup_absorbed_energy_grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,75 +6,109 @@
// ======================================================================
#include "setup_absorbed_energy_grid.h"

void setup_absorbed_energy_grid (geometry_struct& geometry,
runinfo_struct& runinfo,
int doing_emission)
void setup_absorbed_energy_grid(geometry_struct& geometry,
runinfo_struct& runinfo, int doing_emission)

{
if ((!geometry.abs_energy_grid_initialized) && (!doing_emission))
runinfo.absorbed_energy.resize(runinfo.n_waves,0.0);
runinfo.absorbed_energy.resize(runinfo.n_waves, 0.0);

int i,j,k,m,n = 0;
int i, j, k, m, n = 0;
// loop over all the defined grids
for (m = 0; m < int(geometry.grids.size()); m++) {
// loop of the cells in this grid
for (k = 0; k < geometry.grids[m].index_dim[2]; k++)
for (j = 0; j < geometry.grids[m].index_dim[1]; j++)
for (i = 0; i < geometry.grids[m].index_dim[0]; i++) {
if (doing_emission) {
// save the existing value of the radiation field (needed for dust or ere emission)
// but only on the first iteration (i.e., save the stellar RT radiation field)
for (n = 0; n < runinfo.n_waves; n++) {
if (doing_emission == 1) {
geometry.grids[m].grid(i,j,k).save_radiation_field_density[n] =
geometry.grids[m].grid(i,j,k).absorbed_energy[n];
geometry.grids[m].grid(i,j,k).save_radiation_field_density_x2[n] =
geometry.grids[m].grid(i,j,k).absorbed_energy_x2[n];
geometry.grids[m].grid(i,j,k).save_radiation_field_density_num_photons[n] =
geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n];
}
// if ((k == 5) && (j == 5) && (i == 5)) {
// cout << n << " " << geometry.grids[m].grid(i,j,k).save_radiation_field_density[n] << endl;
// }
// zero out the current absorbed energy value
geometry.grids[m].grid(i,j,k).absorbed_energy[n] = 0.0;
geometry.grids[m].grid(i,j,k).absorbed_energy_x2[n] = 0.0;
geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n] = 0;
}
} else {
// if this is the first time or if memory storage requested
// then push zero into the absorbed_energy variable in each grid cell
if (geometry.abs_energy_storage_type == 1) {// disk memory
cout << "need new code for disk storage of absorbed energy grid" << endl;
cout << "ever since limited polychromaticism was implemented (27 Mar 2009 - KDG)" << endl;
exit(8);
// if (!geometry.abs_energy_grid_initialized) {
// geometry.grids[m].grid(i,j,k).absorbed_energy.push_back(0.0);
// geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons.push_back(0);
// geometry.grids[m].grid(i,j,k).save_radiation_field_density.push_back(0.0);
// geometry.grids[m].grid(i,j,k).save_radiation_field_density_num_photons.push_back(0);
// } else {
// geometry.grids[m].grid(i,j,k).absorbed_energy[0] = 0.0;
// geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[0] = 0;
// }
} else
if (!geometry.abs_energy_grid_initialized) {
geometry.grids[m].grid(i,j,k).absorbed_energy.resize(runinfo.n_waves,0.0);
geometry.grids[m].grid(i,j,k).absorbed_energy_x2.resize(runinfo.n_waves,0.0);
geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons.resize(runinfo.n_waves,0);
geometry.grids[m].grid(i,j,k).save_radiation_field_density.resize(runinfo.n_waves,0.0);
geometry.grids[m].grid(i,j,k).save_radiation_field_density_x2.resize(runinfo.n_waves,0.0);
geometry.grids[m].grid(i,j,k).save_radiation_field_density_num_photons.resize(runinfo.n_waves,0);
} else {
for (n = 0; n < runinfo.n_waves; n++) {
geometry.grids[m].grid(i,j,k).absorbed_energy[n] = 0.0;
geometry.grids[m].grid(i,j,k).absorbed_energy_x2[n] = 0.0;
geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n] = 0;
}
}
}
}
for (i = 0; i < geometry.grids[m].index_dim[0]; i++) {
if (doing_emission) {
// save the existing value of the radiation field (needed for dust
// or ere emission) but only on the first iteration (i.e., save the
// stellar RT radiation field)
for (n = 0; n < runinfo.n_waves; n++) {
if (doing_emission == 1) {
geometry.grids[m]
.grid(i, j, k)
.save_radiation_field_density[n] =
geometry.grids[m].grid(i, j, k).absorbed_energy[n];
geometry.grids[m]
.grid(i, j, k)
.save_radiation_field_density_x2[n] =
geometry.grids[m].grid(i, j, k).absorbed_energy_x2[n];
geometry.grids[m]
.grid(i, j, k)
.save_radiation_field_density_num_photons[n] =
geometry.grids[m]
.grid(i, j, k)
.absorbed_energy_num_photons[n];
}
// if ((k == 5) && (j == 5) && (i == 5)) {
// cout << n << " " <<
// geometry.grids[m].grid(i,j,k).save_radiation_field_density[n]
// << endl;
// }
// zero out the current absorbed energy value
geometry.grids[m].grid(i, j, k).absorbed_energy[n] = 0.0;
geometry.grids[m].grid(i, j, k).absorbed_energy_x2[n] = 0.0;
geometry.grids[m].grid(i, j, k).absorbed_energy_num_photons[n] =
0;
}
} else {
// if this is the first time or if memory storage requested
// then push zero into the absorbed_energy variable in each grid
// cell
if (geometry.abs_energy_storage_type == 1) { // disk memory
cout << "need new code for disk storage of absorbed energy grid"
<< endl;
cout << "ever since limited polychromaticism was implemented (27 "
"Mar 2009 - KDG)"
<< endl;
exit(8);
// if (!geometry.abs_energy_grid_initialized) {
// geometry.grids[m].grid(i,j,k).absorbed_energy.push_back(0.0);
// geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons.push_back(0);
// geometry.grids[m].grid(i,j,k).save_radiation_field_density.push_back(0.0);
// geometry.grids[m].grid(i,j,k).save_radiation_field_density_num_photons.push_back(0);
// } else {
// geometry.grids[m].grid(i,j,k).absorbed_energy[0]
// = 0.0;
// geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[0]
// = 0;
// }
} else if (!geometry.abs_energy_grid_initialized) {
geometry.grids[m].grid(i, j, k).absorbed_energy.resize(
runinfo.n_waves, 0.0);
geometry.grids[m].grid(i, j, k).absorbed_energy_x2.resize(
runinfo.n_waves, 0.0);
geometry.grids[m]
.grid(i, j, k)
.absorbed_energy_num_photons.resize(runinfo.n_waves, 0);

if (doing_emission == 1) {
geometry.grids[m]
.grid(i, j, k)
.save_radiation_field_density.resize(runinfo.n_waves, 0.0);
geometry.grids[m]
.grid(i, j, k)
.save_radiation_field_density_x2.resize(runinfo.n_waves,
0.0);
geometry.grids[m]
.grid(i, j, k)
.save_radiation_field_density_num_photons.resize(
runinfo.n_waves, 0);
}
} else {
for (n = 0; n < runinfo.n_waves; n++) {
geometry.grids[m].grid(i, j, k).absorbed_energy[n] = 0.0;
geometry.grids[m].grid(i, j, k).absorbed_energy_x2[n] = 0.0;
geometry.grids[m].grid(i, j, k).absorbed_energy_num_photons[n] =
0;
}
}
}
}
}
// set to know that the absorbed energy grid has been initialized the first time
if (!geometry.abs_energy_grid_initialized) geometry.abs_energy_grid_initialized = 1;
// set to know that the absorbed energy grid has been initialized the first
// time
if (!geometry.abs_energy_grid_initialized)
geometry.abs_energy_grid_initialized = 1;
}
9 changes: 9 additions & 0 deletions DIRTY/store_absorbed_energy_grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,15 @@ void store_absorbed_energy_grid (geometry_struct& geometry,
j_temp /= flux_mult_factor;
j_temp_x2 /= (flux_mult_factor*flux_mult_factor);

//if (j_temp > 0) {
// cout << "emitted_lum = " << runinfo.sed_lum[index] << "; ";
// cout << "num_H = " << geometry.grids[m].grid(i,j,k).num_H << "; ";
// cout << "C_abs = " << runinfo.ave_C_abs[index] << "; ";
// cout << "j_temp" << " ";
// cout << j_temp << " ";
// cout << j_temp_x2 << endl;
//}

// check for roundoff error before converting to double
if (j_temp < 1e-38) {
#ifdef DEBUG_SAEG
Expand Down

0 comments on commit fe96ee1

Please sign in to comment.