diff --git a/solventmapcreator/polynomialanalysis/polynomialplotting.py b/solventmapcreator/polynomialanalysis/polynomialplotting.py
new file mode 100644
index 0000000000000000000000000000000000000000..d7fdbd14fdd4191a260f0b1e6200b8de1f1a09ec
--- /dev/null
+++ b/solventmapcreator/polynomialanalysis/polynomialplotting.py
@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for taking the datapoints for the solvation curve, and generating
+the polynomial fit data and also the plots.
+
+@author: mark
+"""
+
+import logging
+import resultsanalysis.resultsoutput.plottinginput as plottinginput
+import solventmapcreator.io.polynomialdatawriter as polynomialdatawriter
+import solventmapcreator.io.solvationenergyreader as solvationenergyreader
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+def plot_energies_from_file(free_energy_filename, polynomial_order):
+    """This creates in the plotting information, and outputs a plot to file.
+    """
+    plot_data = parse_energies_create_plot_input_data(free_energy_filename, polynomial_order)
+    return plottinginput.create_scatter_graph_with_polynomial(plot_data)
+
+def parse_energies_create_plot_input_data(free_energy_filename, polynomial_order):
+    """This creates the input data for a plot. This overrides the default figure
+    label, so that the figures can be distinguished based on the solvent.
+    """
+    datapoints = parse_free_energy_from_file_with_data_arrays(free_energy_filename)
+    plot_data = datapoints.generateScatterPlotParameters(label_type='', polynomial_order=polynomial_order)
+    plot_data['figure_label'] = free_energy_filename.replace('.xml', '')
+    return plot_data
+
+def parse_energies_and_output_poly_data(free_energy_filename, order_list, poly_filename):
+    """This reads in the energy information, and then outputs the polynomial
+    information to file.
+    """
+    datapoints = parse_free_energy_from_file_with_data_arrays(free_energy_filename)
+    return write_poly_data_to_file(datapoints, order_list, poly_filename)
+
+def write_poly_data_to_file(datapoints, order_list, filename):
+    """This calculates the polynomial information, and then outputs the
+    information to file.
+    """
+    return polynomialdatawriter.write_poly_data_to_file(datapoints, order_list, filename)
+
+def parse_free_energy_from_file_with_data_arrays(filename):
+    """This reads in the information form file to a Datapoints object, and
+    creates the data arrays.
+    """
+    return solvationenergyreader.parse_free_energy_from_file_with_data_arrays(filename)
diff --git a/solventmapcreator/test/polynomialanalysistest/polynomialplottingtest.py b/solventmapcreator/test/polynomialanalysistest/polynomialplottingtest.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a50db705360bd51d50aee4ed6bd6a9841ebe0f0
--- /dev/null
+++ b/solventmapcreator/test/polynomialanalysistest/polynomialplottingtest.py
@@ -0,0 +1,102 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for testing the polynomialplotting script.
+
+@author: mark
+"""
+
+import logging
+import os
+import numpy as np
+import unittest
+import solventmapcreator.io.solvationenergyreader as solvationenergyreader
+import solventmapcreator.polynomialanalysis.polynomialplotting as polynomialplotting
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.DEBUG)
+
+class PolynomialPlottingTestCase(unittest.TestCase):
+    """Test case for polynomialplotting script
+    """
+    def setUp(self):
+        """Set up for tests
+        """
+        self.expected_datapoints = solvationenergyreader.parse_free_energy_from_file_with_data_arrays("resources/water.xml")
+    def tearDown(self):
+        """Clean up after tests.
+        """
+        del self.expected_datapoints
+        if os.path.isfile("actual_water_poly_fit.csv"):
+            os.remove("actual_water_poly_fit.csv")
+    def test_parse_energies_create_plot_input_data(self):
+        """Test to see expected plot data is returned.
+        """
+        expected_dict = {'figure_label': "resources/water",
+                         'plot_axis_range':[-15.4     ,   7.2     , -38.635444,   0.918583],
+                         'polynomial coefficients': np.array([5.18106356e-01,
+                                                              2.54615444e-01,
+                                                              -3.76651666e-01,
+                                                              -2.06376692e-02,
+                                                              -3.76010096e-04]),
+                         'polynomial order': 4,
+                         'x_axis_range':[-15.4, 7.2000000000000002],
+                         'x_data': np.array([[-0.1, -11.1, -14.1, -15.4, -2.4,
+                                              -4.3, -5.4, -9.1, 0.5, 1.2, 7.2]]),
+                         'x_label': "ssip",
+                         'y_axis_range': [-38.635444,   0.918583],
+                         'y_data': np.array([[0.9185830776018105, -26.013780093398356,
+                                              -34.8195308084233, -38.635444108544405,
+                                              -1.4082598726557432, -6.170784422419204,
+                                              -9.315464914185585, -20.143787975013183,
+                                              0.6588682955012527, -0.4429837300591526,
+                                              -25.844114335864262]]),
+                         'y_label': "/$kJmol^{-1}$"}
+        actual_dict = polynomialplotting.parse_energies_create_plot_input_data("resources/water.xml", 4)
+        LOGGER.debug(actual_dict["x_data"].tolist())
+        LOGGER.debug(actual_dict["y_data"].tolist())
+        self.assertListEqual(sorted(actual_dict.keys()), sorted(expected_dict.keys()))
+        for key in actual_dict.keys():
+            LOGGER.debug("Key: %s", key)
+            if key == 'y_data' or key == 'x_data' or key == 'polynomial coefficients':
+                LOGGER.debug("assert equal array x and y data")
+                np.testing.assert_array_almost_equal(actual_dict[key],
+                                                     expected_dict[key])
+            elif key == 'plot_axis_range' or key == 'x_axis_range' or key == 'y_axis_range':
+                np.testing.assert_array_almost_equal(np.array(actual_dict[key]),
+                                                     np.array(expected_dict[key]))
+            else:
+                LOGGER.debug("assert equal string")
+                self.assertEqual(actual_dict[key], expected_dict[key])
+    def test_parse_energies_and_output_poly_data(self):
+        """Test to see if expected polynomial fit information is outputted to file.
+        """
+        expected_file_name = "resources/water_poly_fit.csv"
+        actual_file_name = "actual_water_poly_fit.csv"
+        poly_file_out = polynomialplotting.parse_energies_and_output_poly_data("resources/water.xml", [2, 4], actual_file_name)
+        self.assertEqual(0, poly_file_out)
+        with open(expected_file_name, 'r') as exp_file:
+            exp_file_lines = exp_file.readlines()
+            with open(actual_file_name, 'r') as act_file:
+                act_file_lines = act_file.readlines()
+                self.assertListEqual(act_file_lines, exp_file_lines)
+    def test_write_poly_data_to_file(self):
+        """Test to see if expected polynomial fit data is written to file.
+        """
+        expected_file_name = "resources/water_poly_fit.csv"
+        actual_file_name = "actual_water_poly_fit.csv"
+        poly_file_out = polynomialplotting.write_poly_data_to_file(self.expected_datapoints,
+                                                                   [2, 4], actual_file_name)
+        self.assertEqual(0, poly_file_out)
+        with open(expected_file_name, 'r') as exp_file:
+            exp_file_lines = exp_file.readlines()
+            with open(actual_file_name, 'r') as act_file:
+                act_file_lines = act_file.readlines()
+                self.assertListEqual(act_file_lines, exp_file_lines)
+    def test_parse_free_energy_from_file_with_data_arrays(self):
+        """Test to see if expected datapoints are read in.
+        """
+        actual_datapoints = polynomialplotting.parse_free_energy_from_file_with_data_arrays("resources/water.xml")
+        self.assertEqual(self.expected_datapoints, actual_datapoints)
+        
diff --git a/solventmapcreator/test/resources/water_poly_fit.csv b/solventmapcreator/test/resources/water_poly_fit.csv
new file mode 100644
index 0000000000000000000000000000000000000000..0342fd8575a940e3f3ea75fd0cc8d82a911a87fa
--- /dev/null
+++ b/solventmapcreator/test/resources/water_poly_fit.csv
@@ -0,0 +1,3 @@
+Order	RMSE	Covariance	Coefficients
+2	4.0483148	180.2773802	-3.5005663	-0.7862852	-0.2194883
+4	0.3415296	1.2830673	0.5181064	0.2546154	-0.3766517	-0.0206377	-0.0003760
diff --git a/solventmapcreator/test/solvationmapcreatortests.py b/solventmapcreator/test/solvationmapcreatortests.py
index 2655f1a3b59b3352294ebf01f5ff998ba247f49c..cedc2250c65320f9279b09207e66114399cf110b 100644
--- a/solventmapcreator/test/solvationmapcreatortests.py
+++ b/solventmapcreator/test/solvationmapcreatortests.py
@@ -18,7 +18,7 @@ from solventmapcreator.test.solvationcalculationtest.solvationplotinformationtes
 from solventmapcreator.test.solvationcalculationtest.fractionaloccupancycalculatortest import FractionalOccupancyCalculatorTestCase
 from solventmapcreator.test.polynomialanalysistest.polynomialdataanalysistest import PolynomialDataAnalysisTestCase
 from solventmapcreator.test.polynomialanalysistest.polynomialcomparisontest import PolynomialComparisonTestCase
-
+from solventmapcreator.test.polynomialanalysistest.polynomialplottingtest import PolynomialPlottingTestCase
 logging.basicConfig()
 LOGGER = logging.getLogger(__name__)
 LOGGER.setLevel(logging.WARN)
@@ -32,7 +32,8 @@ SOLVATION_CALCULATION_TEST_CASES = [SolvationCalculatorTestCase,
                                     FractionalOccupancyCalculatorTestCase]
 
 POLYNOMIAL_ANALYSIS_TEST_CASES = [PolynomialDataAnalysisTestCase,
-                                  PolynomialComparisonTestCase]
+                                  PolynomialComparisonTestCase,
+                                  PolynomialPlottingTestCase]
 
 def test_suite():
     """Function creates a test suite and then loads all the tests from the