diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index ea5853d55eeac0440eae8cae539898fe9030ac10..eb55e3c9957a5b9649c07fa17f58e4b60fc1da47 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -215,7 +215,7 @@ def convert_to_lal_binary_black_hole_parameters(parameters):
             converted_parameters['phi_jl'] = 0.0
             converted_parameters['phi_12'] = 0.0
 
-    for angle in ['tilt_1', 'tilt_2', 'iota']:
+    for angle in ['tilt_1', 'tilt_2', 'theta_jn']:
         cos_angle = str('cos_' + angle)
         if cos_angle in converted_parameters.keys():
             converted_parameters[angle] =\
@@ -827,7 +827,7 @@ def generate_component_spins(sample):
     Parameters
     ----------
     sample: A dictionary with the necessary spin conversion parameters:
-    'iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
+    'theta_jn', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
     'mass_2', 'reference_frequency', 'phase'
 
     Returns
@@ -837,7 +837,7 @@ def generate_component_spins(sample):
     """
     output_sample = sample.copy()
     spin_conversion_parameters =\
-        ['iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
+        ['theta_jn', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
          'mass_2', 'reference_frequency', 'phase']
     if all(key in output_sample.keys() for key in spin_conversion_parameters):
         output_sample['iota'], output_sample['spin_1x'],\
@@ -845,7 +845,7 @@ def generate_component_spins(sample):
             output_sample['spin_2x'], output_sample['spin_2y'],\
             output_sample['spin_2z'] =\
             transform_precessing_spins(
-                output_sample['iota'], output_sample['phi_jl'],
+                output_sample['theta_jn'], output_sample['phi_jl'],
                 output_sample['tilt_1'], output_sample['tilt_2'],
                 output_sample['phi_12'], output_sample['a_1'],
                 output_sample['a_2'],
diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py
index 2021a758aad7273035e2dcc956aefb092dc5f2f7..99c353b44a2dd4ebc4b73669a483b2a2438b4e88 100644
--- a/bilby/gw/prior.py
+++ b/bilby/gw/prior.py
@@ -232,7 +232,7 @@ class BBHPriorDict(PriorDict):
         spin_tilt_1_parameters = {'tilt_1', 'cos_tilt_1'}
         spin_tilt_2_parameters = {'tilt_2', 'cos_tilt_2'}
         spin_azimuth_parameters = {'phi_1', 'phi_2', 'phi_12', 'phi_jl'}
-        inclination_parameters = {'iota', 'cos_iota'}
+        inclination_parameters = {'theta_jn', 'cos_theta_jn'}
         distance_parameters = {'luminosity_distance', 'comoving_distance', 'redshift'}
 
         for independent_parameters, parameter_set in \
@@ -317,6 +317,8 @@ Prior._default_latex_labels = {
     'ra': '$\mathrm{RA}$',
     'iota': '$\iota$',
     'cos_iota': '$\cos\iota$',
+    'theta_jn': '$\\theta_{JN}$',
+    'cos_theta_jn': '$\cos\\theta_{JN}$',
     'psi': '$\psi$',
     'phase': '$\phi$',
     'geocent_time': '$t_c$',
diff --git a/bilby/gw/prior_files/GW150914.prior b/bilby/gw/prior_files/GW150914.prior
index 35b12524c498e5a21127e91c74472c9dffe48430..0eef13a9734f8f0e37dc6c578e3d3a0721e28dad 100644
--- a/bilby/gw/prior_files/GW150914.prior
+++ b/bilby/gw/prior_files/GW150914.prior
@@ -10,7 +10,7 @@ phi_jl =  Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=1e3, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
+theta_jn =  Sine(name='theta_jn')
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
 geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time', unit='$s$')
diff --git a/bilby/gw/prior_files/binary_black_holes.prior b/bilby/gw/prior_files/binary_black_holes.prior
index 5314635bd779f09098e3c8c2002f3367f046547d..e41cc73e2f765fbddd821693604ac15e302586f2 100644
--- a/bilby/gw/prior_files/binary_black_holes.prior
+++ b/bilby/gw/prior_files/binary_black_holes.prior
@@ -19,7 +19,7 @@ phi_jl =  Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/bilby/gw/prior_files/binary_neutron_stars.prior b/bilby/gw/prior_files/binary_neutron_stars.prior
index 1bc4d485c1cfdf89ddee3ab4b7b33d6cabd80654..08ec88f2ffb8bd6536563cb20723aabb589c7c13 100644
--- a/bilby/gw/prior_files/binary_neutron_stars.prior
+++ b/bilby/gw/prior_files/binary_neutron_stars.prior
@@ -13,8 +13,8 @@ chi_2 =  bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
 lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index b25aa5d4da232390416a39f10b1b3e37b156d8bd..6c16d3f8564c3fc11b8b57e04ed0d618e4aa8873 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -21,8 +21,8 @@ except ImportError:
 
 
 def lal_binary_black_hole(
-        frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, phi_12, a_2, tilt_2, phi_jl,
-        iota, phase, **kwargs):
+        frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
+        phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs):
     """ A Binary Black Hole waveform model using lalsimulation
 
     Parameters
