Skip to content

Commit

Permalink
merged #21 doxygenization branch into #134, and harmonized the docs f…
Browse files Browse the repository at this point in the history
…or the Shear class. Ellipse is still not well documented though, and I should make some minor changes to the shear.py module to reflect use here.
  • Loading branch information
barnabytprowe committed Jun 19, 2012
2 parents 403eda5 + aa3dafa commit 109a737
Showing 1 changed file with 143 additions and 32 deletions.
175 changes: 143 additions & 32 deletions include/galsim/Shear.h
Original file line number Diff line number Diff line change
@@ -1,12 +1,13 @@
// -*- c++ -*-
#ifndef SHEAR_H
#define SHEAR_H
/**
* @file Shear.h Contains a class definition for Shear and Ellipse.
*
* Shear is used to represent shape distortions; Ellipse includes shear, translation and
* magnification.
*/

#ifndef SHEAR_H
#define SHEAR_H
#include <cmath>
#include "TMV.h"

Expand All @@ -30,17 +31,22 @@ namespace galsim {
* intrinsic shapes, and due to lensing shears. Given semi-major and semi-minor axis lengths
* "a" and "b", there are numerous ways to represent shears:
*
* eta = "conformal shear", a/b = exp(eta)
* g = "reduced shear", g=(a-b)/(a+b)
* e = "distortion", e=(a^2-b^2)/(a^2+b^2)
* q = "axis ratio", q=b/a
* To specify both components, we have a value Beta that is the
* real-space position angle of the major axis, i.e. e1 = e cos(2*Beta) and e2 = e sin(2*Beta).
* eta = "conformal shear", a/b = exp(|eta|)
*
* g = "reduced shear", |g| = (a-b)/(a+b)
*
* e = "distortion", |e| = (a^2-b^2)/(a^2+b^2)
*
* q = "axis ratio", |q|=b/a
*
* To specify both components, we have a value Beta that is the real-space position angle of
* the major axis, i.e. e1 = e cos(2*Beta) and e2 = e sin(2*Beta).
*
* Shears are represented internally by e1 and e2, which relate to second moments via
* e1 = (Mxx - Myy) / (Mxx + Myy)
* e2 = 2 Mxy / (Mxx + Myy)
* however, given that lensing specialists most commonly think in terms of reduced shear, the
* e1 = (Mxx - Myy) / (Mxx + Myy),
* e2 = 2 Mxy / (Mxx + Myy).
*
* However, given that lensing specialists most commonly think in terms of reduced shear, the
* constructor that takes two numbers expects (g1, g2). A user who wishes to specify another
* type of shape should use methods, i.e.
* s = Shear();
Expand All @@ -52,16 +58,23 @@ namespace galsim {
friend Shear operator*(const double, const Shear& );

public:
// Construct w/o variance
/**
* @brief Construct without variance / without initializing transformation matrix.
*
* @param[in] g1 first component (reduced shear definition).
* @param[in] g2 second shear component (reduced shear definition).
*/
explicit Shear(double g1=0., double g2=0.) :
hasMatrix(false), matrixA(0), matrixB(0), matrixC(0)
{ setG1G2(g1, g2); }

/// @brief Copy constructor.
Shear(const Shear& rhs) :
e1(rhs.e1), e2(rhs.e2), hasMatrix(rhs.hasMatrix),
matrixA(rhs.matrixA), matrixB(rhs.matrixB), matrixC(rhs.matrixC)
{}

/// @brief Copy assignment.
const Shear& operator=(const Shear& s)
{
e1 = s.e1; e2=s.e2;
Expand All @@ -70,74 +83,162 @@ namespace galsim {
return *this;
}

/// @brief Set (e1, e2) using distortion definition.
Shear& setE1E2(double e1in=0., double e2in=0.);

/**
* @brief Set (|e|, beta) polar ellipticity representation using distortion definition.
* beta must be an Angle.
*/
Shear& setEBeta(double etain=0., Angle betain=Angle());

/// @brief Set (eta1, eta2) using conformal shear definition.
Shear& setEta1Eta2(double eta1in=0., double eta2in=0.);
Shear& setEtaBeta(double etain=0., Angle betain=Angle());

/// @brief Set (|eta|, beta) using conformal shear definition. beta must be an Angle.
Shear& setEtaBeta(double =0., Angle betain=Angle());

/// @brief set (g1, g2) using reduced shear |g| = (a-b)/(a+b) definition.
Shear& setG1G2(double g1in=0., double g2in=0.);

/// @brief Get e1 using distortion definition.
double getE1() const { return e1; }

/// @brief Get e2 using distortion definition.
double getE2() const { return e2; }

/// @brief Get |e| using distortion definition.
double getE() const { return std::sqrt(e1*e1+e2*e2); }

/// @brief Get |e|^2 using distortion definition.
double getESq() const { return e1*e1+e2*e2; }

/// @brief Get polar angle beta (returns an Angle class object).
Angle getBeta() const { return std::atan2(e2,e1)*0.5 * radians; }

/// @brief Get |eta| using conformal shear definition.
double getEta() const { return atanh(std::sqrt(e1*e1+e2*e2)); } //error checking?

// g = gamma / (1-kappa)
/// @brief Get |g| using reduced shear |g| = (a-b)/(a+b) definition.
double getG() const
{
double e=getE();
return e>0. ? (1-std::sqrt(1-e*e))/e : 0.;
}

/// @brief Get g1 using reduced shear |g| = (a-b)/(a+b) definition.
double getG1() const
{
double esq = getESq();
double scale = (esq>1.e-6) ? (1.-std::sqrt(1.-esq))/esq : 0.5;
return e1*scale;
}

/// @brief Get g1 using reduced shear |g| = (a-b)/(a+b) definition.
double getG2() const
{
double esq = getESq();
double scale = (esq>1.e-6) ? (1.-std::sqrt(1.-esq))/esq : 0.5;
return e2*scale;
}

/**
* @brief Get (eta1, eta2) using conformal shear definition.
*
* @param[in,out] eta1 Reference to eta1 variable.
* @param[in,out] eta2 Reference to eta2 variable.
*
*/
void getEta1Eta2(double& eta1, double& eta2) const;

/**
* @brief Get (g1, g2) using reduced shear |g| = (a-b)/(a+b) definition.
*
* @param[in,out] g1 Reference to g1 variable.
* @param[in,out] g2 Reference to g2 variable.
*
*/
void getG1G2(double& g1, double& g2) const;


//negation
/// @brief Unary negation (both components negated).
Shear operator-() const
{
double esq = getESq();
double scale = (esq>1.e-6) ? (1.-std::sqrt(1.-esq))/esq : 0.5;
return esq>0. ? Shear(-e1*scale, -e2*scale) : Shear(0.0, 0.0);
}

// Composition operation: returns ellipticity of
// circle that is sheared first by RHS and then by
// LHS Shear. Note that this addition is
// ***not commutative***!
// In the += and -= operations, this is LHS
// and the operand is RHS of + or - .
/**
* @brief Composition operation.
*
* Note that this 'addition' is ***not commutative***!
*
* @returns Ellipticity of circle that is sheared first by RHS and then by
* LHS Shear.
*/
Shear operator+(const Shear& ) const;

/**
* @brief Composition (with RHS negation) operation.
*
* Note that this 'subtraction' is ***not commutative***!

This comment has been minimized.

Copy link
@rmjarvis

rmjarvis Jun 19, 2012

Member

Is subtraction ever commutative? That would be an interesting algebra...

This comment has been minimized.

Copy link
@rmandelb

rmandelb via email Jun 19, 2012

Member

This comment has been minimized.

Copy link
@gbernstein

gbernstein Jun 19, 2012

Member

I perpetually need to go back to my notes to decide which order I need to apply shear addition. Neither order is uniquely obvious to me.

This comment has been minimized.

Copy link
@rmandelb

rmandelb Jun 19, 2012

Member

That's the problem I'm having. I go back and forth between thinking A + B should mean that A came along first so it gets applied first, and thinking that if you have (A+B) foo where foo is some galaxy, then B hits foo and gets applied first. Hmmm. I am trying to think through implications for users - e.g., if they make a Shear representing the intrinsic shape, and another Shear representing the lensing shear, then you could imagine them doing foo.applyShear(intrinsic) followed by foo.applyShear(lensing) which would be fine - or you could see someone doing foo.applyShear(intrinsic + lensing) which would be wrong. Obviously config will take care of this for many users, but I could see developers getting tripped up by this. Actually as a developer I might be inclined to say that neither order is obviously right so I will side-step the problem by always using applyShear twice.

This comment has been minimized.

Copy link
@rmjarvis

rmjarvis Jun 19, 2012

Member

I think neither version is consistently intuitive. That's just the nature of this particular beast.

As long as it's well documented, I think we are fine. Users will always need to check the doc and think through which order they mean. And that will be true no matter which way we define it. I agree that the typical usage will be to eschew the + operator and just do two applyShear calls, as you suggested, which will be easier to think through.

In any case, I recommend just leaving it as is.

This comment has been minimized.

Copy link
@rmandelb

rmandelb Jun 19, 2012

Member

I think I should at least add something about this to the python dox.

This comment has been minimized.

Copy link
@barnabytprowe

barnabytprowe via email Jun 19, 2012

Author Member

This comment has been minimized.

Copy link
@barnabytprowe

barnabytprowe Jun 19, 2012

Author Member

Ah, I see you've already been harmonizing, and indeed on a commit I commented on by email. Sorry! :)

*
* @returns Ellipticity of circle that is sheared first by the negative RHS and then by
* LHS Shear.
*/
Shear operator-(const Shear& ) const;

// In the += and -= operations, this is LHS
// and the operand is RHS of + or - .

/**
* @brief Inplace composition operation.
*
* Note that this 'addition' is ***not commutative***!
*
* In the += operation, this is LHS and the operand is RHS of +.
*
* @returns Ellipticity of circle that is sheared first by RHS and then by
* LHS Shear.
*/
Shear& operator+=(const Shear& );

/**
* @brief Inplace composition (with RHS negation) operation.
*
* Note that this 'addition' is ***not commutative***!
*
* In the -= operation, this is LHS and the operand is RHS of -.
*
* @returns Ellipticity of circle that is sheared first by the negative RHS and then by
* LHS Shear.
*/
Shear& operator-=(const Shear& );

// Give the rotation angle for this+rhs:
Angle rotationWith(const Shear& rhs) const;
// Detail on the above: s1 + s2 operation on points in
// the plane induces a rotation as well as a shear.
// Above method tells you what the rotation was for LHS+RHS
/** @brief Give the rotation angle for this+rhs.
*
* Detail on the above: s1 + s2 operation on points in
* the plane induces a rotation as well as a shear.
* Above method tells you what the rotation was for LHS+RHS.
*/
Angle rotationWith(const Shear& rhs) const;

///@brief Test equivalence by comparing e1 and e2.
bool operator==(const Shear& rhs) const
{ return e1==rhs.e1 && e2==rhs.e2; }

///@brief Test non-equivalence by comparing e1 and e2.
bool operator!=(const Shear& rhs) const
{ return e1!=rhs.e1 || e2!=rhs.e2; }

// Classes that treat shear as a point-set map:
// Classes thattreat shear as a point-set map:
/**
* @brief Forward transformation from image to source plane coordinates under shear.
*
* @param[in] p 2D vector Position in image plane.
*
* @returns 2D vector Position in source plane.
*/
template <class T>
Position<T> fwd(Position<T> p) const
{
Expand All @@ -147,6 +248,13 @@ namespace galsim {
return out;
}

/**
* @brief Inverse transformation from source to image plane coordinates under shear.
*
* @param[in] p 2D vector Position in source plane.
*
* @returns 2D vector Position in image plane.
*/
template <class T>
Position<T> inv(Position<T> p) const
{
Expand All @@ -156,10 +264,13 @@ namespace galsim {
return out;
}

// Get matrix representation of the forward transformation.
// Matrix is ... and in limit of small shear:
// ( a c = ( 1+g1 g2
// c b ) g2 1-g1 )
/**
* @brief Get matrix representation of the forward transformation.
*
* Matrix is ... and in limit of small shear:
* ( a c = ( 1+g1 g2
* c b ) g2 1-g1 )
*/
void getMatrix(double& a, double& b, double& c) const
{ calcMatrix(); a=matrixA; b=matrixB; c=matrixC; }

Expand Down

0 comments on commit 109a737

Please sign in to comment.