Merge pull request #35 from karllark/update_rad_field

Output Radiation Field updates

DIRTY/forced_first_scatter.cpp

@@ -200,8 +200,19 @@ int forced_first_scatter (geometry_struct& geometry,
 	// deposit the energy
 	grid_cell& this_cell = geometry.grids[dummy_photon.path_pos_index[0][i]].grid(dummy_photon.path_pos_index[1][i],dummy_photon.path_pos_index[2][i],dummy_photon.path_pos_index[3][i]);
 	this_cell.absorbed_energy[geometry.abs_energy_wave_index] += abs_weight;
-	this_cell.absorbed_energy_x2[geometry.abs_energy_wave_index] += abs_weight*abs_weight;
-	this_cell.absorbed_energy_num_photons[geometry.abs_energy_wave_index]++;
+  if (this_cell.last_photon_number != photon.number) {
+    this_cell.absorbed_energy_num_photons[geometry.abs_energy_wave_index]++;
+	  this_cell.absorbed_energy_x2[geometry.abs_energy_wave_index] += abs_weight*abs_weight;
+    this_cell.last_photon_number = photon.number;
+    this_cell.last_photon_absorbed_energy += abs_weight;
+  } else {
+    // compute the difference in x2 values previously added to what should be added
+    // avoids having to save this information separately and then have a special step at the end to add the total photon's contribution
+    float prev_x2 = pow(this_cell.last_photon_absorbed_energy, 2.0);
+    this_cell.last_photon_absorbed_energy += abs_weight;
+    // add the difference = just like adding the correct x2
+    this_cell.absorbed_energy_x2[geometry.abs_energy_wave_index] += pow(this_cell.last_photon_absorbed_energy, 2.0) - prev_x2;
+  }
 
 	// move to the next grid cell
 	tau_entering = tau_leaving;

DIRTY/include/grid_cell.h

@@ -20,6 +20,8 @@ struct grid_cell {
                                    // can handle multiwavelengths if desired
   vector<float> absorbed_energy_x2;  // sum of the squared absorbed energy - allows uncertainties
   vector<int> absorbed_energy_num_photons;  // number of photons contributing to the absorbed energy in this cell
+  int last_photon_number;  // number of last photon to deposit energy in this cell
+  float last_photon_absorbed_energy;   // energy of last photon, needed to appropriately sum energy once for same photon
   vector<float> save_radiation_field_density;    // save the existing radiation field density (needed for the self-absorption calculation)
   vector<float> save_radiation_field_density_x2;    // ditto
   vector<int> save_radiation_field_density_num_photons;  // number of photons contributing to the absorbed energy in this cell

DIRTY/output_model_grid.cpp

@@ -88,7 +88,7 @@ void output_model_grid (geometry_struct& geometry,
     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<short int> tmp_rad_field_npts;
+    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
@@ -131,15 +131,31 @@ void output_model_grid (geometry_struct& geometry,
 	  			for (n = 0; n < n_waves; n++) {
 	    			tmp_rad_field(i,j,k,n) = geometry.grids[m].grid(i,j,k).absorbed_energy[n];
             tmp_rad_field_npts(i,j,k,n) = geometry.grids[m].grid(i,j,k).absorbed_energy_num_photons[n];
+            // compute the uncertainty on the average contribution from an individual photon
 	    			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)
-								rad_unc = sqrt(rad_unc);
-	      			else
-								rad_unc = 0.0;
+              // using eq. 14 of Camps & Baes (2018), equivalent to eqns in Gordon et al. (2001) with fewer computations
+              // *except* that the number of photons to use is the total, not the number in the cell
+              // do not understand why this is the case, but empirical tests with many independent runs confirm this (KDG 13 May 2024)
+              if (geometry.grids[m].grid(i,j,k).absorbed_energy[n] > 0.0) {
+                rad_unc = geometry.grids[m].grid(i,j,k).absorbed_energy_x2[n]/pow(double(geometry.grids[m].grid(i,j,k).absorbed_energy[n]),double(2.0));
+                rad_unc = sqrt(rad_unc - (1./output.outputs[0].total_num_photons));
+              } else
+                rad_unc = 0.0;
+
+              // scale to the uncertainty on the radiation field
+              rad_unc *= tmp_rad_field(i,j,k,n);
+              // store the result
 	      			tmp_rad_field_unc(i,j,k,n) = rad_unc;
-	    			}
+              // convert from radiation field mean intensity J (needed for dust emission calculations)
+              // to radiation field density U
+              tmp_rad_field(i,j,k,n) *= Constant::U_J;
+              tmp_rad_field_unc(i,j,k,n) *= Constant::U_J;
+	    			} else
+              // store the indexes of the subgrids for cells that are subdivided
+              if (tmp_tau(i,j,k) < 0.0) {
+                tmp_rad_field(i,j,k,n) = tmp_tau(i,j,k);
+                tmp_rad_field_unc(i,j,k,n) = tmp_tau(i,j,k);
+              }
 	  			}
 				}
 