@@ -40,15 +40,16 @@ def lal_binary_black_hole(
     tilt_1: float
         Primary tilt angle
     phi_12: float
-
+        Azimuthal angle between the two component spins
     a_2: float
         Dimensionless secondary spin magnitude
     tilt_2: float
         Secondary tilt angle
     phi_jl: float
-
-    iota: float
-        Orbital inclination
+        Azimuthal angle between the total binary angular momentum and the
+        orbital angular momentum
+    theta_jn: float
+        Angle between the total binary angular momentum and the line of sight
     phase: float
         The phase at coalescence
     kwargs: dict
@@ -59,7 +60,8 @@ def lal_binary_black_hole(
     dict: A dictionary with the plus and cross polarisation strain modes
     """
 
-    waveform_kwargs = dict(waveform_approximant='IMRPhenomPv2', reference_frequency=50.0,
+    waveform_kwargs = dict(waveform_approximant='IMRPhenomPv2',
+                           reference_frequency=50.0,
                            minimum_frequency=20.0)
     waveform_kwargs.update(kwargs)
     waveform_approximant = waveform_kwargs['waveform_approximant']
@@ -80,10 +82,11 @@ def lal_binary_black_hole(
         spin_2x = 0
         spin_2y = 0
         spin_2z = a_2
+        iota = theta_jn
     else:
         iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = (
             lalsim_SimInspiralTransformPrecessingNewInitialConditions(
-                iota, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
+                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
                 mass_2, reference_frequency, phase))
 
     longitude_ascending_nodes = 0.0
@@ -114,7 +117,8 @@ def lal_binary_black_hole(
 
 
 def lal_eccentric_binary_black_hole_no_spins(
-        frequency_array, mass_1, mass_2, eccentricity, luminosity_distance, iota, phase, **kwargs):
+        frequency_array, mass_1, mass_2, eccentricity, luminosity_distance,
+        theta_jn, phase, **kwargs):
     """ Eccentric binary black hole waveform model using lalsimulation (EccentricFD)
 
     Parameters
@@ -129,7 +133,7 @@ def lal_eccentric_binary_black_hole_no_spins(
         The orbital eccentricity of the system
     luminosity_distance: float
         The luminosity distance in megaparsec
-    iota: float
+    theta_jn: float
         Orbital inclination
     phase: float
         The phase at coalescence
@@ -161,6 +165,7 @@ def lal_eccentric_binary_black_hole_no_spins(
     spin_2x = 0.0
     spin_2y = 0.0
     spin_2z = 0.0
+    iota = theta_jn
 
     longitude_ascending_nodes = 0.0
     mean_per_ano = 0.0
@@ -246,7 +251,7 @@ def supernova_pca_model(
 
 def lal_binary_neutron_star(
         frequency_array, mass_1, mass_2, luminosity_distance, chi_1, chi_2,
-        iota, phase, lambda_1, lambda_2, **kwargs):
+        theta_jn, phase, lambda_1, lambda_2, **kwargs):
     """ A Binary Neutron Star waveform model using lalsimulation
 
     Parameters
@@ -263,7 +268,7 @@ def lal_binary_neutron_star(
         Dimensionless aligned spin
     chi_2: float
         Dimensionless aligned spin
-    iota: float
+    theta_jn: float
         Orbital inclination
     phase: float
         The phase at coalescence
@@ -308,6 +313,7 @@ def lal_binary_neutron_star(
     spin_2x = 0
     spin_2y = 0
     spin_2z = chi_2
+    iota = theta_jn
 
     longitude_ascending_nodes = 0.0
     eccentricity = 0.0
@@ -339,7 +345,7 @@ def lal_binary_neutron_star(
 
 
 def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
-        phi_12, a_2, tilt_2, phi_jl, iota, phase, **waveform_arguments):
+        phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **waveform_arguments):
     """
     See https://git.ligo.org/lscsoft/lalsuite/blob/master/lalsimulation/src/LALSimInspiral.c#L1460
 
@@ -365,7 +371,7 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
         Secondary tilt angle
     phi_jl: float
 
-    iota: float
+    theta_jn: float
         Orbital inclination
     phase: float
         The phase at coalescence
@@ -410,11 +416,12 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
         spin_2x = 0
         spin_2y = 0
         spin_2z = a_2
+        iota = theta_jn
     else:
         iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = \
             lalsim_SimInspiralTransformPrecessingNewInitialConditions(
-                iota, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
-                reference_frequency, phase)
+                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
+                mass_2, reference_frequency, phase)
 
     chi_1_l, chi_2_l, chi_p, theta_jn, alpha, phase_aligned, zeta =\
         lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
diff --git a/examples/injection_examples/australian_detector.py b/examples/injection_examples/australian_detector.py
index d770af182ce106e4b0e760a78413cb624bdf9aa7..3cd8c4aa94dd039e77644f9e69680892c3b6a147 100644
--- a/examples/injection_examples/australian_detector.py
+++ b/examples/injection_examples/australian_detector.py
@@ -59,7 +59,7 @@ interferometers.append(AusIFO)
 # signal at 4 Gpc
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=0.2108)
 
 
diff --git a/examples/injection_examples/binary_neutron_star_example.py b/examples/injection_examples/binary_neutron_star_example.py
index f1ad6d2d372ec18049630a2b0d92f199c9fa7ad1..cd540368845817ff1b7162d0cf63f5ba427fddf1 100644
--- a/examples/injection_examples/binary_neutron_star_example.py
+++ b/examples/injection_examples/binary_neutron_star_example.py
@@ -28,7 +28,7 @@ np.random.seed(88170235)
 # aligned spins of both black holes (chi_1, chi_2), etc.
 injection_parameters = dict(
     mass_1=1.5, mass_2=1.3, chi_1=0.02, chi_2=0.02, luminosity_distance=50.,
-    iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
+    theta_jn=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
     ra=1.375, dec=-1.2108, lambda_1=400, lambda_2=450)
 
 # Set the duration and sampling frequency of the data segment that we're going
@@ -66,7 +66,7 @@ interferometers.inject_signal(parameters=injection_parameters,
 # delta_lambda rather than mass_1, mass_2, lambda_1, and lambda_2.
 priors = bilby.gw.prior.BNSPriorDict()
 for key in ['psi', 'geocent_time', 'ra', 'dec', 'chi_1', 'chi_2',
-            'iota', 'luminosity_distance', 'phase']:
+            'theta_jn', 'luminosity_distance', 'phase']:
     priors[key] = injection_parameters[key]
 priors.pop('mass_1')
 priors.pop('mass_2')
diff --git a/examples/injection_examples/calibration_example.py b/examples/injection_examples/calibration_example.py
index ec3f3491e9ad61fda38512557d91a5345cf68f2d..ca58994696f4233ceb60d17bf500e144671e6c72 100644
--- a/examples/injection_examples/calibration_example.py
+++ b/examples/injection_examples/calibration_example.py
@@ -28,7 +28,7 @@ np.random.seed(88170235)
 # spins of both black holes (a, tilt, phi), etc.
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 # Fixed arguments passed into the source model
diff --git a/examples/injection_examples/change_sampled_parameters.py b/examples/injection_examples/change_sampled_parameters.py
index b2c848576bcde36985796e28c206c887411114d3..9f432ef52910fd125b1930ea90420e7eae0b6245 100644
--- a/examples/injection_examples/change_sampled_parameters.py
+++ b/examples/injection_examples/change_sampled_parameters.py
@@ -22,7 +22,7 @@ np.random.seed(151226)
 
 injection_parameters = dict(
     total_mass=66., mass_ratio=0.9, a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000, iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000, theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -63,8 +63,8 @@ priors['redshift'] = bilby.prior.Uniform(
 for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
             'ra', 'dec', 'geocent_time', 'phase']:
     priors[key] = injection_parameters[key]
-priors.pop('iota')
-priors['cos_iota'] = np.cos(injection_parameters['iota'])
+priors.pop('theta_jn')
+priors['cos_theta_jn'] = np.cos(injection_parameters['theta_jn'])
 print(priors)
 
 # Initialise GravitationalWaveTransient
diff --git a/examples/injection_examples/eccentric_inspiral.py b/examples/injection_examples/eccentric_inspiral.py
index ce4635c7b982cc4a44e0f7d976c1b300d0e1bc19..8d3053f353fb769384f25a17c3b1d406860fd88e 100644
--- a/examples/injection_examples/eccentric_inspiral.py
+++ b/examples/injection_examples/eccentric_inspiral.py
@@ -27,7 +27,7 @@ np.random.seed(150914)
 
 injection_parameters = dict(
     mass_1=35., mass_2=30., eccentricity=0.1, luminosity_distance=440.,
-    iota=0.4, psi=0.1, phase=1.2, geocent_time=1180002601.0, ra=45, dec=5.73)
+    theta_jn=0.4, psi=0.1, phase=1.2, geocent_time=1180002601.0, ra=45, dec=5.73)
 
 waveform_arguments = dict(waveform_approximant='EccentricFD',
                           reference_frequency=10., minimum_frequency=10.)
@@ -70,7 +70,7 @@ priors["luminosity_distance"] = bilby.gw.prior.UniformComovingVolume(
 priors["dec"] = bilby.core.prior.Cosine(name='dec')
 priors["ra"] = bilby.core.prior.Uniform(
     name='ra', minimum=0, maximum=2 * np.pi)
-priors["iota"] = bilby.core.prior.Sine(name='iota')
+priors["theta_jn"] = bilby.core.prior.Sine(name='theta_jn')
 priors["psi"] = bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi)
 priors["phase"] = bilby.core.prior.Uniform(
     name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/examples/injection_examples/fast_tutorial.py b/examples/injection_examples/fast_tutorial.py
index c472ba44b62af83b0858ae407bc4653cbb2341f2..6bfd0a0ed8bab100ee78d0fd03476a40b83b04fb 100644
--- a/examples/injection_examples/fast_tutorial.py
+++ b/examples/injection_examples/fast_tutorial.py
@@ -31,7 +31,7 @@ np.random.seed(88170235)
 # spins of both black holes (a, tilt, phi), etc.
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 # Fixed arguments passed into the source model
@@ -61,7 +61,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 # prior is a delta function at the true, injected value.  In reality, the
 # sampler implementation is smart enough to not sample any parameter that has
 # a delta-function prior.
-# The above list does *not* include mass_1, mass_2, iota and luminosity
+# The above list does *not* include mass_1, mass_2, theta_jn and luminosity
 # distance, which means those are the parameters that will be included in the
 # sampler.  If we do nothing, then the default priors get used.
 priors = bilby.gw.prior.BBHPriorDict()
diff --git a/examples/injection_examples/how_to_specify_the_prior.py b/examples/injection_examples/how_to_specify_the_prior.py
index 885ab328b955bc5794cf63bc194dc428972ac1fd..68081b4d2eddabec080e48005c6fc5536bb546e1 100644
--- a/examples/injection_examples/how_to_specify_the_prior.py
+++ b/examples/injection_examples/how_to_specify_the_prior.py
@@ -17,7 +17,7 @@ np.random.seed(151012)
 
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -41,7 +41,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 # This loads in a predefined set of priors for BBHs.
 priors = bilby.gw.prior.BBHPriorDict()
 # These parameters will not be sampled
-for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota', 'ra',
+for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'theta_jn', 'ra',
             'dec', 'geocent_time', 'psi']:
     priors[key] = injection_parameters[key]
 # We can make uniform distributions.
diff --git a/examples/injection_examples/marginalized_likelihood.py b/examples/injection_examples/marginalized_likelihood.py
index 06ef88b6731622146e984d5a750861ac4967bd56..6d84b1893c6b1f08419172c5e3ff77a9f5eba915 100644
--- a/examples/injection_examples/marginalized_likelihood.py
+++ b/examples/injection_examples/marginalized_likelihood.py
@@ -17,7 +17,7 @@ np.random.seed(170608)
 
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -40,7 +40,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 # Set up prior
 priors = bilby.gw.prior.BBHPriorDict()
 # These parameters will not be sampled
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'iota', 'ra',
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'theta_jn', 'ra',
             'dec']:
     priors[key] = injection_parameters[key]
 
diff --git a/examples/injection_examples/plot_skymap.py b/examples/injection_examples/plot_skymap.py
index 6453669ba0bc52396e892f1312351d69562e201a..9ebee7725fb0b913a2329c56622367f51033e129 100644
--- a/examples/injection_examples/plot_skymap.py
+++ b/examples/injection_examples/plot_skymap.py
@@ -12,7 +12,7 @@ outdir = 'outdir'
 label = 'plot_skymap'
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-0.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -33,7 +33,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 priors = bilby.gw.prior.BBHPriorDict()
 for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
             'mass_1', 'mass_2', 'phase', 'geocent_time', 'luminosity_distance',
-            'iota']:
+            'theta_jn']:
     priors[key] = injection_parameters[key]
 
 likelihood = bilby.gw.GravitationalWaveTransient(
diff --git a/examples/injection_examples/plot_time_domain_data.py b/examples/injection_examples/plot_time_domain_data.py
index f3d056b7a9a594a328a200c84ac334ec50775093..eee3c8557bd447714bd0fe9025f838928f9cd542 100644
--- a/examples/injection_examples/plot_time_domain_data.py
+++ b/examples/injection_examples/plot_time_domain_data.py
@@ -15,7 +15,7 @@ label = 'example'
 
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
diff --git a/examples/injection_examples/roq_example.py b/examples/injection_examples/roq_example.py
index cd8475b6c81bc84b4623f37db70a6653026246da..2ad8c99d4d020289e1e148f63829c1dbcc0c6560 100644
--- a/examples/injection_examples/roq_example.py
+++ b/examples/injection_examples/roq_example.py
@@ -33,7 +33,7 @@ sampling_frequency = 2048
 
 injection_parameters = dict(
     chirp_mass=36., mass_ratio=0.9, a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., iota=0.4, psi=0.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., theta_jn=0.4, psi=0.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -63,7 +63,7 @@ search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
     parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
 
 priors = bilby.gw.prior.BBHPriorDict()
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'iota', 'phase', 'psi', 'ra',
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra',
             'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
     priors[key] = injection_parameters[key]
 priors.pop('mass_1')
diff --git a/examples/injection_examples/standard_15d_cbc_tutorial.py b/examples/injection_examples/standard_15d_cbc_tutorial.py
index aa17219acd2da42e20ec878c9a2b13ee03feb831..822c11b1002f1bdb1021e5d0d7daef57dc16a259 100644
--- a/examples/injection_examples/standard_15d_cbc_tutorial.py
+++ b/examples/injection_examples/standard_15d_cbc_tutorial.py
@@ -27,7 +27,7 @@ np.random.seed(88170235)
 # spins of both black holes (a, tilt, phi), etc.
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 # Fixed arguments passed into the source model
diff --git a/examples/injection_examples/using_gwin.py b/examples/injection_examples/using_gwin.py
index 1539195e60cf0e68f1b3d565f35161e91fbce84b..4d9c583ac6fa35f81b038232c566d2fbada5ad8a 100644
--- a/examples/injection_examples/using_gwin.py
+++ b/examples/injection_examples/using_gwin.py
@@ -29,7 +29,7 @@ label = 'using_gwin'
 # Search priors
 priors = dict()
 priors['distance'] = bilby.core.prior.Uniform(500, 2000, 'distance')
-priors['polarization'] = bilby.core.prior.Uniform(0, np.pi, 'iota')
+priors['polarization'] = bilby.core.prior.Uniform(0, np.pi, 'theta_jn')
 
 # Data variables
 seglen = 4
diff --git a/test/gw_likelihood_test.py b/test/gw_likelihood_test.py
index ebe4fe275a5e49350756cec44d6c1c5a30de91be..632aadac94536dcd425e7b2914ef1f63bcf17a94 100644
--- a/test/gw_likelihood_test.py
+++ b/test/gw_likelihood_test.py
@@ -10,7 +10,7 @@ class TestBasicGWTransient(unittest.TestCase):
         np.random.seed(500)
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
         self.interferometers = bilby.gw.detector.InterferometerList(['H1'])
@@ -73,7 +73,7 @@ class TestGWTransient(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
         self.interferometers = bilby.gw.detector.InterferometerList(['H1'])
@@ -143,7 +143,7 @@ class TestTimeMarginalization(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259640, ra=1.375,
             dec=-1.2108)
 
@@ -242,7 +242,7 @@ class TestMarginalizedLikelihood(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
 
@@ -306,7 +306,7 @@ class TestPhaseMarginalization(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
 
@@ -369,7 +369,7 @@ class TestTimePhaseMarginalization(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
 
@@ -472,9 +472,9 @@ class TestROQLikelihood(unittest.TestCase):
         fnodes_quadratic = np.load(fnodes_quadratic_file).T
 
         self.test_parameters = dict(
-            mass_1=36.0, mass_2=36.0, a_1=0.0, a_2=0.0, tilt_1=0.0,
-            tilt_2=0.0, phi_12=1.7, phi_jl=0.3, luminosity_distance=5000.,
-            iota=0.4, psi=0.659, phase=1.3, geocent_time=1.2, ra=1.3, dec=-1.2)
+            mass_1=36.0, mass_2=36.0, a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=5000., theta_jn=0.4,
+            psi=0.659, phase=1.3, geocent_time=1.2, ra=1.3, dec=-1.2)
 
         ifos = bilby.gw.detector.InterferometerList(['H1'])
         ifos.set_strain_data_from_power_spectral_densities(
diff --git a/test/gw_prior_test.py b/test/gw_prior_test.py
index 5679d80000ea04381d2658df040387bce7714938..1eac2e76ca9311323edc6033668611cee5780a31 100644
--- a/test/gw_prior_test.py
+++ b/test/gw_prior_test.py
@@ -46,13 +46,13 @@ class TestBBHPriorDict(unittest.TestCase):
 
     def test_redundant_priors_not_in_dict_before(self):
         for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
-                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_iota',
+                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_theta_jn',
                       'comoving_distance', 'redshift']:
             self.assertTrue(self.bbh_prior_dict.test_redundancy(prior))
 
     def test_redundant_priors_already_in_dict(self):
         for prior in ['mass_1', 'mass_2', 'tilt_1', 'tilt_2',
-                      'phi_1', 'phi_2', 'iota', 'luminosity_distance']:
+                      'phi_1', 'phi_2', 'theta_jn', 'luminosity_distance']:
             self.assertTrue(self.bbh_prior_dict.test_redundancy(prior))
 
     def test_correct_not_redundant_priors_masses(self):
@@ -80,8 +80,8 @@ class TestBBHPriorDict(unittest.TestCase):
             self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
 
     def test_correct_not_redundant_priors_inclination(self):
-        del self.bbh_prior_dict['iota']
-        for prior in ['iota', 'cos_iota']:
+        del self.bbh_prior_dict['theta_jn']
+        for prior in ['theta_jn', 'cos_theta_jn']:
             self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
 
     def test_correct_not_redundant_priors_distance(self):
@@ -96,7 +96,7 @@ class TestBBHPriorDict(unittest.TestCase):
     def test_test_has_redundant_priors(self):
         self.assertFalse(self.bbh_prior_dict.test_has_redundant_keys())
         for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
-                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_iota',
+                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_theta_jn',
                       'comoving_distance', 'redshift']:
             self.bbh_prior_dict[prior] = 0
             self.assertTrue(self.bbh_prior_dict.test_has_redundant_keys())
@@ -139,13 +139,13 @@ class TestBNSPriorDict(unittest.TestCase):
 
     def test_redundant_priors_not_in_dict_before(self):
         for prior in ['chirp_mass', 'total_mass', 'mass_ratio',
-                      'symmetric_mass_ratio', 'cos_iota', 'comoving_distance',
+                      'symmetric_mass_ratio', 'cos_theta_jn', 'comoving_distance',
                       'redshift', 'lambda_tilde', 'delta_lambda']:
             self.assertTrue(self.bns_prior_dict.test_redundancy(prior))
 
     def test_redundant_priors_already_in_dict(self):
         for prior in ['mass_1', 'mass_2', 'chi_1', 'chi_2',
-                      'iota', 'luminosity_distance',
+                      'theta_jn', 'luminosity_distance',
                       'lambda_1', 'lambda_2']:
             self.assertTrue(self.bns_prior_dict.test_redundancy(prior))
 
@@ -159,8 +159,8 @@ class TestBNSPriorDict(unittest.TestCase):
         self.assertFalse(self.bns_prior_dict.test_redundancy('chi_2'))
 
     def test_correct_not_redundant_priors_inclination(self):
-        del self.bns_prior_dict['iota']
-        for prior in ['iota', 'cos_iota']:
+        del self.bns_prior_dict['theta_jn']
+        for prior in ['theta_jn', 'cos_theta_jn']:
             self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
 
     def test_correct_not_redundant_priors_distance(self):
@@ -180,7 +180,7 @@ class TestBNSPriorDict(unittest.TestCase):
     def test_test_has_redundant_priors(self):
         self.assertFalse(self.bns_prior_dict.test_has_redundant_keys())
         for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
-                      'cos_iota', 'comoving_distance', 'redshift']:
+                      'cos_theta_jn', 'comoving_distance', 'redshift']:
             self.bns_prior_dict[prior] = 0
             self.assertTrue(self.bns_prior_dict.test_has_redundant_keys())
             del self.bns_prior_dict[prior]
diff --git a/test/make_standard_data.py b/test/make_standard_data.py
index fcab591cd11b14313fa87aa7d74483f58ac9823f..458ce69b4652ba7cd90fd90a3c746d055cc0b60f 100644
--- a/test/make_standard_data.py
+++ b/test/make_standard_data.py
@@ -22,7 +22,7 @@ simulation_parameters = dict(
     phi_12=0.,
     phi_jl=0.,
     luminosity_distance=100.,
-    iota=0.4,
+    theta_jn=0.4,
     phase=1.3,
     waveform_approximant='IMRPhenomPv2',
     reference_frequency=50.,
diff --git a/test/prior_files/binary_black_holes.prior b/test/prior_files/binary_black_holes.prior
index 5314635bd779f09098e3c8c2002f3367f046547d..e41cc73e2f765fbddd821693604ac15e302586f2 100644
--- a/test/prior_files/binary_black_holes.prior
+++ b/test/prior_files/binary_black_holes.prior
@@ -19,7 +19,7 @@ phi_jl =  Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/test/prior_files/binary_neutron_stars.prior b/test/prior_files/binary_neutron_stars.prior
index 1bc4d485c1cfdf89ddee3ab4b7b33d6cabd80654..08ec88f2ffb8bd6536563cb20723aabb589c7c13 100644
--- a/test/prior_files/binary_neutron_stars.prior
+++ b/test/prior_files/binary_neutron_stars.prior
@@ -13,8 +13,8 @@ chi_2 =  bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
 lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
diff --git a/test/prior_test.py b/test/prior_test.py
index 4b85a66c625b43f1194fa79213278e807b65f399..d73d8d5b9001debc72c7cc7379c864a4619f90d1 100644
--- a/test/prior_test.py
+++ b/test/prior_test.py
@@ -380,7 +380,7 @@ class TestPriorDict(unittest.TestCase):
             dec=bilby.core.prior.Cosine(name='dec'),
             ra=bilby.core.prior.Uniform(
                 name='ra', minimum=0, maximum=2 * np.pi),
-            iota=bilby.core.prior.Sine(name='iota'),
+            theta_jn=bilby.core.prior.Sine(name='theta_jn'),
             psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
             phase=bilby.core.prior.Uniform(
                 name='phase', minimum=0, maximum=2 * np.pi)
@@ -436,7 +436,7 @@ class TestPriorDict(unittest.TestCase):
             dec=bilby.core.prior.Cosine(name='dec'),
             ra=bilby.core.prior.Uniform(
                 name='ra', minimum=0, maximum=2 * np.pi),
-            iota=bilby.core.prior.Sine(name='iota'),
+            theta_jn=bilby.core.prior.Sine(name='theta_jn'),
             psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
             phase=bilby.core.prior.Uniform(
                 name='phase', minimum=0, maximum=2 * np.pi)