Skip to content

Commit

Permalink
StokesVector and MuellerMatrix Improvements
Browse files Browse the repository at this point in the history
* Add StokesVector constructor from polarization angle and efficiency
* Add MuellerMatrix construct from StokesVector components
* Add MuellerMatrix.from_vector method for populating an existing MuellerMatrix
  from a StokesVector
* Use these features in MapBinner and MapMockObserver
  • Loading branch information
arahlin committed Oct 24, 2023
1 parent b02c746 commit 61e6256
Show file tree
Hide file tree
Showing 4 changed files with 40 additions and 54 deletions.
14 changes: 14 additions & 0 deletions maps/include/maps/G3SkyMap.h
Original file line number Diff line number Diff line change
Expand Up @@ -241,6 +241,9 @@ class StokesVector {
// you initialize one of these statically and then use it in arithmetic.
StokesVector() : t(backing[0]), q(backing[1]), u(backing[2]) {}

// Constructor for polarization coupling
StokesVector(double pol_ang, double pol_eff);

double &t, &q, &u;

StokesVector &operator +=(const StokesVector &r) {
Expand Down Expand Up @@ -292,8 +295,19 @@ class MuellerMatrix {
MuellerMatrix() : tt(backing[0]), tq(backing[1]),
tu(backing[2]), qq(backing[3]), qu(backing[4]), uu(backing[5]) {}

// Constructor from Stokes vector
MuellerMatrix(const StokesVector &r) : tt(backing[0]), tq(backing[1]),
tu(backing[2]), qq(backing[3]), qu(backing[4]), uu(backing[5]) {
from_vector(r);
}

double &tt, &tq, &tu, &qq, &qu, &uu;

void from_vector(const StokesVector &r) {
tt = r.t * r.t; tq = r.t * r.q; tu = r.t * r.u;
qq = r.q * r.q; qu = r.q * r.u; uu = r.u * r.u;
}

MuellerMatrix &operator +=(const MuellerMatrix &r) {
tt += r.tt; tq += r.tq; tu += r.tu;
qq += r.qq; qu += r.qu; uu += r.uu;
Expand Down
16 changes: 16 additions & 0 deletions maps/src/G3SkyMap.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -1240,6 +1240,22 @@ double MuellerMatrix::Cond() const
return c;
}

StokesVector::StokesVector(double pol_ang, double pol_eff) :
t(backing[0]), q(backing[1]), u(backing[2])
{
double phi = pol_ang / G3Units::rad * 2.0;
double gamma = pol_eff / (2.0 - pol_eff);

t = 1.0;
q = cos(phi) * gamma;
u = sin(phi) * gamma;

if (fabs(q) < 1e-12)
q = 0.0;
if (fabs(u) < 1e-12)
u = 0.0;
}

StokesVector & StokesVector::operator /=(const MuellerMatrix &r)
{
MuellerMatrix ir = r.Inv();
Expand Down
34 changes: 2 additions & 32 deletions maps/src/MapBinner.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -335,35 +335,6 @@ MapBinner::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
out.push_back(frame);
}

// Following two utility functions copied from the old map-maker and need to
// be re-tooled -- either moved into the constructors of the relevant objects
// and/or improved to handle things like boresight rotation.
static void
fill_mueller_matrix_from_stokes_coupling(
const StokesVector &stokes_coupling, MuellerMatrix &pcm)
{
pcm.tt = stokes_coupling.t * stokes_coupling.t;
pcm.tq = stokes_coupling.t * stokes_coupling.q;
pcm.tu = stokes_coupling.t * stokes_coupling.u;
pcm.qq = stokes_coupling.q * stokes_coupling.q;
pcm.qu = stokes_coupling.q * stokes_coupling.u;
pcm.uu = stokes_coupling.u * stokes_coupling.u;
}

static void
set_stokes_coupling(double pol_ang, double pol_eff,
StokesVector &stokes_coupling)
{
stokes_coupling.t = 1.0;
stokes_coupling.q = cos(pol_ang/G3Units::rad *2.)*pol_eff/(2.0-pol_eff);
stokes_coupling.u = sin(pol_ang/G3Units::rad *2.)*pol_eff/(2.0-pol_eff);

if (fabs(stokes_coupling.q) < 1e-12)
stokes_coupling.q = 0;
if (fabs(stokes_coupling.u) < 1e-12)
stokes_coupling.u = 0;
}

void
MapBinner::BinTimestream(const G3Timestream &det, double weight,
const BolometerProperties &bp, const G3VectorQuat &pointing,
Expand All @@ -382,12 +353,11 @@ MapBinner::BinTimestream(const G3Timestream &det, double weight,
// And polarization coupling
// XXX: does not handle focal-plane rotation, since it assumes
// polarization coupling is constant for the whole scan.
StokesVector pcoupling;
set_stokes_coupling(bp.pol_angle, bp.pol_efficiency, pcoupling);
StokesVector pcoupling(bp.pol_angle, bp.pol_efficiency);

MuellerMatrix mueller;
if (W) {
fill_mueller_matrix_from_stokes_coupling(pcoupling, mueller);
mueller.from_vector(pcoupling);
mueller *= weight;
}

Expand Down
30 changes: 8 additions & 22 deletions maps/src/MapMockObserver.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ class MapMockObserver : public G3Module {

bool interp_;
bool error_on_zero_;
int pol_sign_;

SET_LOGGER("MapMockObserver");
};
Expand Down Expand Up @@ -109,26 +110,12 @@ MapMockObserver::MapMockObserver(std::string pointing, std::string timestreams,
{
if ((Q_ && !U_) || (U_ && !Q_))
log_fatal("If simulating polarized maps, pass both Q and U.");
}

// Following utility function copied from old map-maker and needs re-tooling
// to handle boresight rotation.
static void
set_stokes_coupling(double pol_ang, double pol_eff,
StokesVector &stokes_coupling, G3SkyMap::MapPolConv pol_conv)
{
stokes_coupling.t = 1.0;
stokes_coupling.q = cos(pol_ang/G3Units::rad *2.)*pol_eff/(2.0-pol_eff);
stokes_coupling.u = sin(pol_ang/G3Units::rad *2.)*pol_eff/(2.0-pol_eff);
if (pol_conv == G3SkyMap::COSMO)
stokes_coupling.u *= -1.0;
else if (pol_conv == G3SkyMap::ConvNone)
log_fatal("Missing pol_conv");

if (fabs(stokes_coupling.q) < 1e-12)
stokes_coupling.q = 0;
if (fabs(stokes_coupling.u) < 1e-12)
stokes_coupling.u = 0;
if (U_) {
if (U_->GetPolConv() == G3SkyMap::ConvNone)
log_fatal("Missing pol_conv");
pol_sign_ = (U_->GetPolConv() == G3SkyMap::COSMO) ? -1 : 1;
}
}

void
Expand Down Expand Up @@ -210,9 +197,8 @@ MapMockObserver::Process(G3FramePtr frame, std::deque<G3FramePtr> &out)
detpointing = T_->AnglesToPixels(alpha, delta);

if (Q_) {
StokesVector pcoupling;
set_stokes_coupling(bp.pol_angle, bp.pol_efficiency,
pcoupling, U_->GetPolConv());
// needs re-tooling to handle boresight rotation.
StokesVector pcoupling(bp.pol_angle * pol_sign_, bp.pol_efficiency);
for (size_t i = 0; i < det.size(); i++) {
if (interp_) {
std::vector<size_t> pixels;
Expand Down

0 comments on commit 61e6256

Please sign in to comment.