Deprecate simulator.models.lorenzian and rename to lorentzian
In simulator/models.py there is a function called lorenzian(). This function is addressed normally by passing its name as a string to the simulator. The problem is the function name has a typo in it: It should be lorentzian.

I suggest to deprecate the rename lorenzian to lorentzian and introduce a new function lorenzian that (a) prints a deprecation warning  and (b) calls lorentzian.

But if you're more fearless than me you can of course also simply rename.
diff --git a/stingray/simulator/models.py b/stingray/simulator/models.py
index 23b876d4..51eda4eb 100644
--- a/stingray/simulator/models.py
+++ b/stingray/simulator/models.py
@@ -1,10 +1,12 @@
 from __future__ import division, print_function, absolute_import
 from astropy.modeling.models import custom_model
 
+
+#TODO: Add Jacobian
 @custom_model
 def GeneralizedLorentz1D(x, x_0=1., fwhm=1., value=1., power_coeff=1.):
     """
-    Generalized lorenzian function, 
+    Generalized Lorentzian function,
     implemented using astropy.modeling.models custom model
 
     Parameters
@@ -24,8 +26,10 @@ def GeneralizedLorentz1D(x, x_0=1., fwhm=1., value=1., power_coeff=1.):
         generalized Lorentzian psd model
     """
     assert power_coeff > 0., "The power coefficient should be greater than zero."
-    return value * (fwhm / 2)**power_coeff * 1./(abs(x - x_0)**power_coeff + (fwhm / 2)**power_coeff) 
-    
+    return value * (fwhm / 2)**power_coeff * 1./(abs(x - x_0)**power_coeff + (fwhm / 2)**power_coeff)
+
+
+#TODO: Add Jacobian
 @custom_model
 def SmoothBrokenPowerLaw(x, norm=1., gamma_low=1., gamma_high=1., break_freq=1.):
     """
@@ -34,7 +38,6 @@ def SmoothBrokenPowerLaw(x, norm=1., gamma_low=1., gamma_high=1., break_freq=1.)
 
     Parameters
     ----------
-
     x: numpy.ndarray
         non-zero frequencies
 
@@ -51,9 +54,9 @@ def SmoothBrokenPowerLaw(x, norm=1., gamma_low=1., gamma_high=1., break_freq=1.)
     return norm * x**(-gamma_low) / (1. + (x / break_freq)**2)**(-(gamma_low - gamma_high) / 2)
     
 
-def lorenzian(x, p):
+def generalized_lorentzian(x, p):
     """
-    Generalized lorenzian function.
+    Generalized Lorentzian function.
 
     Parameters
     ----------
@@ -70,11 +73,12 @@ def lorenzian(x, p):
     Returns
     -------
     model: numpy.ndarray
-        generalized lorenzian psd model
+        generalized lorentzian psd model
     """
 
     assert p[3] > 0., "The power coefficient should be greater than zero."
-    return p[2] * (p[1] / 2)**p[3] * 1./(abs(x - p[0])**p[3] + (p[1] / 2)**p[3]) 
+    return p[2] * (p[1] / 2)**p[3] * 1./(abs(x - p[0])**p[3] + (p[1] / 2)**p[3])
+
 
 def smoothbknpo(x, p):
     """
diff --git a/stingray/simulator/tests/test_simulator.py b/stingray/simulator/tests/test_simulator.py
index 5130f3e1..7bed5536 100644
--- a/stingray/simulator/tests/test_simulator.py
+++ b/stingray/simulator/tests/test_simulator.py
@@ -55,14 +55,17 @@ class TestSimulator(object):
         """
         Simulate an energy channel.
         """
-        self.simulator.simulate_channel('3.5-4.5', 'lorenzian', [1, 2, 3, 4])
+        self.simulator.simulate_channel('3.5-4.5',
+                                        'generalized_lorentzian', [1, 2, 3, 4])
         self.simulator.delete_channel('3.5-4.5')
 
     def test_simulate_channel_odd(self):
         """
         Simulate an energy channel.
         """
-        self.simulator_odd.simulate_channel('3.5-4.5', 'lorenzian', [1, 2, 3, 4])
+        self.simulator_odd.simulate_channel('3.5-4.5',
+                                            'generalized_lorentzian',
+                                            [1, 2, 3, 4])
         self.simulator_odd.delete_channel('3.5-4.5')
 
     def test_incorrect_simulate_channel(self):
@@ -164,33 +167,35 @@ class TestSimulator(object):
         s = np.random.rand(1024)
         assert len(self.simulator.simulate(s)), 1024
 
-    def test_simulate_lorenzian(self):
+    def test_simulate_lorentzian(self):
         """
-        Simulate light curve using lorenzian model.
+        Simulate light curve using lorentzian model.
         """
