c++filterarmadilloastronomyssp

Computing apparent magnitude with Armadillo (C++)

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!

Solution

  • There are a few problems with this code and maybe a few operations that are not doing what you want.

    1. You are coping the the input vectors, which can be very costly and is completely unnecessary. Change
    mat compute_ab(vector <long double> waves,vector <long double> sed, vector <long double> fwaves,vector <long double> fresp){

    to

    mat compute_ab(const vector <long double>& waves,const vector <long double>& sed, const vector <long double>& fwaves,const vector <long double>& fresp){
    1. colvec and vec are the same thing. Column vectors of doubles. You can use just vec, if you want. You are converting the std::vector<double> to arma::colvec using conv_to<colvec>::from. This is fine, but again you are copying the elements. Since you just want to perform linear algebra operation with these std::vector<double> inputs and they are only used inside the function, you can do better.

    For instance, consider the code below

    std::vector<double> v;
v.push_back(2.4);
v.push_back(-1.4);
v.push_back(10);

arma::vec c = arma::conv_to<arma::colvec>::from(v);
arma::colvec n;

std::cout << &c[0] << std::endl;
std::cout << &v[0] << std::endl;

    The couts will print the memory address of the first element in the original vector and in the converted vec, which are indeed diferente. A better approach is using a special constructor from arma::vec that initializes the vector from a memory address and tell it to not copy the elements. That is, change the arma::vec c line to

    arma::vec c = arma::vec(v.data(), v.size(), false);

    The data member of std::vector gives a pointer to its elements that we can pass to arma::vec constructor along the number of elements. Now the couts will print the same address for &c[0] and &v[0] confirming that the elements are not being copied.

    1. Take a look in the output of your arma::conv calls. For instance, passing two input vectors each with 3 elements will result in a output vector with 5 elements. Look at the documentation of arma::conv. Particularly, there is a third argument called "shape" that might be useful to you.
    1. arma::trapz( X, Y ) documentation says

    Compute the trapezoidal integral of Y with respect to spacing in X. X must be a vector; its length must equal to the number of rows in Y.

    The dimensions of the inputs you are passing to arma::trapz don't match and you get the mentioned error. The reason is that the output of arma::conv is likely not what you expected and you probably want to pass "same" as the third argument of arma::conv.

    TLDR;

    My suggestion is that you print the result of each line in compute_ab to see if it is what you expect. If it is not, then look into that particular armadillo function's documentation. Do this and you will quickly get more confident using the armadillo library. Often a function in armadillo has several overloads.