diff --git a/cosipy/polarization/__init__.py b/cosipy/polarization/__init__.py
new file mode 100644
index 00000000..eadc83f2
--- /dev/null
+++ b/cosipy/polarization/__init__.py
@@ -0,0 +1,2 @@
+from .conventions import PolarizationConvention, OrthographicConvention, StereographicConvention, IAUPolarizationConvention
+from .polarization_angle import PolarizationAngle
diff --git a/cosipy/polarization/conventions.py b/cosipy/polarization/conventions.py
new file mode 100644
index 00000000..a44eea01
--- /dev/null
+++ b/cosipy/polarization/conventions.py
@@ -0,0 +1,296 @@
+import numpy as np
+from astropy.coordinates import SkyCoord, Angle
+import astropy.units as u
+import inspect
+from scoords import Attitude, SpacecraftFrame
+
+# Base class for polarization conventions
+class PolarizationConvention:
+
+    def __init__(self):
+        """
+        Base class for polarization conventions
+        """
+
+    _registered_conventions = {}
+        
+    @classmethod
+    def register(cls, name):
+
+        name = name.lower()
+        
+        def _register(convention_class):
+            cls._registered_conventions[name] = convention_class
+            return convention_class 
+        return _register
+
+    @classmethod
+    def get_convention(cls, name, *args, **kwargs):
+
+        if inspect.isclass(name):
+            if issubclass(name, PolarizationConvention):
+                return name(*args, **kwargs)
+            else:
+                raise TypeError("Class must be subclass of PolarizationConvention")
+
+        if isinstance(name, PolarizationConvention):
+            return name
+
+        if not isinstance(name, str):
+            raise TypeError("Input must be str, or PolarizationConvention subclass or object")
+        
+        name = name.lower()
+
+        try:
+            return cls._registered_conventions[name](*args, **kwargs)
+        except KeyError as e:
+            raise Exception(f"No polarization convention by name '{name}'") from e
+
+    @property
+    def frame(self):
+        """
+        Astropy coordinate frame
+        """
+        return None
+        
+    def get_basis(self, source_direction: SkyCoord):
+        """
+        Get the px,py unit vectors that define the polarization plane on 
+        this convention. Polarization angle increments from px to py.
+
+        Parameters
+        ----------
+        source_direction : SkyCoord
+            The direction of the source
+
+        Return
+        ------
+        px,py : SkyCoord
+            Polarization angle increaes from px to py. pz is always
+            the opposite of the source direction --i.e. in the direction of the
+            particle.
+        """
+
+        
+# Orthographic projection convention
+class OrthographicConvention(PolarizationConvention):
+
+    def __init__(self,
+                 ref_vector: SkyCoord = None,
+                 clockwise: bool = False):
+        """
+        The local polarization x-axis points towards an arbitrary reference vector, 
+        and the polarization angle increasing counter-clockwise when looking 
+        at the source.
+        
+        Parameters
+        ----------
+        ref_vector : SkyCoord
+            Set the reference vector, defaulting to celestial north if not provided 
+            (IAU convention)
+        clockwise : bool
+            Direction of increasing PA, when looking at the source. Default is false 
+            --i.e. counter-clockwise when looking outwards.
+            
+        """
+        if ref_vector is None:
+            self.ref_vector = SkyCoord(ra=0 * u.deg, dec=90 * u.deg, frame="icrs")
+        else:
+            self.ref_vector = ref_vector
+
+        self._sign = 1 if clockwise else -1
+
+    def __repr__(self):
+        return f"<OrthographicConvention(starting from {self.ref_vector} and {'clockwise' if self.is_clockwise else 'counter-clockwise'} when looking at the source)>"
+
+    @property
+    def is_clockwise(self):
+        """
+        When looking at the source
+        """
+        return True if self._sign == 1 else False
+    
+    @property
+    def frame(self):
+        return self.ref_vector.frame
+        
+    def get_basis(self, source_direction: SkyCoord):
+        # Extract Cartesian coordinates for the source direction.
+        pz = self._sign * source_direction.transform_to(self.frame).cartesian.xyz
+
+        # Broadcast reference vector
+        ref = np.expand_dims(self.ref_vector.cartesian.xyz,
+                             axis = tuple(np.arange(1,pz.ndim, dtype = int)))
+
+        # Get py. Normalize because pz and ref dot not make 90deg angle
+        py = np.cross(pz, ref, axisa = 0, axisb = 0, axisc = 0)
+        py /= np.linalg.norm(py, axis = 0, keepdims = True)
+
+        # Get px
+        px = np.cross(py, pz, axisa = 0, axisb = 0, axisc = 0)
+        
+        # To SkyCoord
+        px = SkyCoord(*px, representation_type='cartesian', frame = self.frame)
+        py = SkyCoord(*py, representation_type='cartesian', frame = self.frame)
+        
+        return px, py
+
+
+#https://github.com/zoglauer/megalib/blob/1eaad14c51ec52ad1cb2399a7357fe2ca1074f79/src/cosima/src/MCSource.cc#L3452
+class MEGAlibRelative(OrthographicConvention):
+
+    def __init__(self, axis, attitude = None):
+        """
+        Use a polarization vector which is created the following way: 
+        Create an initial polarization vector which is orthogonal on the 
+        initial flight direction vector of the particle and the given axis vector 
+        (e.g. x-axis for RelativeX). This is a simple crossproduct. Then rotate 
+        the polarization vector (right-hand-way) around the initial flight 
+        direction vector of the particle by the given rotation angle. 
+        """
+
+        if not isinstance(axis, str):
+            raise TypeError("Axis must be a string. 'x', 'y' or 'z'.")
+
+        axis = axis.lower()
+        
+        if axis == 'x':
+            ref_vector = SkyCoord(lon=0 * u.deg, lat=0 * u.deg,
+                                  frame = SpacecraftFrame(attitude = attitude))
+        elif axis == 'y':
+            ref_vector = SkyCoord(lon=90 * u.deg, lat=0 * u.deg,
+                                  frame = SpacecraftFrame(attitude = attitude))
+        elif axis == 'z':
+            ref_vector = SkyCoord(lon=0 * u.deg, lat=90 * u.deg,
+                                  frame = SpacecraftFrame(attitude = attitude))
+        else:
+            raise ValueError("Axis must be 'x', 'y' or 'z'.")
+            
+        super().__init__(ref_vector, clockwise = False)
+        
+    def get_basis(self, source_direction: SkyCoord):
+
+        # The MEGAlib and orthographic definitions are prett much the same, but
+        # they differ on the order of the cross products
+        
+        # In MEGAlib definition
+        # pz = -source_direction = particle_direction
+        # px = particle_direction x ref_vector  = pz x ref_vector
+        # py = pz x px
+
+        # In the base orthographic definition
+        # pz = -source_direction = particle_direction
+        # px = py x pz
+        # py = source_direction x ref_vector = -pz x ref_vector
+
+        # Therefore
+        # px = py_base
+        # py = -px_base
+
+        # MEGAlib's PA is counter-clockwise when looking at the sourse
+
+        # Flip px <-> py
+        py,px = super().get_basis(source_direction)
+
+        # Sign of px
+        py = SkyCoord(-py.cartesian,
+                      representation_type = 'cartesian',
+                      frame = py.frame)
+        
+        return px,py
+
+@PolarizationConvention.register("RelativeX")    
+class MEGAlibRelativeX(MEGAlibRelative):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__('x', *args, **kwargs)
+        
+@PolarizationConvention.register("RelativeY")    
+class MEGAlibRelativeY(MEGAlibRelative):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__('y', *args, **kwargs)
+
+@PolarizationConvention.register("RelativeZ")    
+class MEGAlibRelativeZ(MEGAlibRelative):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__('z', *args, **kwargs)
+
+    
+# https://lambda.gsfc.nasa.gov/product/about/pol_convention.html
+# https://www.iau.org/static/resolutions/IAU1973_French.pdf
+@PolarizationConvention.register("IAU")
+class IAUPolarizationConvention(OrthographicConvention):
+
+    def __init__(self):
+        """
+        The following resolution was adopted by Commissions 25 and 40: 
+        'RESOLVED, that the frame of reference for the Stokes parameters 
+        is that of Right Ascension and Declination with the position 
+        angle of electric-vector maximum, e, starting from North and 
+        increasing through East.
+        """
+        super().__init__(ref_vector = SkyCoord(ra=0 * u.deg, dec=90 * u.deg,
+                                               frame="icrs"),
+                         clockwise = False)
+    
+    
+# Stereographic projection convention
+class StereographicConvention(PolarizationConvention):
+
+    def __init__(self,
+                 clockwise: bool = False,
+                 attitude: Attitude = None):
+        """
+        Basis vector follow the steregraphic projection lines. Meant to describe
+        polarization in spacecraft coordinate by minimizing the number of undefined location withing the field of view.
+        
+        Near the boresight --i.e. on axis, center of the FoV, north pole-- it is
+        similar to 
+        ``OrthographicConvention(ref_vector = SkyCoord(lon = 0*u.deg, lat = 0*u.deg, frame = SpacecraftFrame())``
+        however, it has a single undefined point on the opposite end --i.e. south pole,
+        back of the detector---
+        
+        
+        Parameters
+        ----------
+        clockwise : bool
+            Direction of increasing PA, when looking at the source. Default is false 
+            --i.e. counter-clockwise when looking outwards.
+        attitude : Attitude
+            Spacecraft orientation
+        """
+
+        self._attitude = attitude
+        
+        self._sign = 1 if clockwise else -1
+
+    @property 
+    def frame(self):
+        return SpacecraftFrame(attitude = self._attitude)
+        
+    def get_basis(self, source_direction: SkyCoord):
+        # Extract Cartesian coordinates for the source direction
+        x, y, z = source_direction.cartesian.xyz
+
+        # Calculate the projection of the reference vector in stereographic coordinates
+        px_x = 1 - (x**2 - y**2) / (z + 1) ** 2
+        px_y = -2 * x * y / (z + 1) ** 2
+        px_z = -2 * x / (z + 1)
+
+        # Combine the components into the projection vector px
+        px = np.array([px_x, px_y, px_z])
+
+        # Normalize the projection vector
+        norm = np.linalg.norm(px, axis=0)
+        px /= norm
+
+        # Calculate the perpendicular vector py using the cross product
+        py = self._sign*np.cross([x, y, z], px, axis=0)
+
+        # To SkyCoord
+        px = SkyCoord(*px, representation_type='cartesian', frame = self.frame)
+        py = SkyCoord(*py, representation_type='cartesian', frame = self.frame)
+        
+        return px, py
diff --git a/cosipy/polarization/polarization_angle.py b/cosipy/polarization/polarization_angle.py
new file mode 100644
index 00000000..07aa04ac
--- /dev/null
+++ b/cosipy/polarization/polarization_angle.py
@@ -0,0 +1,103 @@
+import numpy as np
+from astropy.coordinates import SkyCoord, Angle
+import astropy.units as u
+
+from .conventions import PolarizationConvention
+
+class PolarizationAngle:
+
+    def __init__(self, angle, skycoord ,
+                 convention = 'iau',
+                 *args, **kwargs):
+        """
+        Defines a polarization angle in the context of a source direction and
+        polarization angle convention.
+
+        Parameters:
+        angle : :py:class:`astropy.coordinates.Angle
+            Polarization angle
+        skycoord : :py:class:`astropy.coordinates.SkyCoord`
+            Source direction
+        convention : PolarizationConvention
+            Convention the defined the polarization basis and direction in 
+            the polarization plane (for which the source direction is normal)
+        *args, **kwargs
+            Passed to convention class.
+        """
+
+        # Ensure pa is an Angle object
+        self._angle = Angle(angle)
+
+        self._convention = PolarizationConvention.get_convention(convention,
+                                                                 *args, **kwargs)
+
+        self._skycoord = skycoord
+
+    def __repr__(self):
+        return f"<PolarizationAngle({self._angle.degree} deg at {self._skycoord} using convention {self._convention})>"
+        
+    @property
+    def angle(self):
+        return self._angle
+
+    @property
+    def convention(self):
+        return self._convention
+
+    @property
+    def skycoord(self):
+        return self._skycoord
+
+    @property
+    def vector(self):
+        """
+        Direction of the electric field vector
+        """
+
+        # Get the projection vectors for the source direction in the current convention
+        px, py = self.convention.get_basis(self.skycoord)
+
+        px = px.cartesian.xyz
+        py = py.cartesian.xyz
+
+        # Calculate the cosine and sine of the polarization angle
+        cos_pa = np.cos(self.angle.radian)
+        sin_pa = np.sin(self.angle.radian)
+
+        # Calculate the polarization vector
+        pol_vec = px * cos_pa + py * sin_pa
+                
+        return SkyCoord(*pol_vec,
+                        representation_type = 'cartesian',
+                        frame = self.convention.frame)
+    
+    def transform_to(self, convention, *args, **kwargs):
+
+        # Standarize convention 2
+        convention2 = PolarizationConvention.get_convention(convention, *args, **kwargs)
+        
+        # Calculate the polarization vector in the current convention
+        pol_vec = self.vector.transform_to(convention2.frame).cartesian.xyz
+
+        # Get the projection vectors for the source direction in the new convention
+        (px2, py2) = convention2.get_basis(self.skycoord)
+
+        px2 = px2.cartesian.xyz
+        py2 = py2.cartesian.xyz
+        
+        # Compute the dot products for the transformation
+        a = np.sum(pol_vec * px2, axis=0)
+        b = np.sum(pol_vec * py2, axis=0)
+
+        # Calculate the new polarization angle in the new convention
+        pa_2 = Angle(np.arctan2(b, a), unit=u.rad)
+
+        # Normalize the angle to be between 0 and pi
+        if pa_2 < 0:
+            pa_2 += Angle(np.pi, unit=u.rad)
+
+        return PolarizationAngle(pa_2,
+                                 self.skycoord,
+                                 convention = convention2)
+
+        
diff --git a/docs/api/figures/OrthographicZ-IAU.png b/docs/api/figures/OrthographicZ-IAU.png
new file mode 100644
index 00000000..b2487860
Binary files /dev/null and b/docs/api/figures/OrthographicZ-IAU.png differ
diff --git a/docs/api/figures/RelativeX.png b/docs/api/figures/RelativeX.png
new file mode 100644
index 00000000..4df5ed83
Binary files /dev/null and b/docs/api/figures/RelativeX.png differ
diff --git a/docs/api/figures/RelativeY.png b/docs/api/figures/RelativeY.png
new file mode 100644
index 00000000..15626208
Binary files /dev/null and b/docs/api/figures/RelativeY.png differ
diff --git a/docs/api/figures/RelativeZ.png b/docs/api/figures/RelativeZ.png
new file mode 100644
index 00000000..980b2380
Binary files /dev/null and b/docs/api/figures/RelativeZ.png differ
diff --git a/docs/api/figures/Stereographic.png b/docs/api/figures/Stereographic.png
new file mode 100644
index 00000000..48beb5ff
Binary files /dev/null and b/docs/api/figures/Stereographic.png differ
diff --git a/docs/api/index.rst b/docs/api/index.rst
index 8754fba6..9b28ddb0 100644
--- a/docs/api/index.rst
+++ b/docs/api/index.rst
@@ -18,7 +18,7 @@ If you are instead interested in an overview on how to use cosipy, see out `tuto
    threeml
    ts_map
    image_deconvolution
-   util
+   polarization
    source_injector          
-   
+   util
 
diff --git a/docs/api/polarization.rst b/docs/api/polarization.rst
new file mode 100644
index 00000000..396142f4
--- /dev/null
+++ b/docs/api/polarization.rst
@@ -0,0 +1,74 @@
+Polarization
+============
+
+Conventions
+-----------
+
+When specifying a polarization angle (PA) it is important to specify the direction of the source and the convention used. All of this is necessary to reconstruct the polarization vector --i.e. the direction of the electric field.
+
+A polarization convention is defined by a set of (px,py) unit vectors (a basis) that define a plane at each location on the sphere s = [x,y,z]. The location s is always normal to the polarization plane defined by (px,py). The polarization angle increases from px to py.
+
+Cosipy supports two major convention constructions: orthographic and stereographic. These are named after the corresponding conformal projection transformation.
+
+In addition for the overall construction prescription, a polarization convention fully defined by the reference frame and whether the polarization angle increases clockwise or counter-clockwise when looking at the source.
+
+Orthographic
+~~~~~~~~~~~~
+
+This the regular polarization convention prescription found in the literature, where the px,py are aligned with the meridians and parallels, and one of them always points towards an arbitrary reference vector. For eample, the following figure corresponds to orthographic convention in ICRS coordinates, pointing towards the North-pole and counter-clockwise when looking at the source (only one quadrant is plotted for easier visualization)::
+
+  OrthographicConvention(ref_vector = SkyCoord(ra = 0*u.deg, dec = 90*u.deg, frame = 'icrs'), clockwise = False)
+
+.. image:: figures/OrthographicZ-IAU.png
+  :width: 400
+  :alt: OrthographicZ convention / IAU
+
+This matches the `IAU convention <https://www.iau.org/static/resolutions/IAU1973_French.pdf>`_ and can be also called by name::
+
+  PolarizationConvention.get_convention("IAU")
+
+Other special conventions of the orthographic family that can also be called by name are MEGAlibs RelativeX/Y/Z, where the -py vector points towards the reference vector (as opposed to px):
+
+.. image:: figures/RelativeZ.png
+  :width: 400
+  :alt: MEGAlib RelativeZ
+
+.. image:: figures/RelativeX.png
+  :width: 400
+  :alt: MEGAlib RelativeX
+
+.. image:: figures/RelativeY.png
+  :width: 400
+  :alt: MEGAlib RelativeY
+
+
+Stereographic
+~~~~~~~~~~~~~
+
+An issue of the orthographic projection is that it is not well-defined at the location of the reference vector and its antipodal direction. This can cause issues due to numerical errors when construction a detector response near those locations. The stereographic convention is meant to solve this issue by describing the polarization angle in spacecraft coordinates with only a single undefined location at -z, where the effective area is almost null. Due to the hairy ball theorem it is impossible to obtain a convention where the polarization basis vector have a smooth transition and with undefined locations.
+
+Near the z-axis, it is very similar to using an orthographic convention with +x as the reference vector, and it deviates near and below the equator:
+
+.. image:: figures/Stereographic.png
+  :width: 400
+  :alt: Stereographic
+
+
+Converting from one convention and/or frame to another
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Use the :code:`transform_to` function::
+
+ pa_inertial =  PolarizationAngle(20*u.deg, source_direction, convention = 'IAU')
+
+ pa2_sc = pa.transform_to('RelativeX', attitude = Attitude.from_rotvec([0,10,0]*u.deg))
+
+ print(pa2_sc.angle.degree)
+
+Results in 161.95. Note that in addition to accounting for the difference polarization reference vector, it also transforms from one reference frame to another on the fly.
+
+
+.. automodule:: cosipy.polarization
+   :imported-members:
+   :members:
+   :undoc-members:
diff --git a/tests/polarization/__init__.py b/tests/polarization/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/polarization/test_conventions.py b/tests/polarization/test_conventions.py
new file mode 100644
index 00000000..04c6d7bd
--- /dev/null
+++ b/tests/polarization/test_conventions.py
@@ -0,0 +1,28 @@
+import numpy as np
+import pytest
+from astropy.coordinates import SkyCoord, Angle
+import astropy.units as u
+from cosipy.polarization import OrthographicConvention, StereographicConvention
+from scoords import SpacecraftFrame
+
+# Define common test data
+source_direction = SkyCoord(ra=-90*u.deg, dec=0*u.deg, frame='icrs') # -y
+
+def test_orthographic_projection_default():
+
+    ortho_convention = OrthographicConvention()
+    px, py = ortho_convention.get_basis(source_direction)
+
+    assert np.isclose(px.separation(SkyCoord(ra = 0*u.deg, dec = 90*u.deg)), 0)
+    assert np.isclose(py.separation(SkyCoord(ra = 0*u.deg, dec = 0*u.deg)), 0)
+
+def test_stereographic_projection_default():
+    
+    stereo_convention = StereographicConvention()
+    px, py = stereo_convention.get_basis(source_direction)
+    
+    assert np.isclose(px.separation(SkyCoord(lon = 0*u.deg, lat = 0*u.deg,
+                                             frame = SpacecraftFrame())), 0)
+    assert np.isclose(py.separation(SkyCoord(lon = 0*u.deg, lat = -90*u.deg,
+                                             frame = SpacecraftFrame())), 0)
+
diff --git a/tests/polarization/test_polarization_angle.py b/tests/polarization/test_polarization_angle.py
new file mode 100644
index 00000000..4b6b3641
--- /dev/null
+++ b/tests/polarization/test_polarization_angle.py
@@ -0,0 +1,39 @@
+from cosipy.polarization import OrthographicConvention, StereographicConvention, PolarizationAngle
+
+import numpy as np
+import pytest
+from astropy.coordinates import SkyCoord, Angle
+import astropy.units as u
+from cosipy.polarization import OrthographicConvention, StereographicConvention
+from scoords import SpacecraftFrame, Attitude
+
+# Define common test data
+source_direction = SkyCoord(ra = -36*u.deg, dec = 30*u.deg, frame = 'icrs')
+
+def test_pa_transformation():
+
+    pa = PolarizationAngle(20*u.deg, source_direction, convention = 'IAU')
+
+    pa2 = pa.transform_to(StereographicConvention(attitude = Attitude.identity()))
+    
+    assert np.isclose(pa2.angle, 56*u.deg)
+
+    pa2 = pa.transform_to('RelativeZ', attitude = Attitude.identity())
+
+    assert np.isclose(pa2.angle, 110*u.deg)
+
+    pa2 = pa.transform_to('RelativeZ', attitude = Attitude.from_rotvec([0,0,10]*u.deg))
+
+    assert np.isclose(pa2.angle, 110*u.deg)
+
+    pa2 = pa.transform_to('RelativeZ', attitude = Attitude.from_rotvec([30,0,0]*u.deg))
+
+    assert np.isclose(pa2.angle, 143.852403963*u.deg)
+
+    
+    pa2 = pa.transform_to('RelativeX', attitude = Attitude.identity())
+    
+    assert np.isclose(pa2.angle, 165.46460289540*u.deg)
+    
+    
+