rotate_pol methods, lowercase names

CMB-S4 · Oct 24, 2023 · a5d8685 · a5d8685
1 parent 61e6256
commit a5d8685
Show file tree

Hide file tree

Showing 3 changed files with 54 additions and 14 deletions.
diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h
@@ -271,6 +271,8 @@ class StokesVector {
 
 	StokesVector &operator /=(const MuellerMatrix &r);
 
+	void rotate_pol(double pol_ang);
+
 private:
 	double backing[3];
 };
@@ -356,14 +358,16 @@ class MuellerMatrix {
 		return s;
 	}
 
-	double Det() const {
+	double det() const {
 		return (tt * (qq * uu - qu * qu) -
 			tq * (tq * uu - qu * tu) +
 			tu * (tq * qu - qq * tu));
 	}
 
-	MuellerMatrix Inv() const;
-	double Cond() const;
+	MuellerMatrix inv() const;
+	double cond() const;
+
+	void rotate_pol(double pol_ang);
 
 private:
 	double backing[6];

diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx
@@ -1160,18 +1160,18 @@ pyskymapweights_imultma(G3SkyMapWeightsPtr a, const G3SkyMapMask &b)
 	return a;
 }
 
-MuellerMatrix MuellerMatrix::Inv() const
+MuellerMatrix MuellerMatrix::inv() const
 {
 	MuellerMatrix m;
-	double c = Cond();
+	double c = cond();
 	if (tt == 0 || c != c || c > 1e12) {
 		if ((c != c || c > 1e12) && tt != 0)
 			log_trace("Singular matrix found when inverting!  Cond is %lE\n", c);
 		m.tt = m.tq = m.tu = m.qq = m.qu = m.uu = 0.0 / 0.0;
 		return m;
 	}
 
-	double d = Det();
+	double d = det();
 	m.tt = (qq * uu - qu * qu) / d;
 	m.tq = (tu * qu - tq * uu) / d;
 	m.tu = (tq * qu - tu * qq) / d;
@@ -1182,7 +1182,7 @@ MuellerMatrix MuellerMatrix::Inv() const
 	return m;
 }
 
-double MuellerMatrix::Cond() const
+double MuellerMatrix::cond() const
 {
 	// Compute eigenvalues of a symmetrix 3x3 matrix
 	// See https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3%C3%973_matrices
@@ -1222,7 +1222,7 @@ double MuellerMatrix::Cond() const
 	B.tu = tu;
 	B.qu = qu;
 	B /= p;
-	double r = B.Det() / 2.0;
+	double r = B.det() / 2.0;
 
 	double phi;
 	if (r <= -1)
@@ -1256,9 +1256,45 @@ StokesVector::StokesVector(double pol_ang, double pol_eff) :
 		u = 0.0;
 }
 
+void StokesVector::rotate_pol(double pol_ang)
+{
+	double phi = pol_ang * 2.0 / G3Units::rad;
+	double c = cos(phi);
+	double s = sin(phi);
+
+	double qr = q * c - u * s;
+	double ur = q * s + u * c;
+
+	q = qr;
+	u = ur;
+}
+
+void MuellerMatrix::rotate_pol(double pol_ang)
+{
+	double phi = pol_ang * 2.0 / G3Units::rad;
+	double c = cos(phi);
+	double s = sin(phi);
+	double cc = c * c;
+	double ss = s * s;
+	double cs = c * s;
+	double xx = 2.0 * qu * cs;
+
+	double tqr = tq * c - tu * s;
+	double tur = tq * s + tu * c;
+	double qqr = qq * cc - xx + uu * ss;
+	double qur = (qq - uu) * cs + qu * (cc - ss);
+	double uur = qq * ss + xx + uu * cc;
+
+	tq = tqr;
+	tu = tur;
+	qq = qqr;
+	qu = qur;
+	uu = uur;
+}
+
 StokesVector & StokesVector::operator /=(const MuellerMatrix &r)
 {
-	MuellerMatrix ir = r.Inv();
+	MuellerMatrix ir = r.inv();
 	if (ir.tt != ir.tt)
 		t = q = u = 0.0 / 0.0;
 	else
@@ -1271,7 +1307,7 @@ G3SkyMapPtr G3SkyMapWeights::Det() const
 	G3SkyMapPtr D = TT->Clone(false);
 
 	for (size_t i = 0; i < TT->size(); i++) {
-		double det = this->at(i).Det();
+		double det = this->at(i).det();
 		if (det != 0)
 			(*D)[i] = det;
 	}
@@ -1285,7 +1321,7 @@ G3SkyMapPtr G3SkyMapWeights::Cond() const
 	C->ConvertToDense();
 
 	for (size_t i = 0; i < TT->size(); i++)
-		(*C)[i] = this->at(i).Cond();
+		(*C)[i] = this->at(i).cond();
 
 	return C;
 }
@@ -1306,7 +1342,7 @@ G3SkyMapWeightsPtr G3SkyMapWeights::Inv() const
 		out->UU->ConvertToDense();
 
 	for (size_t i = 0; i < TT->size(); i++)
-		(*out)[i] = this->at(i).Inv();
+		(*out)[i] = this->at(i).inv();
 
 	return out;
 }

diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx
@@ -65,12 +65,12 @@ void RemoveWeights(G3SkyMapPtr T, G3SkyMapPtr Q, G3SkyMapPtr U, G3SkyMapWeightsC
 				if (empty && v == 0)
 					continue;
 			} else {
-				double c = m.Cond();
+				double c = m.cond();
 				empty = (c != c) || (c > 1e12);
 				if (empty && v == 0 && Q->at(pix) == 0 && U->at(pix) == 0)
 					continue;
 				if (!empty)
-					empty |= (m.Det() == 0);
+					empty |= (m.det() == 0);
 			}
 
 			// set bad pixels to 0