Computing apparent magnitude with Armadillo (C++)

  armadillo, astronomy, c++, filter, ssp

I’m looking for someone very expert in using Armadillo (C++). I am a wannabe astronomer and I’m trying to compute the ab apparent magnitude m_ab of a SSP (simple stellar population, a group of chemically homogeneus and coeval stars) in a specific filter of a telescope.

The input of the function compute_ab, which does this calculation, are the wavelength wavelength and the corresponding Spectral Energy Distribution sed of the SSP (basically it is the luminoisty of the SSP per unit wavelength, over a range of wavelength); the input fresp and fwaves are the throughput of the telescope (basically the filter in a certain band) and the corresponding wavelength range over which the througput is distributed. They are std::vector<long double> in 1D. What I need to output is a number, possibly a double or long double.

What I surely need to compute m_ab from these information is an interpolation, because the SED and the throughput can be at very different wavelengths; it is a convolution and an integration. Physically speaking, the passages I make in this function are correct, so I’m asking some help to use Armadillo, am I doing it correctly? How can I set the output as a double? Moreover, I’m getting this error right now, when I run my code:

    error: trapz(): length of X must equal the number of rows in Y when dim=0
terminate called after throwing an instance of 'std::logic_error'
  what():  trapz(): length of X must equal the number of rows in Y when dim=0

Here it is the function:

mat compute_ab(vector <long double> waves,vector <long double> sed, vector <long double> fwaves,vector <long double> fresp){

  colvec filter_interp;

  colvec csed = conv_to< colvec >::from(sed);
  colvec cwaves = conv_to< colvec >::from(waves);
  colvec cfwaves = conv_to< colvec >::from(fwaves);
  colvec cfresp = conv_to< colvec >::from(fresp);

  arma::interp1(cfwaves,cfresp,cwaves,filter_interp);

  colvec filterSpec = arma::conv(filter_interp,csed);

  colvec i1 = arma::conv(filterSpec, cwaves);
  colvec i2 = arma::conv(filter_interp, 1./cwaves); 

  mat I1=arma::trapz(i1,cwaves);
  mat I2=arma::trapz(i2,cwaves);

  mat fnu=I1/I2/speedcunitas;
  mat mAB= -2.5 * log10(fnu) -48.6;
  return mAB;

}

Thank you all for your help!

Source: Windows Questions C++

LEAVE A COMMENT