diff --git a/solventmapcreator/polynomialanalysis/polynomialcomparison.py b/solventmapcreator/polynomialanalysis/polynomialcomparison.py
new file mode 100644
index 0000000000000000000000000000000000000000..9d2e8133c2144c93746b4810bd96a7aeaba33179
--- /dev/null
+++ b/solventmapcreator/polynomialanalysis/polynomialcomparison.py
@@ -0,0 +1,107 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for the comparison of the polynomial fits of different solvents.
+
+@author: mark
+"""
+
+import logging
+import numpy as np
+import numpy.polynomial.polynomial as poly
+import scipy.spatial.distance as sciptdist
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+def calculate_hausdorff_distance_matrix(x_values, polynomial_dict_by_solvent_id):
+    """This calculates the hausdorff distance between all polynomial curves.
+    """
+    sorted_id_list = get_sorted_solvent_id_list(polynomial_dict_by_solvent_id)
+    hausdorff_matrix = np.zeros((len(sorted_id_list), len(sorted_id_list)))
+    for i in range(len(sorted_id_list)):
+        poly_values_i = polynomial_dict_by_solvent_id[sorted_id_list[i]]
+        for j in range(i):
+            poly_values_j = polynomial_dict_by_solvent_id[sorted_id_list[j]]
+            hausdorff_value = calculate_hausdorff_distance(x_values, poly_values_i, poly_values_j)
+            if i == j:
+                hausdorff_matrix[i][i] = hausdorff_value
+            else:
+                hausdorff_matrix[i][j] = hausdorff_value
+                hausdorff_matrix[j][i] = hausdorff_value
+    return hausdorff_matrix
+
+def calculate_rmsd_matrix(polynomial_dict_by_solvent_id):
+    """This calculates the rmsd between all curves, and outputs a matrix.
+    """
+    sorted_id_list = get_sorted_solvent_id_list(polynomial_dict_by_solvent_id)
+    rmsd_matrix = np.zeros((len(sorted_id_list), len(sorted_id_list)))
+    for i in range(len(sorted_id_list)):
+        poly_values_i = polynomial_dict_by_solvent_id[sorted_id_list[i]]
+        for j in range(i):
+            poly_values_j = polynomial_dict_by_solvent_id[sorted_id_list[j]]
+            rmsd_value = calculate_rmsd_between_curves(poly_values_i, poly_values_j)
+            #Noting that the matrix is symmetric
+            if i == j:
+                rmsd_matrix[i][i] = rmsd_value
+            else:
+                rmsd_matrix[i][j] = rmsd_value
+                rmsd_matrix[j][i] = rmsd_value
+    return rmsd_matrix
+
+def get_sorted_solvent_id_list(polynomial_dict_by_solvent_id):
+    """This returns a sorted list ofthe solvent IDs.
+    """
+    return sorted(polynomial_dict_by_solvent_id.keys())
+
+def calculate_hausdorff_distance(x_values, poly_values_1, poly_values_2):
+    """This calculates the Hausdorff distance between the 2 polynomial curves.
+    """
+    poly_values_1_with_x = create_array_for_hausdorff_analysis(x_values, poly_values_1)
+    poly_values_2_with_x = create_array_for_hausdorff_analysis(x_values, poly_values_2)
+    return max(sciptdist.directed_hausdorff(poly_values_1_with_x, poly_values_2_with_x)[0],
+               sciptdist.directed_hausdorff(poly_values_2_with_x, poly_values_1_with_x)[0])
+
+def calculate_rmsd_between_curves(poly_values_1, poly_values_2):
+    """This calculates the RMSD between the Assoicaiton energies along
+    the both curves. Note that both arrays must have the same dimensionality
+    for this analysis.
+    """
+    return np.sqrt(((poly_values_1 - poly_values_2)**2).mean())
+
+def create_array_for_hausdorff_analysis(x_values, polynomial_values):
+    """This creates a array of the correct dimensionality for hausdorff analysis.
+    """
+    array_size_2n = np.array([x_values, polynomial_values])
+    array_size_n2 = array_size_2n.transpose()
+    return array_size_n2
+
+def calculate_polynomial_values_all_solvents(x_values, polynomial_order,
+                                             polynomial_data_by_solvent_id):
+    """This calculates the polynomial vlaues for each solvent ID, and creates
+    a new dictionary, where the key is the solvent ID, and te value is the
+    polynomial values at the given x points.
+    """
+    poly_values_by_solvent_id = {}
+    for solvent_id, solvent_info_dict in polynomial_data_by_solvent_id.items():
+        poly_values = calculate_polynomial_values_from_dict(x_values, solvent_info_dict, polynomial_order)
+        poly_values_by_solvent_id[solvent_id] = poly_values
+    return poly_values_by_solvent_id
+
+def calculate_polynomial_values_from_dict(x_values, solvent_info_dict, polynomial_order):
+    """This calculates the polynomial values for the given info dict.
+    """
+    polynomial_coefficients = get_polynomial_coefficients(solvent_info_dict, polynomial_order)
+    return calculate_polynomial_values(x_values, polynomial_coefficients)
+
+def get_polynomial_coefficients(solvent_info_dict, polynomial_order):
+    """This returns the polynomial coefficients for the selected order.
+    """
+    return solvent_info_dict[polynomial_order]["coefficients"]
+
+
+def calculate_polynomial_values(x_values, polynomial_coefficients):
+    """This returns the values for the given polynomial function coefficients.
+    """
+    return poly.polyval(x_values, polynomial_coefficients)
diff --git a/solventmapcreator/test/polynomialanalysistest/polynomialcomparisontest.py b/solventmapcreator/test/polynomialanalysistest/polynomialcomparisontest.py
new file mode 100644
index 0000000000000000000000000000000000000000..1047d4c7c0fd421cbf154792aede42166b0ffba9
--- /dev/null
+++ b/solventmapcreator/test/polynomialanalysistest/polynomialcomparisontest.py
@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for tests of the polynomial comparison.
+
+@author: mark
+"""
+
+import logging
+import unittest
+import numpy as np
+import solventmapcreator.polynomialanalysis.polynomialcomparison as polynomialcomparison
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+class PolynomialComparisonTestCase(unittest.TestCase):
+    """Test case for the polynomialcomparison script.
+    """
+    def setUp(self):
+        """Set up before tests
+        """
+        self.poly_data_by_solvent_id = {"poly_fit":{2:{"coefficients":np.array([0.0, 
+                                                                                             0.0, 1.0]),
+                                                                "RMSE":0.0254303,
+                                                                "order":1,
+                                                                "covar":0.0008143}},
+                                        "poly_fit2":{2:{"coefficients":np.array([1.0, 
+                                                                                          0.0, 1.0]),
+                                                                 "RMSE":0.0354303,
+                                                                 "order":1,
+                                                                 "covar":0.0008143}},
+                                        "poly_fit3":{2:{"coefficients":np.array([0.0, 
+                                                                                          0.0, 2.0]),
+                                                                 "RMSE":0.0354303,
+                                                                 "order":1,
+                                                                 "covar":0.0008143}}}
+        self.x_values = np.array([-2.0, -1.0, 0.0, 1.0, 2.0])
+        self.poly_values_dict = {"poly_fit":np.array([4.0, 1.0, 0.0, 1.0, 4.0]),
+                                 "poly_fit2":np.array([5.0, 2.0, 1.0, 2.0, 5.0]), 
+                                 "poly_fit3":np.array([8.0, 2.0, 0.0, 2.0, 8.0])}
+    def tearDown(self):
+        """Tear down after tests.
+        """
+        del self.poly_data_by_solvent_id
+        del self.x_values
+        del self.poly_values_dict
+    def test_calculate_hausdorff_distance_matrix(self):
+        """
+        """
+        expected_matrix = np.array([[0.0, 1.0, 4.0],
+                                    [1.0, 0.0, 3.0],
+                                    [4.0, 3.0, 0.0]])
+        actual_matrix = polynomialcomparison.calculate_hausdorff_distance_matrix(self.x_values,
+                                                                                 self.poly_values_dict)
+        np.testing.assert_array_almost_equal(expected_matrix, actual_matrix)
+    def test_calculate_rmsd_matrix(self):
+        """Test to see if expected RMSD matrix is produced.
+        """
+        expected_matrix = np.array([[0.0, 1.0, 2.607681],
+                                    [1.0, 0.0, 1.949359],
+                                    [2.607681, 1.949359, 0.0]])
+        actual_matrix = polynomialcomparison.calculate_rmsd_matrix(self.poly_values_dict)
+        np.testing.assert_array_almost_equal(expected_matrix, actual_matrix)
+    def test_get_sorted_solvent_id_list(self):
+        """Test to see if expected list is returned.
+        """
+        expected_list = ["poly_fit", "poly_fit2", "poly_fit3"]
+        actual_list = polynomialcomparison.get_sorted_solvent_id_list(self.poly_values_dict)
+        self.assertListEqual(expected_list, actual_list)
+    def test_calculate_hausdorff_distance(self):
+        """Test to see if expected Hausdorff distance is calculated.
+        """
+        expected_hausdorff_distance = 1.0
+        actual_hausdorff_distance = polynomialcomparison.calculate_hausdorff_distance(self.x_values, np.array([4.0, 1.0, 0.0, 1.0, 4.0]),np.array([5.0, 2.0, 1.0, 2.0, 5.0]))
+        self.assertAlmostEqual(expected_hausdorff_distance, actual_hausdorff_distance)
+    def test_calculate_rmsd_between_curves(self):
+        """test to see if expected RMSD is calculated.
+        """
+        expected_rmsd = 1.0
+        actual_rmsd = polynomialcomparison.calculate_rmsd_between_curves(np.array([4.0, 1.0, 0.0, 1.0, 4.0]),np.array([5.0, 2.0, 1.0, 2.0, 5.0]))
+        self.assertAlmostEqual(expected_rmsd, actual_rmsd)
+    def test_create_array_for_hausdorff_analysis(self):
+        """Test to see if expected array is produced.
+        """
+        expected_array = np.array([[-2.0, 4.0], [-1.0, 1.0], [0.0, 0.0], [1.0, 1.0], [2.0, 4.0]])
+        actual_array = polynomialcomparison.create_array_for_hausdorff_analysis(self.x_values, np.array([4.0, 1.0, 0.0, 1.0, 4.0]))
+        np.testing.assert_array_almost_equal(expected_array, actual_array)
+    def test_calculate_polynomial_values_all_solvents(self):
+        """Test to see if expected dictionary is produced.
+        """
+        expected_dict = {"poly_fit":np.array([4.0, 1.0, 0.0, 1.0, 4.0]),
+                         "poly_fit2":np.array([5.0, 2.0, 1.0, 2.0, 5.0]), 
+                         "poly_fit3":np.array([8.0, 2.0, 0.0, 2.0, 8.0])}
+        actual_dict = polynomialcomparison.calculate_polynomial_values_all_solvents(self.x_values, 2, self.poly_data_by_solvent_id)
+        self.assertListEqual(sorted(expected_dict.keys()), sorted(actual_dict.keys()))
+        for key in expected_dict.keys():
+            np.testing.assert_array_almost_equal(expected_dict[key], actual_dict[key])
+    def test_calculate_polynomial_values_from_dict(self):
+        """Test to see if expected array is produced.
+        """
+        expected_array = np.array([4.0, 1.0, 0.0, 1.0, 4.0])
+        actual_array = polynomialcomparison.calculate_polynomial_values_from_dict(self.x_values, self.poly_data_by_solvent_id["poly_fit"], 2)
+        np.testing.assert_array_almost_equal(expected_array, actual_array)
+    def test_get_polynomial_coefficients(self):
+        """Test to see if expected coefficients are returned.
+        """
+        expected_coefficients = np.array([0.0, 0.0, 1.0])
+        actual_coefficients = polynomialcomparison.get_polynomial_coefficients(self.poly_data_by_solvent_id["poly_fit"], 2)
+        np.testing.assert_array_almost_equal(expected_coefficients, actual_coefficients)
+    def test_calculate_polynomial_values(self):
+        """Test to see if expected values are returned.
+        """
+        expected_values = np.array([4.0, 1.0, 0.0, 1.0, 4.0])
+        poly_coefficients = self.poly_data_by_solvent_id["poly_fit"][2]["coefficients"]
+        actual_values = polynomialcomparison.calculate_polynomial_values(self.x_values, poly_coefficients)
+        np.testing.assert_array_almost_equal(expected_values, actual_values)
\ No newline at end of file
diff --git a/solventmapcreator/test/solvationmapcreatortests.py b/solventmapcreator/test/solvationmapcreatortests.py
index 1cdf9d885aee53f68bb61d6c4678285f640e7590..3b64ce011c1059b050c68584f0be3318d144e9bf 100644
--- a/solventmapcreator/test/solvationmapcreatortests.py
+++ b/solventmapcreator/test/solvationmapcreatortests.py
@@ -13,6 +13,7 @@ from solventmapcreator.test.iotest.polynomialdatawritertest import PolynomialDat
 from solventmapcreator.test.iotest.polynomialdatareadertest import PolynomialDataReaderTestCase
 from solventmapcreator.test.solvationcalculationtest.solvationcalculatortest import SolvationCalculatorTestCase
 from solventmapcreator.test.polynomialanalysistest.polynomialdataanalysistest import PolynomialDataAnalysisTestCase
+from solventmapcreator.test.polynomialanalysistest.polynomialcomparisontest import PolynomialComparisonTestCase
 
 logging.basicConfig()
 LOGGER = logging.getLogger(__name__)
@@ -23,7 +24,8 @@ IO_TEST_CASES = [SolvationEnergyReaderTestCase, PolynomialDataWriterTestCase,
 
 SOLVATION_CALCULATION_TEST_CASES = [SolvationCalculatorTestCase]
 
-POLYNOMIAL_ANALYSIS_TEST_CASES = [PolynomialDataAnalysisTestCase]
+POLYNOMIAL_ANALYSIS_TEST_CASES = [PolynomialDataAnalysisTestCase,
+                                  PolynomialComparisonTestCase]
 
 def test_suite():
     """Function creates a test suite and then loads all the tests from the