@@ -178,10 +194,10 @@ void output_model_grid (geometry_struct& geometry,
 #endif
 
     // create and output each grid (rad_field npts)
-    fits_create_img(out_npts_ptr, 16, 4, four_vector_size, &status);
+    fits_create_img(out_npts_ptr, 32, 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,
+    fits_write_img(out_npts_ptr, TINT, 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

DIRTY/scatter_photon.cpp

@@ -191,8 +191,19 @@ void scatter_photon (geometry_struct& geometry,
       // deposit the energy
       grid_cell& this_cell = geometry.grids[dummy_photon.path_pos_index[0][i]].grid(dummy_photon.path_pos_index[1][i],dummy_photon.path_pos_index[2][i],dummy_photon.path_pos_index[3][i]);
       this_cell.absorbed_energy[geometry.abs_energy_wave_index] += abs_weight;
-      this_cell.absorbed_energy_x2[geometry.abs_energy_wave_index] += abs_weight*abs_weight;
-      this_cell.absorbed_energy_num_photons[geometry.abs_energy_wave_index]++;
+      if (this_cell.last_photon_number != photon.number) {
+        this_cell.absorbed_energy_num_photons[geometry.abs_energy_wave_index]++;
+        this_cell.absorbed_energy_x2[geometry.abs_energy_wave_index] += abs_weight*abs_weight;
+        this_cell.last_photon_number = photon.number;
+        this_cell.last_photon_absorbed_energy = abs_weight;
+      } else {
+        // compute the difference in x2 values previously added to what should be added
+        // avoids having to save this information separately and then have a special step at the end to add the total photon's contribution
+      float prev_x2 = pow(this_cell.last_photon_absorbed_energy, 2.0);
+      this_cell.last_photon_absorbed_energy += abs_weight;
+      // add the difference = just like adding the correct x2
+      this_cell.absorbed_energy_x2[geometry.abs_energy_wave_index] += pow(this_cell.last_photon_absorbed_energy, 2.0) - prev_x2;
+      }
 
       // move to the next grid cell
       tau_entering = tau_leaving;

DIRTY/setup_absorbed_energy_grid.cpp

@@ -82,6 +82,8 @@ void setup_absorbed_energy_grid(geometry_struct& geometry,
               geometry.grids[m]
                   .grid(i, j, k)
                   .absorbed_energy_num_photons.resize(runinfo.n_waves, 0);
+              geometry.grids[m].grid(i, j, k).last_photon_number = -1;
+              geometry.grids[m].grid(i, j, k).last_photon_absorbed_energy = 0.0;
 
               if (doing_emission == 1) {
                 geometry.grids[m]

DIRTY/store_absorbed_energy_grid.cpp

@@ -154,6 +154,9 @@ void store_absorbed_energy_grid (geometry_struct& geometry,
 // 	    exit(8);
 	    geometry.grids[m].grid(i,j,k).absorbed_energy[geometry.abs_energy_wave_index] = float(j_temp);
 	    geometry.grids[m].grid(i,j,k).absorbed_energy_x2[geometry.abs_energy_wave_index] = float(j_temp_x2);
+		// reset the number of photon to last contribute to this cell
+		geometry.grids[m].grid(i,j,k).last_photon_number = -1;
+		geometry.grids[m].grid(i,j,k).last_photon_absorbed_energy = 0.0;
 
 #ifdef DEBUG_SAEG
 	    cout << "J_temp = " << j_temp << endl;