-        assert len(self.simulator.simulate('lorenzian', [1, 2, 3, 4])), 1024
+        assert len(self.simulator.simulate('generalized_lorentzian',
+                                           [1, 2, 3, 4])), 1024
 
-    def test_simulate_lorenzian_odd(self):
+    def test_simulate_lorentzian_odd(self):
         """
-        Simulate light curve using lorenzian model.
+        Simulate light curve using lorentzian model.
         """
-        assert len(self.simulator_odd.simulate('lorenzian',
+        assert len(self.simulator_odd.simulate('generalized_lorentzian',
                                                [1, 2, 3, 4])), 1024
 
-    def test_compare_lorenzian(self):
+    def test_compare_lorentzian(self):
         """
-        Compare simulated lorenzian spectrum with original spectrum.
+        Compare simulated lorentzian spectrum with original spectrum.
         """
         N, red_noise, dt = 1024, 10, 1
 
         self.simulator = simulator.Simulator(N=N, dt=dt, mean=0.1,
                                              rms=0.4, red_noise=red_noise)
-        lc = [self.simulator.simulate('lorenzian', [0.3, 0.9, 0.6, 0.5])
+        lc = [self.simulator.simulate('generalized_lorentzian',
+                                      [0.3, 0.9, 0.6, 0.5])
               for i in range(1, 30)]
         simulated = self.simulator.powerspectrum(lc, lc[0].tseg)
 
         w = np.fft.rfftfreq(N, d=dt)[1:]
-        actual = models.lorenzian(w, [0.3, 0.9, 0.6, 0.5])[:-1]
+        actual = models.generalized_lorentzian(w, [0.3, 0.9, 0.6, 0.5])[:-1]
 
         actual_prob = actual/float(sum(actual))
         simulated_prob = simulated/float(sum(simulated))
@@ -259,14 +264,18 @@ class TestSimulator(object):
         Simulate a light curve using SmoothBrokenPowerLaw model
         called as a string
         """
-        assert len(self.simulator.simulate('SmoothBrokenPowerLaw', {'norm':1., 'gamma_low':1., 'gamma_high':2., 'break_freq':1.})), 1024
+        assert len(
+            self.simulator.simulate('SmoothBrokenPowerLaw',
+                                    {'norm':1., 'gamma_low':1.,
+                                     'gamma_high':2., 'break_freq':1.})), 1024
     
     def test_simulate_SmoothBrokenPowerLaw(self):
         """
         Simulate a light curve using SmoothBrokenPowerLaw model
         called as a astropy.modeling.Model class
         """
-        mod = models.SmoothBrokenPowerLaw(norm=1., gamma_low=1., gamma_high=2., break_freq=1.)
+        mod = models.SmoothBrokenPowerLaw(norm=1., gamma_low=1., gamma_high=2.,
+                                          break_freq=1.)
         assert len(self.simulator.simulate(mod)), 1024
         
         
@@ -275,7 +284,8 @@ class TestSimulator(object):
         Simulate a light curve using a generic model
         called as a astropy.modeling.Model class
         """
-        mod = astropy.modeling.models.Gaussian1D(amplitude=10., mean=1., stddev=2.)
+        mod = astropy.modeling.models.Gaussian1D(amplitude=10., mean=1.,
+                                                 stddev=2.)
         assert len(self.simulator.simulate(mod)), 1024
 
     def test_simulate_generic_model_odd(self):
@@ -298,8 +308,11 @@ class TestSimulator(object):
         m = 30000.
         
         self.simulator = simulator.Simulator(N=N, mean=m, dt=dt)
-        smoothbknpo = models.SmoothBrokenPowerLaw(norm=1., gamma_low=1., gamma_high=2., break_freq=1.)
-        lorentzian = models.GeneralizedLorentz1D(x_0=10, fwhm=1., value=10., power_coeff=2.)
+        smoothbknpo = \\
+            models.SmoothBrokenPowerLaw(norm=1., gamma_low=1., gamma_high=2.,
+                                        break_freq=1.)
+        lorentzian = models.GeneralizedLorentz1D(x_0=10, fwhm=1., value=10.,
+                                                 power_coeff=2.)
         myModel = smoothbknpo + lorentzian
         
         lc = [self.simulator.simulate(myModel) for i in range(1, 50)]
@@ -312,7 +325,8 @@ class TestSimulator(object):
         actual_prob = actual/float(sum(actual))
         simulated_prob = simulated/float(sum(simulated))
 
-        assert np.all(np.abs(actual_prob - simulated_prob) < 3*np.sqrt(actual_prob))
+        assert np.all(np.abs(actual_prob - simulated_prob) <
+                      3*np.sqrt(actual_prob))
         
         
     def test_simulate_wrong_model(self):
['stingray/simulator/models.py', 'stingray/simulator/tests/test_simulator.py']