diff --git a/doc/conf.py b/doc/conf.py
index f9ed3c4..1790f75 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -23,6 +23,9 @@
 
 # The full version, including alpha/beta/rc tags
 import dsgp4
+import sys
+import os
+sys.path.insert(0, os.path.abspath('../dsgp4'))
 
 release = dsgp4.__version__
 
@@ -32,7 +35,11 @@
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
-extensions = ["myst_nb", "sphinx.ext.autodoc", "sphinx.ext.doctest", "sphinx.ext.intersphinx", "sphinx.ext.autosummary"]
+extensions = ["myst_nb", "sphinx.ext.autodoc", "sphinx.ext.doctest", "sphinx.ext.intersphinx", "sphinx.ext.autosummary","sphinx.ext.napoleon"]
+
+autosummary_generate = True
+napoleon_google_docstring = True
+napoleon_numpy_docstring = False
 
 intersphinx_mapping = {
     "numpy": ("https://numpy.org/doc/stable/", None),
diff --git a/dsgp4/newton_method.py b/dsgp4/newton_method.py
index 024ecf0..87928ad 100644
--- a/dsgp4/newton_method.py
+++ b/dsgp4/newton_method.py
@@ -34,8 +34,8 @@ def initial_guess(tle_0, time_mjd, target_state=None):
     Args:
         - tle_0 (``dsgp4.tle.TLE``): starting TLE at time0
         - new_date (``datetime.datetime``): new date of the TLE, as a datetime object
-        - target_state (``torch.tensor``): a 2x3 tensor for the position and velocity of the state. If None, then this is computed
-                                           by propagating `tle_0` at the `new_date`.
+        - target_state (``torch.tensor``): a 2x3 tensor for the position and velocity of the state. If None, then this is computed 
+                                            by propagating `tle_0` at the `new_date`.
 
     Returns:
         - y0 (``torch.tensor``): initial guess for the TLE elements. In particular, y0 contains
@@ -152,21 +152,20 @@ def update_TLE(old_tle,y0):
 
 def newton_method(tle_0, time_mjd, target_state=None, new_tol=1e-12,max_iter=50):
     """
-    This method performs Newton method starting from an initial TLE and a given propagation time. The objective
+    This method performs Newton method starting from an initial TLE and a given propagation time. The objective 
     is to find a TLE that accurately reconstructs the propagated state, at observation time.
 
     Args:
         - tle_0 (``dsgp4.tle.TLE``): starting TLE (i.e., TLE at a given initial time)
-        - new_date (``datetime.datetime``): time (as a datetime object) at which we want the state to be propagated, and we want the TLE at that time
-        - target_state (``torch.tensor``): 2x3 tensor representing target state. If None, then this is directly computed by propagating the TLE at `new_date`
-        - new_tol (``float``): newton tolerance
-        - max_iter (``int``): maximum iterations for Newton's method
+        - time_mjd (``float``): time at which we want to propagate the TLE (in Modified Julian Date)
+        - target_state (``torch.tensor``): a 2x3 tensor for the position and velocity of the state. If None, then this is computed by propagating `tle_0` at the `new_date`.
+        - new_tol (``float``): tolerance for the Newton method
+        - max_iter (``int``): maximum number of iterations for the Newton method
 
     Returns:
-        - tle (``dsgp4.tle.TLE``): found TLE
-        - y (``torch.tensor``): returns the solution that satisfied F(y)=0 (with a given tolerance)
-                                or (in case no convergence is reached within the tolerance) the best
-                                guess found in `max_iter` iterations
+        - next_tle (``dsgp4.tle.TLE``): TLE corresponding to the propagated state
+        - y0 (``torch.tensor``): initial guess for the TLE elements. In particular, y0 contains the following elements (see SGP4 for a thorough description of these parameters):
+
     """
     i=0
     tol=1e9
diff --git a/dsgp4/plot.py b/dsgp4/plot.py
index 3a82ba8..7f8e129 100644
--- a/dsgp4/plot.py
+++ b/dsgp4/plot.py
@@ -11,7 +11,7 @@ def plot_orbit(states, r_earth=6378.137, elevation_azimuth=None, ax=None, *args,
                                      position (in km) and the second the spacecraft velocity (in km/s). Reference frame is TEME.
         - r_earth (``float``): Earth radius in km (used for the plot). Default value is 6378.137 km.
         - elevation_azimuth (``tuple``): tuple of two floats, representing the elevation and azimuth angles of the plot. If None, the default values are used.
-        - ax (``matplotlib.axes._subplots.Axes3DSubplot``): 3D axis object
+        - ax (``matplotlib.axes._subplots.Axes3DSubplot``): 3D axis object. If None, a new figure is created.
         - *args, **kwargs: additional arguments to be passed to the plot function
     
     Returns:
@@ -53,20 +53,30 @@ def plot_orbit(states, r_earth=6378.137, elevation_azimuth=None, ax=None, *args,
     plt.tight_layout()
     return ax
 
-def plot_tles(tles, file_name=None, figsize = (36,18), show=True, axs=None, return_axs=False, log_yscale=False, *args, **kwargs):
+def plot_tles(tles, 
+              file_name=None, 
+              figsize = (36,18), 
+              show=True, 
+              axs=None, 
+              return_axs=False, 
+              log_yscale=False, 
+              *args, 
+              **kwargs):
     """
     This function takes a list of tles as input and plots the histograms of some of their elements.
-
-    Inputs
-    ----------
-    - tles (`list`): list of tles, where each element is a `dsgp4.tle.TLE` object.
-    - save_fig (`bool`): boolean variable, if True, figure is saved to a file.
-    - file_name (`str`): name of the file (including path) where the plot is to be saved.
-    - figsize (`tuple`): figure size.
-
-    Outputs
-    ----------
-    - ax (`numpy.ndarray`): array of AxesSubplot objects
+    
+    Args:    
+        - tles (``list``): list of tles, where each element is a ``dsgp4.tle.TLE`` object.
+        - file_name (``str``): name of the file (including path) where the plot is to be saved.
+        - figsize (``tuple``): figure size.
+        - show (``bool``): if True, the plot is shown.
+        - axs (``numpy.ndarray``): array of AxesSubplot objects.
+        - return_axs (``bool``): if True, the function returns the array of AxesSubplot objects.
+        - log_yscale (``bool``): if True, the y-scale is logarithmic.
+        - *args, **kwargs: additional arguments to be passed to the hist function.
+    
+    Returns:    
+        - ax (``numpy.ndarray``): array of AxesSubplot objects
     """
     #I collect all the six variables from the TLEs:
     mean_motion, eccentricity, inclination, argument_of_perigee, raan, b_star, mean_anomaly, mean_motion_first_derivative, mean_motion_second_derivative = [], [], [], [], [], [], [], [], []
diff --git a/dsgp4/tle.py b/dsgp4/tle.py
index 1d4dbde..cea6e28 100644
--- a/dsgp4/tle.py
+++ b/dsgp4/tle.py
@@ -28,10 +28,10 @@ def read_satellite_catalog_number(string):
     the corresponding satellite catalog number.
 
     Args:
-        string (`str`): string line
+        string (``str``): string line
     
     Returns:
-        `int`: satellite catalog number
+        ``int``: satellite catalog number
     """
     if not string[0].isalpha():
         return int(string)
@@ -47,11 +47,11 @@ def load_from_lines(lines, opsmode='i'):
     """
     This function takes a TLE as a list of strings and returns both itself and its dictionary representation.
     Args:
-        lines (`list`): TLE data in the form of a list
-        opsmode (`str`): operation mode, either 'i' or 'a'
+        lines (``list``): TLE data in the form of a list
+        opsmode (``str``): operation mode, either 'i' or 'a'
     Returns:
-        `list`: TLE data in the form of a list
-        `dict`: TLE data in the form of a dictionary
+        ``list``: TLE data in the form of a list
+        ``dict``: TLE data in the form of a dictionary
     """
     if isinstance(lines, str):
         lines = util.get_non_empty_lines(lines)
diff --git a/dsgp4/util.py b/dsgp4/util.py
index 9933b9d..7ea94ca 100644
--- a/dsgp4/util.py
+++ b/dsgp4/util.py
@@ -40,9 +40,9 @@ def propagate_batch(tles, tsinces, initialized=True):
     This function takes a list of TLEs and a tensor of times (which must be of same length), and returns the corresponding states.
 
     Args:
-        - tle_sat (``list``): list of TLE objects or TLE object to be propagated
-        - tsince (``torch.tensor``): propagation time in minutes (it has to be a tensor of the same size of the list, in case a list of TLEs is passed)
-        - initialized (``bool``): whether the TLEs have been initialized or not (default: True)
+        - tles (``list`` of ``dsgp4.tle.TLE``): list of TLE objects to be propagated
+        - tsinces (``torch.tensor``): propagation times in minutes (it has to be a tensor of the same size of the list of TLEs)
+        - initialized (``bool``): whether the TLEs have been initialized or not (default: True
 
     Returns:
         - state (``torch.tensor``): (Nx2x3) tensor representing position and velocity in km and km/s, where the first dimension is the batch size.
@@ -63,8 +63,8 @@ def propagate(tle, tsinces, initialized=True):
     In the second case, the length of the list of TLEs must be equal to the length of the tensor of times.
 
     Args:
-        - tle_sat (``list`` or ``dsgp4.tle.TLE``): list of TLE objects or TLE object to be propagated
-        - tsince (``torch.tensor``): propagation time in minutes (it has to be a tensor of the same size of the list, in case a list of TLEs is passed)
+        - tle (``dsgp4.tle.TLE`` or ``list`` of ``dsgp4.tle.TLE``): TLE object or list of TLE objects to be propagated
+        - tsinces (``torch.tensor``): propagation times in minutes
         - initialized (``bool``): whether the TLEs have been initialized or not (default: True)
 
     Returns:
@@ -167,9 +167,28 @@ def initialize_tle(tles,
         return tle_elements
 
 def from_year_day_to_date(y,d):
+    """
+    Converts a year and day of the year to a date.
+
+    Args:
+        - y (``int``): year
+        - d (``int``): day of the year
+
+    Returns:
+        - ``datetime.datetime``: date
+    """
     return (datetime.datetime(y, 1, 1) + datetime.timedelta(d - 1))
 
 def gstime(jdut1):
+    """
+    This function computes the Greenwich Sidereal Time (GST) at the given Julian Date (UT1).
+
+    Args:
+        - jdut1 (``float``): Julian Date (UT1)
+
+    Returns:
+        - ``float``: Greenwich Sidereal Time (GST)
+    """
     tut1 = (jdut1 - 2451545.0) / 36525.0
     temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 + \
          (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841  #  sec
@@ -180,6 +199,15 @@ def gstime(jdut1):
     return temp
 
 def clone_w_grad(y):
+    """
+    This function takes a tensor and returns a copy of it with gradients.
+
+    Args:
+        - y (``torch.tensor``): tensor to clone
+
+    Returns:
+        - ``torch.tensor``: tensor with gradients
+    """
     return y.clone().detach().requires_grad_(True)
 
 def jday(year, mon, day, hr, minute, sec):
@@ -187,15 +215,15 @@ def jday(year, mon, day, hr, minute, sec):
     Converts a date and time to a Julian Date. The Julian Date is the number of days since noon on January 1st, 4713 BC.
 
     Args:
-        year (`int`): year
-        mon (`int`): month
-        day (`int`): day
-        hr (`int`): hour
-        minute (`int`): minute
-        sec (`float`): second
+        - year (``int``): year
+        - mon (``int``): month
+        - day (``int``): day
+        - hr (``int``): hour
+        - minute (``int``): minute
+        - sec (``float``): second
 
     Returns:
-        `float`: Julian Date
+        - ``tuple``: Julian Date as integer and fractional part of the day
     """
     jd=(367.0 * year -
             7.0 * (year + ((mon + 9.0) // 12.0)) * 0.25 // 1.0 +
@@ -209,10 +237,10 @@ def invjday(jd):
     Converts a Julian Date to a date and time. The Julian Date is the number of days since noon on January 1st, 4713 BC.
 
     Args:
-        jd (`float`): Julian Date
+        - jd (``float``): Julian Date
 
     Returns:
-        `tuple`: (year, month, day, hour, minute, second)
+        - ``tuple``: year, month, day, hour, minute, second
     """
     temp    = jd - 2415019.5
     tu      = temp / 365.25
@@ -232,11 +260,11 @@ def days2mdhms(year, fractional_day):
     Converts a number of days to months, days, hours, minutes, and seconds.
 
     Args:
-        year (`int`): year
-        fractional_day (`float`): number of days
+        - year (``int``): year
+        - fractional_day (``float``): number of days
 
     Returns:
-        `tuple`: (month, day, hour, minute, second)
+        - ``tuple``: month, day, hour, minute, second
     """
     d=datetime.timedelta(days=fractional_day)
     datetime_obj=datetime.datetime(year-1,12,31)+d
@@ -247,10 +275,10 @@ def from_string_to_datetime(string):
     Converts a string to a datetime object.
 
     Args:
-        string (`str`): string to convert
+        - string (``str``): string to convert
 
     Returns:
-        `datetime.datetime`: datetime object
+        - ``datetime.datetime``: datetime object
     """
     if string.find('.')!=-1:
         return datetime.datetime.strptime(string, '%Y-%m-%d %H:%M:%S.%f')
@@ -262,10 +290,10 @@ def from_mjd_to_epoch_days_after_1_jan(mjd_date):
     Converts a Modified Julian Date to the number of days after 1 Jan 2000.
 
     Args:
-        mjd_date (`float`): Modified Julian Date
+        - mjd_date (``float``): Modified Julian Date
 
     Returns:
-        `float`: number of days after 1 Jan 2000
+        - ``float``: number of days after 1 Jan 2000
     """
     d = from_mjd_to_datetime(mjd_date)
     dd = d - datetime.datetime(d.year-1, 12, 31)
@@ -278,10 +306,10 @@ def from_mjd_to_datetime(mjd_date):
     Converts a Modified Julian Date to a datetime object. The Modified Julian Date is the number of days since midnight on November 17, 1858.
 
     Args:
-        mjd_date (`float`): Modified Julian Date
+        - mjd_date (``float``): Modified Julian Date
 
     Returns:
-        `datetime.datetime`: datetime object
+        - ``datetime.datetime``: datetime
     """
     jd_date=mjd_date+2400000.5
     return from_jd_to_datetime(jd_date)
@@ -291,10 +319,10 @@ def from_jd_to_datetime(jd_date):
     Converts a Julian Date to a datetime object. The Julian Date is the number of days since noon on January 1st, 4713 BC.
 
     Args:
-        jd_date (`float`): Julian Date
+        - jd_date (``float``): Julian Date
 
     Returns:
-        `datetime.datetime`: datetime object
+        - ``datetime.datetime``: datetime
     """
     year, month, day, hour, minute, seconds=invjday(jd_date)
     e_1=datetime.datetime(year=int(year), month=int(month), day=int(day), hour=int(hour), minute=int(minute), second=0)
@@ -305,10 +333,10 @@ def get_non_empty_lines(lines):
     This function returns the non-empty lines of a list of lines.
 
     Args:
-        lines (`list`): list of lines
+        - lines (``list``): list of lines
 
     Returns:
-        `list`: non-empty lines
+        - ``list``: list of non-empty lines
     """
     if not isinstance(lines, str):
         raise ValueError('Expecting a string')
@@ -321,10 +349,10 @@ def from_datetime_to_fractional_day(datetime_object):
     Converts a datetime object to a fractional day. The fractional day is the number of days since the beginning of the year. For example, January 1st is 0.0, January 2nd is 1.0, etc.
 
     Args:
-        datetime_object (`datetime.datetime`): datetime object to convert
+        - datetime_object (``datetime.datetime``): datetime object to convert
 
     Returns:
-        `float`: fractional day
+        - ``float``: fractional day
     """
     d = datetime_object-datetime.datetime(datetime_object.year-1, 12, 31)
     fractional_day = d.days + d.seconds/60./60./24 + d.microseconds/60./60./24./1e6
@@ -335,10 +363,10 @@ def from_datetime_to_mjd(datetime_obj):
     Converts a datetime object to a Modified Julian Date. The Modified Julian Date is the number of days since midnight on November 17, 1858.
 
     Args:
-        datetime_obj (`datetime.datetime`): datetime object to convert
+        - datetime_obj (``datetime.datetime``): datetime object to convert
 
     Returns:
-        `float`: Modified Julian Date
+        - ``float``: Modified Julian Date
     """
     return from_datetime_to_jd(datetime_obj)-2400000.5
 
@@ -347,10 +375,10 @@ def from_datetime_to_jd(datetime_obj):
     Converts a datetime object to a Julian Date. The Julian Date is the number of days since noon on January 1, 4713 BC.
 
     Args:
-        datetime_obj (`datetime.datetime`): datetime object to convert
+        - datetime_obj (``datetime.datetime``): datetime object to convert
 
     Returns:
-        `float`: Julian Date
+        - ``float``: Julian Date
     """
     return sum(jday(year=datetime_obj.year, mon=datetime_obj.month, day=datetime_obj.day, hr=datetime_obj.hour, minute=datetime_obj.minute, sec=datetime_obj.second+float('0.'+str(datetime_obj.microsecond))))
 
@@ -359,11 +387,11 @@ def from_cartesian_to_tle_elements(state, gravity_constant_name='wgs-72'):
     This function converts the provided state from Cartesian to TLE elements.
 
     Args:
-        state (`np.ndarray`): state to convert
-        gravity_constant_name (`str`): name of the central body (default: 'wgs-72')
-
+        - state (``np.array``): state in Cartesian coordinates
+        - gravity_constant_name (``str``): name of the gravity constant to be used (default: 'wgs-72')
+    
     Returns:
-        tuple: tuple containing: - `float`: semi-major axis - `float`: eccentricity - `float`: inclination - `float`: right ascension of the ascending node - `float`: argument of perigee - `float`: mean anomaly
+        - ``dict``: dictionary of TLE elements
     """
     _,mu_earth,_,_,_,_,_,_=get_gravity_constants(gravity_constant_name)
     mu_earth=float(mu_earth)*1e9
@@ -384,12 +412,12 @@ def from_cartesian_to_keplerian(state, mu):
     parameter of the central body, and returns the state in Keplerian elements.
 
     Args:
-        state (`np.array`): numpy array of 2 rows and 3 columns, where
+        - state (``np.array``): numpy array of 2 rows and 3 columns, where
                                     the first row represents position, and the second velocity.
-        mu (`float`): gravitational parameter of the central body
+        - mu (``float``): gravitational parameter of the central body
 
     Returns:
-        `np.array`: numpy array of the six keplerian elements: (a,e,i,omega,Omega,mean_anomaly)
+        - ``np.array``: numpy array of the six keplerian elements: (a,e,i,omega,Omega,mean_anomaly)
                                              (i.e., semi major axis, eccentricity, inclination,
                                              right ascension of ascending node, argument of perigee,
                                              mean anomaly). All the angles are in radiants, eccentricity is unitless
@@ -432,12 +460,12 @@ def from_cartesian_to_keplerian_torch(state, mu):
     Same as from_cartesian_to_keplerian, but for torch tensors.
 
     Args:
-        state (`np.array`): numpy array of 2 rows and 3 columns, where
+        - state (``torch.tensor``): torch tensor of 2 rows and 3 columns, where
                                     the first row represents position, and the second velocity.
-        mu (`float`): gravitational parameter of the central body
+        - mu (``float``): gravitational parameter of the central body
 
     Returns:
-        `np.array`: numpy array of the six keplerian elements: (a,e,i,omega,Omega,mean_anomaly)
+        - ``torch.tensor``: torch tensor of the six keplerian elements: (a,e,i,omega,Omega,mean_anomaly)
                                              (i.e., semi major axis, eccentricity, inclination,
                                              right ascension of ascending node, argument of perigee,
                                              mean anomaly). All the angles are in radiants, eccentricity is unitless