Skip to content

Commit

Permalink
Merge pull request #135 from riclarsson/master
Browse files Browse the repository at this point in the history 
Work towards ECS model
  • Loading branch information
@olemke
olemke committed Mar 2, 2020
2 parents 4ea92b7 + 6fe3c42 commit 08de9d4
Show file tree
Hide file tree
Showing 11 changed files with 254 additions and 82 deletions.
12 changes: 5 additions & 7 deletions src/absorptionlines.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -503,15 +503,13 @@ Absorption::SingleLineExternal Absorption::ReadFromArtscat4Stream(istream& is) {
if (q[1] > -1)
data.quantumidentity.UpperQuantumNumbers().Set(QuantumNumberType::J, Rational(q[1]));
if (q[2] > -1)
data.quantumidentity.UpperQuantumNumbers().Set(QuantumNumberType::F,
q[2] - Rational(1, 2));
data.quantumidentity.UpperQuantumNumbers().Set(QuantumNumberType::F, q[2] - 1_2);
if (q[3] > -1)
data.quantumidentity.LowerQuantumNumbers().Set(QuantumNumberType::N, Rational(q[3]));
if (q[4] > -1)
data.quantumidentity.LowerQuantumNumbers().Set(QuantumNumberType::J, Rational(q[4]));
if (q[5] > -1)
data.quantumidentity.LowerQuantumNumbers().Set(QuantumNumberType::F,
q[5] - Rational(1, 2));
data.quantumidentity.LowerQuantumNumbers().Set(QuantumNumberType::F, q[5] - 1_2);
}
}
}
Expand Down Expand Up @@ -2955,17 +2953,17 @@ bool Absorption::line_is_id(const Absorption::Lines& band, const QuantumIdentifi

Numeric Absorption::reduced_rovibrational_dipole(Rational Jf, Rational Ji, Rational lf, Rational li, Rational k) {
const Numeric val = sqrt(2 * Jf + 1) * wigner3j(Jf, k, Ji, li, lf - li, -lf);
if ((Jf + lf + 1) % 2)
if (not even(Jf + lf + 1))
return -val;
else
return +val;
}

Numeric Absorption::reduced_magnetic_quadrapole(Rational Jf, Rational Ji, Rational N) {
constexpr Rational one(1, 1);
constexpr Rational one=1;
const Numeric val = sqrt(6 * (2 * Jf + 1) * (2 * Ji + 1)) *
wigner6j(one, one, one, Ji, Jf, N);
if ((Jf + N) % 2)
if (not even(Jf + N))
return -val;
else
return +val;
Expand Down
32 changes: 16 additions & 16 deletions src/lin_alg.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2384,22 +2384,11 @@ void linreg(Vector& p, ConstVectorView x, ConstVectorView y) {
p[0] = s3 / Numeric(n) - p[1] * xm;
}

/*!
* Least squares fitting by solving x for known A and y
*
* (A^T A)x = A^T y
*
* Returns the squared residual, i.e., <(A^T A)x-A^T y|(A^T A)x-A^T y>.
*
* \param x Out: As equation
* \param A In: As equation
* \param y In: As equation
*/
Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y) noexcept {
Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y, bool residual) noexcept {
// Size of the problem
const Index n = x.nelem();
Matrix AT, ATA(n, n);
Vector ATy(n), r(n);
Vector ATy(n);

// Solver
AT = transpose(A);
Expand All @@ -2408,7 +2397,18 @@ Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y) noexcept {
solve(x, ATA, ATy);

// Residual
mult(r, ATA, x);
r -= ATy;
return r * r;
if (residual) {
Vector r(n);
mult(r, ATA, x);
r -= ATy;
return r * r;
} else {
return 0;
}
}

Eigen::ComplexEigenSolver<Eigen::MatrixXcd> eig(const Eigen::Ref<Eigen::MatrixXcd> A)
{
Eigen::ComplexEigenSolver<Eigen::MatrixXcd> ces;
return ces.compute(A);
}
23 changes: 22 additions & 1 deletion src/lin_alg.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,27 @@ Numeric det(ConstMatrixView A);

void linreg(Vector& p, ConstVectorView x, ConstVectorView y);

Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y) noexcept;

