使用犰狳 (C++) 计算视星等

Question

我正在寻找在 telescope 的特定过滤器中使用 Armadillo (C++). I am a wannabe astronomer and I'm trying to compute the ab apparent magnitude m_ab SSP（简单星族，一组化学同质和同时代的恒星）的专家。

进行此计算的函数 compute_ab 的输入是波长 wavelength 和 SSP 的相应光谱能量分布 sed（基本上是在一定波长范围内每单位波长的 SSP）；输入 fresp 和 fwaves 是 telescope 的吞吐量（基本上是某个波段的滤波器）以及吞吐量分布的相应波长范围。它们在 1D 中是 std::vector<long double>。我需要输出的是一个数字，可能是double或者long double。

我肯定需要根据这些信息计算 m_ab 的是插值，因为 SED 和吞吐量可能处于非常不同的波长；它是一个卷积和一个集成。从物理上讲，我在这个函数中所做的段落是正确的，所以我正在寻求一些帮助来使用犰狳，我做得对吗？如何将输出设置为双精度？此外，当我运行我的代码：

时，我现在收到此错误

    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

这里是函数：

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;

}

谢谢大家的帮助！

Answer 1

这段代码有一些问题，可能有一些操作没有按照您的要求进行。

您正在处理输入向量，这可能非常昂贵并且完全没有必要。变化

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

到

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

colvec 和 vec 是一回事。双精度的列向量。如果需要，您可以只使用 vec。您正在使用 conv_to<colvec>::from 将 std::vector<double> 转换为 arma::colvec。这很好，但是您又在复制元素。由于您只想使用这些 std::vector<double> 输入执行线性代数运算，并且它们仅在函数内部使用，因此您可以做得更好。

例如，考虑下面的代码

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;

cout会打印原始vector和转换后的vec中第一个元素的内存地址，这确实是不同的。更好的方法是使用来自 arma::vec 的特殊构造函数，它从内存地址初始化向量并告诉它不要复制元素。即，将 arma::vec c 行更改为

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

std::vector 的 data member 给出了一个指向它的元素的指针，我们可以沿着元素的数量传递给 arma::vec 构造函数。现在 cout 将为 &c[0] 和 &v[0] 打印相同的地址，确认元素未被复制。

查看 arma::conv 调用的输出。例如，传递两个每个具有 3 个元素的输入向量将导致具有 5 个元素的输出向量。查看 arma::conv 的文档。特别是，还有第三个参数称为“形状”，可能对您有用。

arma::trapz( X, Y ) 文档说

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.

您传递给 arma::trapz 的输入的尺寸不匹配，您会收到上述错误。原因是 arma::conv 的输出可能不是您所期望的，您可能希望传递“相同”作为 arma::conv.

的第三个参数

TLDR;

我的建议是你把compute_ab中每一行的结果都打印出来，看看是不是你期望的那样。如果不是，则查看该特定犰狳函数的文档。这样做，您将很快对使用 armadillo 库更有信心。犰狳中的一个函数通常有多个重载。

使用犰狳 (C++) 计算视星等

Computing apparent magnitude with Armadillo (C++)

c++

ssp

filter

astronomy

armadillo