generate_riaf_image.cpp File Reference

Example illustrating the creationg of a RIAF model. More...

#include "model_image_sed_fitted_riaf.h"
#include <mpi.h>

Include dependency graph for generate_riaf_image.cpp:

Functions
int	main (int argc, char *argv[])

Detailed Description

Author

Avery E. Broderick

Date

July, 2017

Themis allows a vast variety of models to be compared to EHT data. This example shows how to query the SED-fitted RIAF model (an explicit model_image class), and generate an intrinsic image for the best fit parameters in Broderick et al. (2016). The output file, riaf_image.d, is in the standard pmap format and may be plotted with the python script vrt2_image.py via:

$ python ../../analysis/plot_vrt2_images.py riaf_image.d

The resulting image is shown below.

Image generated by plot_vrt2_image.py from the output pmap file.

// Initialize MPI
MPI_Init(&argc, &argv);

// Set the model to the SED-fitted RIAF class.  To do this, read in the
// file containing the disk parameter SED-fit results, located in
// ../../src/VRT2/DataFiles/2010_combined_fit_parameters.d
Themis::model_image_sed_fitted_riaf image("../../src/VRT2/DataFiles/2010_combined_fit_parameters.d");

// Make a low-resolution image if that is more convenient.  Commenting this
// line out makes the default resolution images.
// image.use_small_images();

// Choose a specific set of physical parameters at which to create an image.
std::vector<double> parameters;
parameters.push_back( 0.15 ); // Spin parameter (0-1)
parameters.push_back( std::cos(60*M_PI/180.0) ); // Cos(Inclination)
parameters.push_back( 0.0 ); // Position angle (which doesn't impact the intrinsic image).

// Generate a model.
image.generate_model(parameters);

// Access the image data (not usually required)
std::vector<std::vector<double> > alpha, beta, I;
image.get_image(alpha,beta,I);

// Convert from radians to projected GM/c^2 just for convenience
double rad2M = VRT2::VRT2_Constants::D_SgrA_cm/VRT2::VRT2_Constants::M_SgrA_cm;

// Return the image (not usually required) and output it to a file in the
// standard pmap format -- a rastered array of ASCII values with the limits
// and dimensions at the top.
std::ofstream imout("riaf_image.d");
imout << "Nx:"
  << std::setw(15) << rad2M*alpha[0][0]
  << std::setw(15) << rad2M*alpha[alpha.size()-1].back()
  << std::setw(15) << alpha.size()
  << '\n';
imout << "Ny:"
  << std::setw(15) << rad2M*beta[0][0]
  << std::setw(15) << rad2M*beta[beta.size()-1].back()
  << std::setw(15) << beta.size()
  << '\n';
imout << std::setw(5) << "i"
  << std::setw(5) << "j"
  << std::setw(15) << "I (Jy/sr)"
  << '\n';
for (size_t j=0; j<alpha[0].size(); ++j)
  for (size_t i=0; i<alpha.size(); ++i)
    imout << std::setw(5) << i
      << std::setw(5) << j
      << std::setw(15) << I[i][j]
      << '\n';

// Finalize MPI
MPI_Finalize();

return 0;