/** Least squares fitting by solving x for known A and y
*
* (A^T A)x = A^T y
*
* Returns the squared residual, i.e., <(A^T A)x-A^T y|(A^T A)x-A^T y>.
*
* @param[in] x As equation
* @param[in] A As equation
* @param[in] y As equation
* @param[in] residual (optional) Returns the residual if true
* @return Squared residual or 0
*/
Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y, bool residual=true) noexcept;


/** Return the Eigen decomposition of the eigen matrix
*
* @param[in] A a matrix to eigen value decompose
* @return Object with eigenvalues and eigenvectors computed
*/
Eigen::ComplexEigenSolver<Eigen::MatrixXcd> eig(const Eigen::Ref<Eigen::MatrixXcd> A);

#endif // linalg_h
7 changes: 3 additions & 4 deletions src/linescaling.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,10 +173,9 @@ Numeric dboltzman_ratio_dT(const Numeric& boltzmann_ratio,
}

// Boltzmann factor at T
Numeric boltzman_factor(Numeric T, Numeric E0) {
using namespace Constant;
static constexpr Numeric c1 = -1 / k;
return std::exp(c1 * E0 / T);
Numeric boltzman_factor(Numeric T, Numeric E0)
{
return std::exp(- E0 / (Constant::k*T));
}

// Boltzmann factor at T
Expand Down
148 changes: 125 additions & 23 deletions src/predefined_absorption_models.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,13 @@
*/

#include "predefined_absorption_models.h"
#include "lin_alg.h"
#include "linescaling.h"
#include "wigner_functions.h"


constexpr auto nlines_mpm2020 = 38 + 6;
constexpr auto necs2020 = 38;
constexpr auto nlines_mpm2020 = necs2020 + 6;


constexpr LineShape::SingleSpeciesModel init_mpm2020_slsm(Numeric g00,
Expand Down Expand Up @@ -504,53 +507,152 @@ void Absorption::PredefinedModel::makarov2020_o2_lines_mpm(Matrix& xsec,
}


void Absorption::PredefinedModel::makarov2020_o2_lines_ecs(Numeric P, Numeric T, Numeric water_vmr)
void normalize_relaxation_matrix(Eigen::Ref<Eigen::MatrixXcd> W,
const Eigen::Ref<Eigen::ArrayXd> rho,
const Eigen::Ref<Eigen::ArrayXd> d)
{
constexpr auto necs2020 = nlines_mpm2020 - 6;
// FIXME:
// rho[i] / rho[j] or rho[j] / rho[i] in the next for-loop. [See {Prop. 1}]
// The current order is the same as for the second-to-next for-loop. [See {Prop. 2}]
// If the first loop changes, the renormalization in the second loop
// must change...

// Population density balanced parameters
for (Index i=0; i<necs2020; i++) {
for (Index j=0; j<i; j++) {
if (i not_eq j) {
W(i, j) = W(j, i) * (rho[i] / rho[j]); // {Prop. 1} --- division-order should change?
}
}
}

// Balance by the reduced dipole
for (Index i=0; i<necs2020; i++) {
auto norm = - std::imag(d[i] * W(i, i)) / std::imag(W.row(i) * d.abs().matrix() - std::abs(d[i]) * W(i, i)); // {Prop. 2} --- dot-product axis should to column?
for (Index j=0; j<necs2020; j++) {
if (i not_eq j) {
W(i, j) *= norm; // {Prop. 2} --- If yes, then change W(i, j) to W(j, i)
}
}
}
}


void Absorption::PredefinedModel::makarov2020_o2_lines_ecs(ComplexVector& I, const Vector& f, Numeric P, Numeric T, Numeric water_vmr)
{
using Constant::h;
using Constant::k;
using Constant::pow3;

constexpr auto qids = init_mpm2020_qids(0, 0);
constexpr auto lsm = init_mpm2020_lsm();
constexpr auto T0 = 300;

Vector dl(necs2020, 9);
Matrix W(necs2020, necs2020, 0);
// Central frequency [Hz]
constexpr std::array<Numeric, nlines_mpm2020> f0 = {
1.18750334E+11, 5.6264774E+10, 6.2486253E+10, 5.8446588E+10, 6.0306056E+10,
5.9590983E+10, 5.9164204E+10, 6.0434778E+10, 5.8323877E+10, 6.1150562E+10,
5.7612486E+10, 6.1800158E+10, 5.6968211E+10, 6.2411220E+10, 5.6363399E+10,
6.2997984E+10, 5.5783815E+10, 6.3568526E+10, 5.5221384E+10, 6.4127775E+10,
5.4671180E+10, 6.4678910E+10, 5.4130025E+10, 6.5224078E+10, 5.3595775E+10,
6.5764779E+10, 5.3066934E+10, 6.6302096E+10, 5.2542418E+10, 6.6836834E+10,
5.2021429E+10, 6.7369601E+10, 5.1503360E+10, 6.7900868E+10, 5.0987745E+10,
6.8431006E+10, 5.0474214E+10, 6.8960312E+10, 3.68498246E+11, 4.24763020E+11,
4.87249273E+11, 7.15392902E+11, 7.73839490E+11, 8.34145546E+11, };

// Intensity [1 / Pa] (rounded to 10 digits because at most 9 digits exist in f0)
constexpr std::array<Numeric, nlines_mpm2020> intens = {
1.591521878E-21, 1.941172240E-21, 4.834543970E-21, 4.959264029E-21, 7.010386457E-21,
7.051673348E-21, 8.085012578E-21, 8.108262250E-21, 8.145673278E-21, 8.149757320E-21,
7.396406085E-21, 7.401923754E-21, 6.162286575E-21, 6.168475265E-21, 4.749226167E-21,
4.754435107E-21, 3.405982896E-21, 3.408455562E-21, 2.282498656E-21, 2.283934341E-21,
1.432692459E-21, 1.433513473E-21, 8.439995690E-22, 8.443521837E-22, 4.672706507E-22,
4.676049313E-22, 2.435008301E-22, 2.437304596E-22, 1.195038747E-22, 1.196873412E-22,
5.532759045E-23, 5.537261239E-23, 2.416832398E-23, 2.418989865E-23, 9.969285671E-24,
9.977543709E-24, 3.882541154E-24, 3.888101811E-24, 3.676253816E-23, 3.017524005E-22,
9.792882227E-23, 2.756166168E-23, 1.486462215E-22, 4.411918954E-23, };

// Temperature intensity modifier
constexpr std::array<Numeric, nlines_mpm2020> a2 = {
0.01, 0.014, 0.083, 0.083, 0.207, 0.207, 0.387, 0.386,
0.621, 0.621, 0.910, 0.910, 1.255, 1.255, 1.654, 1.654,
2.109, 2.108, 2.618, 2.617, 3.182, 3.181, 3.800, 3.800,
4.474, 4.473, 5.201, 5.200, 5.983, 5.982, 6.819, 6.818,
7.709, 7.708, 8.653, 8.652, 9.651, 9.650, 0.048, 0.044,
0.049, 0.145, 0.141, 0.145};

// Computed arrays values
Eigen::Array<Numeric, necs2020, 1> d;
Eigen::Array<Numeric, necs2020, 1> rho;
Eigen::Matrix<Complex, necs2020, necs2020> W;

// Diagonal and upper triangular values
for (Index i=0; i<necs2020; i++) {
auto Ji = qids[i].UpperQuantumNumber(QuantumNumberType::J);
auto Jf = qids[i].LowerQuantumNumber(QuantumNumberType::J);
auto Ni = qids[i].UpperQuantumNumber(QuantumNumberType::N);
auto Nf = qids[i].LowerQuantumNumber(QuantumNumberType::N);

// Compute reduced dipole
dl[i] = o2_makarov2013_reduced_dipole(Ji, Jf, Ni);

// Do virtual transitions
for (Index j=0; j<necs2020; j++) {
for (Index j=i; j<necs2020; j++) {
auto Ji_p = qids[j].UpperQuantumNumber(QuantumNumberType::J);
auto Jf_p = qids[j].LowerQuantumNumber(QuantumNumberType::J);
auto Ni_p = qids[j].UpperQuantumNumber(QuantumNumberType::N);
auto Nf_p = qids[j].LowerQuantumNumber(QuantumNumberType::N);
if (i not_eq j) {
W(i, j) = o2_ecs_wigner_symbol_tran(Ji, Jf, Ni, Nf, 1, 1, Ji_p, Jf_p, Ni_p, Nf_p, 1, T);
W(i, j) = 1i * o2_ecs_wigner_symbol_tran(Ji, Jf, Ni, Nf, 1, 1, Ji_p, Jf_p, Ni_p, Nf_p, 1, T);
} else {
W(i, j) = (1 + 0.1*water_vmr) * P * lsm[i].compute(T, T0, LineShape::Variable::G0);
W(i, j) = 1i * (1 + 0.1*water_vmr) * P * lsm[i].compute(T, T0, LineShape::Variable::G0);
}
}
}

std::cout << "W = np.array([";
std::transform(qids.cbegin(), qids.cbegin()+necs2020, d.data(), [](auto& qns){return
o2_makarov2013_reduced_dipole (qns.UpperQuantumNumber(QuantumNumberType::J),
qns.LowerQuantumNumber(QuantumNumberType::J),
qns.UpperQuantumNumber(QuantumNumberType::N));});

std::transform(qids.cbegin(), qids.cbegin()+necs2020, rho.data(), [T](auto& qns){return
boltzman_factor(T, o2_ecs_erot_jn_same(qns.LowerQuantumNumber(QuantumNumberType::J)));});

const Numeric theta = T0 / T;
const Numeric theta_m1 = theta - 1;
const Numeric theta_3 = pow3(theta);
for (Index i=0; i<necs2020; i++) {
std::cout << "[";
for (Index j=0; j<necs2020; j++) {
std::cout << W(i, j) << ", ";
}
std::cout << "], ";
const Numeric ST = theta_3 * P * intens[i] * std::exp(-a2[i] * theta_m1);
const Numeric sgn = std::abs(d[i]) / d[i];
d[i] = sgn * f0[i] * std::sqrt(ST / rho[i]);
}
std::cout << "])\n";

normalize_relaxation_matrix(W, rho, d);

/*if (not full)*/
for (Index i=0; i<necs2020; i++) {
W(i, i) += f0[i];
}

Eigen::ComplexEigenSolver<Eigen::Matrix<Complex, necs2020, necs2020>> eV;
eV.compute(W);
auto& D = eV.eigenvalues();
auto& V = eV.eigenvectors();
auto& Vinv = W = eV.eigenvectors().inverse();

Eigen::Array<Complex, necs2020, 1> B; B *= 0;
for (Index m=0; m<necs2020; m++) {
for (Index i=0; i<necs2020; i++) {
for (Index j=0; j<necs2020; j++) {
B[m] += rho[i] * d[i] * d[j] * V(j, m) * Vinv(m, i);
}
}
}

// // Lorentz profile!
// const Numeric div = Constant::inv_pi / rho.cwiseProduct(d.cwiseAbs2()).sum();
// for (Index i=0; i<f.nelem(); i++) {
// I[i] = div * (B.array() / (f[i] - D.array())).sum();

std::cout << "dl = np.array([";
auto GD0 = Linefunctions::DopplerConstant(T, SpeciesTag("O2-66").SpeciesMass());
for (Index i=0; i<necs2020; i++) {
std::cout << dl[i] << ", ";
for (Index iv=0; iv<f.nelem(); iv++) {
I[iv] += (Constant::inv_sqrt_pi / (GD0 * D[i].real())) * Faddeeva::w((f[iv] - std::conj(D[i])) / (GD0 * D[i].real())) * B[i];
}
}
std::cout << "])\n";
}
2 changes: 1 addition & 1 deletion src/predefined_absorption_models.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ void makarov2020_o2_lines_mpm(Matrix& xsec,
const ArrayOfRetrievalQuantity& jacs,
const ArrayOfIndex& jacs_pos);

void makarov2020_o2_lines_ecs(Numeric P, Numeric T, Numeric water_vmr);
void makarov2020_o2_lines_ecs(ComplexVector& I, const Vector& f, Numeric P, Numeric T, Numeric water_vmr);
}; //PredefinedModel
}; //Absorption

Expand Down
Loading

0 comments on commit 08de9d4

Please sign in to comment.