--- a/app/nenzfr/nenzfr_stats.py	Tue May 05 11:13:34 2020 +1000
+++ b/app/nenzfr/nenzfr_stats.py	Wed Jun 03 20:53:51 2020 +1000
@@ -8,6 +8,11 @@
     Area and mass weighted statistics are also calculated and written out.
 22-Aug-2012
     Brute-force approach to estimating the CMME 1D flow conditions.
+04-May-2020
+    Function add_extra_variables (which is useful for uncertainty analysis) 
+    is now included here. This is in preparation for the removal of
+    nenzfr_perturbed.py, nenzfr_sensitivity.py. These have been superceded
+    with perturb.py and sensitivity.py
 """
 
 import sys, os
@@ -66,6 +71,32 @@
     return variable_list, data
 
 
+def add_extra_variables(data, var):
+    """
+    We add some additional variables that of are interest to
+    the input data dictionary and list. 
+
+    data: dictionary of data
+    var: list of variables
+    """
+    U_squared = data['vel.x']*data['vel.x'] + data['vel.y']*data['vel.y']
+    U = U_squared**0.5
+    # Dynamic Pressure
+    data['q'] = 0.5 * data['rho'] * U_squared
+    var.append('q')
+    # Mass flow rate per unit area
+    data['m_dot'] = data['rho'] * U
+    var.append('m_dot')
+    # Unit Reynolds number
+    data['Re_u'] = data['m_dot'] / data['mu']
+    var.append('Re_u')
+    # Pressure coefficient
+    data['p/q'] = data['p'] / data['q']
+    var.append('p/q')
+
+    return data, var
+
+
 def print_stats_MoA(sliceFileName,jobName,coreRfraction,weight):
     """
     (MoA = Mass- or Area-weighted)

--- a/app/nenzfr/nenzfr_utils.py	Tue May 05 11:13:34 2020 +1000
+++ b/app/nenzfr/nenzfr_utils.py	Wed Jun 03 20:53:51 2020 +1000
@@ -3,9 +3,11 @@
 """
 
 import sys, os
-import shlex, subprocess, string
+import shlex, subprocess, string, copy
 E3BIN = os.path.expandvars("$HOME/e3bin")
 sys.path.append(E3BIN)
+from nenzfr_stats import add_extra_variables
+
 
 #---------------------------------------------------------------
 
@@ -67,32 +69,33 @@
 #   nenzfr_sensitivity.py
 #   nenzfr_RSA.py
 
-def read_case_summary(FileToRead=None):
-    """
-    Reads the file "perturbation_cases.dat" to determine which variables
-    have been perturbed and values for the perturbed variables for each
-    case consituting the sensitivity calculation.
-
-    :returns: perturbedVariables - a list of variable names
-            : DictOfCases - a dictionary with the case names as keys
-                            and the values of each perturbed variable.
-    """
-    if FileToRead is None:
-        fp = open('perturbation_cases.dat','r')
-    else:
-        fp = open(FileToRead,'r')
-
-    varList = fp.readline().strip().split(" ")
-    perturbedVariables = [k for k in varList if k!="#" and k!=""]
-    fp.readline()
-    DictOfCases = {}
-    for line in fp.readlines():
-        caseData = line.strip().split(" ")
-        caseName = caseData[0]
-        DictOfCases[caseName] = [float(k) for k in caseData \
-                                 if k!=caseName and k!=""]
-    fp.close()
-    return perturbedVariables, DictOfCases
+# 04-May-2020: read_case_summary now included in perturb_utils.py
+#def read_case_summary(FileToRead=None):
+#    """
+#    Reads the file "perturbation_cases.dat" to determine which variables
+#    have been perturbed and values for the perturbed variables for each
+#    case consituting the sensitivity calculation.
+#
+#    :returns: perturbedVariables - a list of variable names
+#            : DictOfCases - a dictionary with the case names as keys
+#                            and the values of each perturbed variable.
+#    """
+#    if FileToRead is None:
+#        fp = open('perturbation_cases.dat','r')
+#    else:
+#        fp = open(FileToRead,'r')
+#
+#    varList = fp.readline().strip().split(" ")
+#    perturbedVariables = [k for k in varList if k!="#" and k!=""]
+#    fp.readline()
+#    DictOfCases = {}
+#    for line in fp.readlines():
+#        caseData = line.strip().split(" ")
+#        caseName = caseData[0]
+#        DictOfCases[caseName] = [float(k) for k in caseData \
+#                                 if k!=caseName and k!=""]
+#    fp.close()
+#    return perturbedVariables, DictOfCases
 
 def read_nenzfr_outfile(FileToRead,inclStats=0):
     """
@@ -164,6 +167,47 @@
     fp.close()
     return supplyDict
 
+
+def read_output_for_sensitivity_analysis(FileToRead):
+    """
+    A simple function that reads both nenzfr outfile and estcj
+    outfile and combines the results, returning a dictionary of 
+    values and list of variables.
+
+    It is intended that this function be imported into 'read_outfile.py',
+    which must be created by the user in order to run a sensitivity
+    analysis (using perturb.py and sensitivity.py
+    
+    It reads files: 
+        FileToRead-exit.stats
+        FileToRead-estcj.dat
+    """
+    dataDict = {}
+    varList  = {}
+    
+    exitData, exitVar = \
+              read_nenzfr_outfile(FileToRead+'-exit.stats')
+    supplyData = read_estcj_outfile(FileToRead+'-estcj.dat')
+
+    dataDict = copy.deepcopy(exitData)
+    varList  = copy.deepcopy(exitVar)
+    
+    # Add the relevant nozzle supply data
+    dataDict['supply_rho'] = supplyData['rho']
+    dataDict['supply_T']   = supplyData['T']
+    dataDict['supply_h']   = supplyData['h']
+
+    # Add the variables to the list
+    varList.insert(0,'supply_rho')
+    varList.insert(1,'supply_T')
+    varList.insert(2,'supply_h')
+
+    # Add extra freestream properties of interest
+    dataDict, varList = add_extra_variables(dataDict, varList)
+
+    return dataDict, varList
+
+
 #---------------------------------------------------------------
 # The following is required by:
 #    nenzfr_perturbed.py
@@ -204,6 +248,7 @@
     
     return scriptFileName, cfgFileName
 
+
 def build_cfg_file(cfg, filename):
     """
     Takes a nenzfr config dictionary and uses it to build a nenzfr config file.close

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/app/nenzfr/perturb.py	Wed Jun 03 20:53:51 2020 +1000
@@ -0,0 +1,175 @@
+#!/usr/bin/env python
+""" perturb.py
+
+This script coordinates the running of a series of calculations 
+that are perturbations around a nominal condition.
+
+.. Author: Luke Doherty (luke.doherty@eng.ox.ac.uk)
+           Oxford Thermofluids Institute
+           The University of Oxford
+"""
+
+VERSION_STRING = "04-May-2020"
+
+import shlex, subprocess, string
+from subprocess import PIPE
+import sys, os, gzip
+import optparse
+from numpy import array, mean, logical_and, zeros
+from perturb_input_utils import perturb_input_checker
+from perturb_utils import run_command, \
+                          set_case_running, set_perturbed_values,\
+                          write_case_config, write_case_summary
+import copy
+E3BIN = os.path.expandvars("$HOME/e3bin")
+sys.path.append(E3BIN)
+
+#---------------------------------------------------------------
+def run_perturb(cfg):
+    """
+    Function that accepts a config dictionary and runs perturb.
+    """
+    #
+    cfg['bad_input'] = False
+    #
+    # First check our input and assign any defaults values
+    cfg = perturb_input_checker(cfg)
+    #bail out here if there is an issue
+    if cfg['bad_input']:
+        return -2
+    #
+    # Unpack the input dictionary of nominal values into the main 
+    # configuration dictionary
+    for var in cfg['nominalValues'].keys():
+        cfg[var] = cfg['nominalValues'][var]
+    #
+    # Re-assign the dictionary of perturbed values (for later convenience)
+    perturbedDict = cfg['perturbedDict']
+    #
+    if not cfg['createRSA']: # Perturbing for a sensitivity calculation
+        # Calculate Nominal condition
+        caseString = 'case'+"{0:02}".format(0)+"{0:01}".format(0)
+        textString = "Nominal Condition"
+        caseDict = copy.copy(cfg)
+        caseDict['caseName'] = caseString
+        write_case_config(caseDict)
+        # Run the nominal case and write the values of the perturbed 
+        # variables to a summary file
+        set_case_running(caseString, caseDict, textString,'')
+        write_case_summary(cfg['variablesToPerturb'],caseDict,caseString,1)
+        #
+        # Now run all the perturbed conditions
+        for k in range(len(cfg['variablesToPerturb'])):
+            var = cfg['variablesToPerturb'][k]
+            perturbCount = cfg['levels']
+
+            for kk in range(perturbCount):
+                if kk != 0:
+                    caseString = 'case'+"{0:02}".format(k)+\
+                                 "{0:01}".format(kk)
+                    textString = var+" perturbed to "+\
+                                 str(perturbedDict[var][kk])
+                    caseDict = copy.copy(cfg)
+                    caseDict[var] = perturbedDict[var][kk]
+                    caseDict['caseName'] = caseString
+                    # Run the current case
+                    set_case_running(caseString, caseDict, \
+                                     textString, var)
+                    # Write current case to the summary file
+                    write_case_summary(cfg['variablesToPerturb'],\
+                                       caseDict,caseString,0)
+    #
+    else: # Perturbing to create a Response Surface
+        # Currently we can only fit a response surface through
+        # two variables, so we take the first two in the input list
+        var1 = cfg['variablesToPerturb'][0] # 'Vs'
+        var2 = cfg['variablesToPerturb'][1] # 'pe'
+
+        if cfg['levels'] == 2.5:
+            casesToRun = [(2,1),      (1,1),
+                                (0,0),
+                          (2,2),      (1,2)]
+        elif cfg['levels'] == 3:
+            casesToRun = [(2,1),(0,1),(1,1),
+                          (2,0),(0,0),(1,0),
+                          (2,2),(0,2),(1,2)]
+        elif cfg['levels'] == 5:
+            casesToRun = [(4,3),      (0,3),      (3,3),
+                                (2,1),      (1,1),
+                          (4,0),      (0,0),      (3,0),
+                                (2,2),      (1,2),
+                          (4,4),      (0,4),      (3,4)]
+            #casesToRun = [      (2,3),      (1,3),
+            #              (4,1),      (0,1),      (3,1),
+            #                    (2,0),      (1,0),
+            #              (4,2),      (0,2),      (3,2),
+            #                    (2,4),      (1,4)      ]
+        #
+        # Run the nominal case first
+        caseString = 'case'+"{0:01}{0:01}".format(0,0)
+        caseDict = copy.copy(paramDict)
+        caseDict['caseName'] = caseString
+        textString = "Nominal Case: "+var1+"="+str(perturbedDict[var1][0])+\
+                                 "; "+var2+"="+str(perturbedDict[var2][0])
+        write_case_config(caseDict)
+        set_case_running(caseString, caseDict, textString, '')
+        write_case_summary(cfg['variablesToPerturb'], caseDict, caseString, 1)
+        #
+        # Now run all other cases
+        for case in casesToRun:
+            if case != (0,0):
+                caseString = 'case'+"{0:01}{1:01}".format(case[0],case[1])
+                textString = var1+" perturbed to "+\
+                             str(perturbedDict[var1][case[0]])+\
+                        "\n"+var2+" perturbed to "+\
+                             str(perturbedDict[var2][case[1]])
+                caseDict = copy.copy(paramDict)
+                caseDict['caseName'] = caseString
+                caseDict[var1] = perturbedDict[var1][case[0]]
+                caseDict[var2] = perturbedDict[var2][case[1]]
+                #
+                set_case_running(caseString, caseDict, textString, '')
+                write_case_summary(cfg['variablesToPerturb'], \
+                                   caseDict, caseString, 0)
+    #
+    return
+
+
+def main():
+    """
+    Examine the command-line options to decide the what to do
+    and then coordinate a series of calculations using inputs
+    that are perturbed around specified nominal values.
+    """
+    op = optparse.OptionParser(version=VERSION_STRING)
+    op.add_option('-c', '--config_file', dest='config_file',
+                  help=("filename for the config file"))
+    opt, args = op.parse_args()
+    config_file = opt.config_file
+    #
+    cfg = {}
+    #   
+    if not cfg: #if the configuration dictionary has not been filled up already, load it from a file
+        try: #from Rowan's onedval program
+            execfile(config_file, globals(), cfg)
+        except IOError as e:
+            print "Error {0}".format(str(e))
+            print "There was a problem reading the config file: '{0}'".format(config_file)
+            print "Check that it conforms to Python syntax."
+            print "Bailing out!"
+            sys.exit(1)
+    #
+    run_perturb(cfg)
+    #
+    return
+
+#---------------------------------------------------------------
+
+if __name__ == '__main__':
+    if len(sys.argv) <= 1:
+        print "perturb:\n run simulations that are perturbations about a nominal"
+        print "   Version:", VERSION_STRING
+        print "   To get some useful hints, invoke the program with option --help."
+        sys.exit(0)
+    return_flag = main()
+    sys.exit(return_flag)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/app/nenzfr/perturb_input_utils.py	Wed Jun 03 20:53:51 2020 +1000
@@ -0,0 +1,139 @@
+"""
+perturb_input_utils.py -- Function to check the perturb input dictionary
+    that is pulled in from a config file. Default values are added 
+    where appropriate.
+
+.. Author: Luke Doherty (luke.doherty@eng.ox.ac.uk)
+           Oxford Thermofluids Institute
+           The University of Oxford
+
+  Version: 04-May-2020
+"""
+
+from perturb_utils import set_perturbed_values
+
+def perturb_input_checker(cfg):
+    """Takes the config dictionary and checks that all the necessary
+    inputs are in there. Will bail out if important variables are missing, 
+    and set default values for those that are not crucial.
+       
+    Returns the config dictionary after it has been checked over and will
+    tell perturb to bail out if it finds an issue.
+    
+    """
+
+    print "Checking perturb inputs."
+    #
+    if 'runCMD' not in cfg:
+        cfg['runCMD'] = './'
+        print "Run command not specified. We will assume you are running"
+        print "    the simulations on a local machine that does not have"
+        print "    a que. Setting to a default run command ('{0}')".format(cfg['runCMD'])
+    #
+    if 'levels' not in cfg:
+        cfg['levels'] = 3
+        print "The number of perturbation levels to use is not specified."
+        print "    Setting it to default value of {0}".format(cfg['levels'])
+    #
+    if 'createRSA' not in cfg:
+        print "Response surface approximation variable not set."
+        print "    Setting it to False."
+    #
+    if cfg['levels'] == '3-reduced':
+        if cfg['createRSA']:
+            cfg['levels'] = 2.5
+        else:
+            print "Setting levels to '3-reduced' is only valid when creating a"
+            print "    response surface. Changing this to levels = 3."
+            cfg['levels'] = 3
+    else:
+        cfg['levels'] = int(cfg['levels'])
+    #
+    if 'variablesToPerturb' not in cfg:
+        print "You have not specified a list of which variables are to be perturbed."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+    #
+    if 'nominalValues' not in cfg:
+        print "You have not specified a dictionary of default values for each variable"
+        print "    that will be perturbed."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+    else:
+        for var in cfg['variablesToPerturb']:
+            if var not in cfg['nominalValues']:
+                print "{0} is not present in the nominalValues dictionary."
+                print "    Bailing out."
+                cfg['bad_input'] = True
+    # 
+    if 'typeOfPerturbation' not in cfg and 'perturbationMagnitudes' not in cfg:
+        print "Neither the typeOfPerturbation or perturbationMagnitudes dictionaries"
+        print "    have been specified. Setting a default 'relative' perturbation of"
+        print "    1.0 percent for each variable."
+        cfg['typeOfPerturbation'] = {}
+        cfg['perturbationMagnitudes'] = {}
+        for var in cfg['variablesToPerturb']:
+            cfg['typeOfPerturbation'][var] = 'relative'
+            cfg['perturbationMagnitudes'][var] = [1.0] # This must a list...
+    elif 'typeOfPerturbation' not in cfg and 'perturbationMagnitudes' in cfg:
+        print "You must specify either both or neither of the typeOfPerturbation"
+        print "    and perturbationMagnitudes dictionaries."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+    elif 'typeOfPerturbation' in cfg and 'perturbationMagnitudes' not in cfg:
+        print "You must specify either both or neither of the typeOfPerturbation"
+        print "    and perturbationMagnitudes dictionaries."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+    else:
+        for var in cfg['variablesToPerturb']:
+            if var not in cfg['typeOfPerturbation']:
+                print "{0} is not present in typeOfPerturbation dictionary.".format(var)
+                print "    Bailing out."
+                cfg['bad_input'] = True
+            if var not in cfg['perturbationMagnitudes']:
+                print "{0} is not present in perturbationMagnitudes dictionary.".format(var)
+                print "    Bailing out."
+                cfg['bad_input'] = True
+    # 
+    #if 'inputFileTemplate' not in cfg:
+    #    print "You have not specified the inputFileTemplate file name."
+    #    print "    Bailing out."
+    #    cfg['bad_input'] = True
+    #
+    if 'runScriptTemplate' not in cfg:
+        print "You have not specified the runScriptTemplate file name."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+    #
+    # Warn the user that only the first two elements of the input 
+    # variablesToPerturb list will be perturbed when we want to 
+    # create a response surface
+    if cfg['createRSA']:
+        print "Create Response Surface (RS) option selected."
+        if len(cfg['variablesToPerturb'])>2:
+            print "NOTE: We can currently only fit a response surface to"
+            print "two variables."
+
+    if 'extraFilesToCopy' not in cfg:
+        print "No extra files to copy have specified."
+        cfg['extraFilesToCopy'] = []
+
+    # Now go through and calculate the perturbed values for each variable
+    # of interest.
+    if not cfg['bad_input']:
+        cfg['perturbedDict'] = {} #this will get populated below
+        #
+        for variable in cfg['variablesToPerturb']: 
+            cfg['perturbedDict'][variable] = set_perturbed_values( \
+                                                 cfg['nominalValues'][variable], \
+                                                 cfg['perturbationMagnitudes'][variable],\
+                                                 cfg['typeOfPerturbation'][variable], \
+                                                 cfg['levels'])
+
+            print "Done with perturbed variable {0}.".format(variable)
+
+    print "Done checking perturb inputs."
+
+    return cfg
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/app/nenzfr/perturb_utils.py	Wed Jun 03 20:53:51 2020 +1000
@@ -0,0 +1,282 @@
+"""
+perturb_utils.py -- Small untility functions needed by the
+                    perturb and sensitivity  programs.
+
+.. Author: Luke Doherty (luke.doherty@eng.ox.ac.uk)
+           Oxford Thermofluids Institute
+           The University of Oxford
+
+  Version: 04-May-2020
+"""
+
+import sys, os
+import shlex, subprocess, string
+E3BIN = os.path.expandvars("$HOME/e3bin")
+sys.path.append(E3BIN)
+
+#-------------------------------------------------------------------------
+def run_command(cmdText):
+    """
+    Run the command as a subprocess. Copied verbatim from nenzfr_utils.py
+    """
+    
+    # Flush before using subprocess to ensure output is in the right order.
+    sys.stdout.flush()
+
+    if (type(cmdText) is list):
+        args = cmdText
+    else:
+        args = shlex.split(cmdText)
+    print "About to run cmd:", string.join(args)
+    # p = subprocess.Popen(args)
+    # wait until the subprocess is finished
+    # stdoutData, stderrData = p.communicate()
+    return subprocess.check_call(args)
+
+
+def set_case_running(caseString, caseDict, textString, var):
+    """
+    A short function to set a given case running in an appropriately
+    named sub-directory.
+    """
+    print 60*"-"
+    print caseString
+    print textString
+    #
+    if var not in ['coreRfraction','CoreRadiusFraction']:
+        # Create sub-directory for the current case
+        run_command('mkdir ./'+caseString)
+        #
+        # Set up the run script
+        scriptFileName = caseDict['runScriptTemplate']
+        create_file_from_template(scriptFileName, './'+caseString+\
+                                  '/'+scriptFileName, caseDict)
+        #
+        # Set up the input file
+        if 'inputFileTemplate' in caseDict:
+            inputFileName = caseDict['inputFileTemplate']
+            create_file_from_template(inputFileName,'./'+caseString+\
+                                      '/'+inputFileName, caseDict)
+        #
+        # If required, copy additional files to the sub-directory
+        for file in caseDict['extraFilesToCopy']:
+            command_text = 'cp ./'+file+' ./'+caseString+'/'
+            run_command(command_text)
+        #
+        # Change into the sub-directory, ensure the run script and 
+        # input file are exectuable, then execute the run script
+        os.chdir(caseString)
+        run_command('chmod u+x '+scriptFileName)
+        if 'inputFileTemplate' in caseDict:
+            run_command('chmod u+x '+inputFileName)
+        print ""
+        print caseDict['runCMD']+scriptFileName
+        print ""
+        # I am not sure how to replace the following line with the 
+        # run_command function
+        os.system(caseDict['runCMD']+scriptFileName)
+        os.chdir('../')
+    else:
+        # Below is very specific to nenzfr. I don't know how best to remove
+        # this so I've left it. LukeD 05-05-2020
+        if not os.path.exists(caseString):
+            # Copy nominal case data into a new case directory
+            command_text = 'cp -r case000 '+caseString
+            run_command(command_text)
+            #
+            # Change into new directory and call nenzfr --just-stats
+            os.chdir(caseString)
+            command_text = \
+                'nenzfr.py --just-stats --CoreRadiusFraction='+\
+                str(caseDict[var])
+                #str(DictOfCases[caseString][perturbedVariables.index(var)])
+            print ""
+            print command_text
+            run_command(command_text)
+            print ""
+            os.chdir('../')
+    #
+    return
+
+def create_file_from_template(templateFileName, outPutFileName, substituteDict):
+    """Function that will write out a new file by undertaking a dictionary 
+    substitution on a nominated template file.
+    """
+    #
+    fp = open(templateFileName, 'r')
+    text = fp.read()
+    fp.close()
+    template = string.Template(text)
+    text = template.safe_substitute(substituteDict)
+    fp = open(outPutFileName, 'w')
+    fp.write(text)
+    fp.close()
+    #
+    return
+
+def set_perturbed_values(NominalValue, PerturbationMagnitude, TypeOfPerturbation, Levels):
+    """
+    
+    """
+    if TypeOfPerturbation == "relative":
+        delta = [k/100.0*NominalValue for k in PerturbationMagnitude]
+    else:
+        delta = PerturbationMagnitude
+    #
+    if Levels == 3:
+        if len(delta) == 1:
+            # have [delta]
+            values = [NominalValue, NominalValue+delta[0],\
+                                    NominalValue-delta[0]]
+        elif len(delta) == 2:
+            # have [+delta, -delta]
+            values = [NominalValue, NominalValue+delta[0],\
+                                    NominalValue-delta[1]]
+    elif Levels == 5:
+        if len(delta) == 1:
+            # have [delta]
+            values = [NominalValue, NominalValue+delta[0],\
+                                    NominalValue-delta[0],\
+                                    NominalValue+2.0*delta[0],\
+                                    NominalValue-2.0*delta[0]]
+        elif len(delta) == 2:
+            # have [delta, "2"*delta]
+            values = [NominalValue, NominalValue+delta[0],\
+                                    NominalValue-delta[0],\
+                                    NominalValue+delta[1],\
+                                    NominalValue-delta[1]]
+        elif len(delta) == 3:
+            # have [+delta, -delta, "2"*delta]
+            values = [NominalValue, NominalValue+delta[0],\
+                                    NominalValue-delta[1],\
+                                    NominalValue+delta[2],\
+                                    NominalValue-delta[2]]
+        elif len(data) == 5:
+            # have [+delta, -delta, +"2"*delta, -"2"*delta]
+            values = [NominalValue, NominalValue+delta[0],\
+                                    NominalValue-delta[1],\
+                                    NominalValue+delta[2],\
+                                    NominalValue-delta[3]]
+        
+    return values
+
+def write_case_summary(varList,caseDict,caseString,newfile):
+    """
+    A short function to write the values of the perturbed variables 
+    for the current case to a summary file.
+    """
+    # Define some format strings
+    #l = max( max([len(k) for k in varList]), 15 ) 
+    #dataFormat = '{0:>'+str(l)+'.6g}'
+    #titleFormat = '{0:{fill}>'+str(l)+'}'
+    dataFormat = {k : '{0:>'+"{0}".format(max(11,len(k))+2)+".5g}" \
+                      for k in varList}
+    titleFormat = {k : '{0:{fill}>'+"{0}".format(max(11,len(k))+2)+"}" \
+                       for k in varList}
+     
+    # For the first time, we create a new file and write the header information.
+    # Each following case is appended to the existing file.
+    if newfile == 1:
+        fout = open('perturbation_cases.dat','w')
+        # Write title line
+        fout.write('{0:>7}'.format('#'))
+        for k in tuple(varList):
+            fout.write(titleFormat[k].format(k,fill=''))
+        fout.write('\n')
+        # Underline the title
+        for k in varList:
+            fout.write(titleFormat[k].format('-',fill='-'))
+        fout.write('{0:->7}'.format('-'))
+        fout.write('\n')
+    else:
+        fout = open('perturbation_cases.dat','a')
+    # Now write out the data for the current case
+    fout.write('{0:>7}'.format(caseString))
+    for k in varList:
+        fout.write(dataFormat[k].format(caseDict[k]))
+    fout.write('\n')
+    fout.close()
+
+def read_case_summary(FileToRead=None):
+    """
+    Reads the file "perturbation_cases.dat" to determine which variables
+    have been perturbed and values for the perturbed variables for each
+    case consituting the sensitivity calculation.
+
+    :returns: perturbedVariables - a list of variable names
+            : DictOfCases - a dictionary with the case names as keys
+                            and the values of each perturbed variable.
+    """
+    if FileToRead is None:
+        fp = open('perturbation_cases.dat','r')
+    else:
+        fp = open(FileToRead,'r')
+
+    varList = fp.readline().strip().split(" ")
+    perturbedVariables = [k for k in varList if k!="#" and k!=""]
+    fp.readline()
+    DictOfCases = {}
+    for line in fp.readlines():
+        caseData = line.strip().split(" ")
+        caseName = caseData[0]
+        caseValues = [float(k) for k in caseData if k!=caseName and k!=""]
+        DictOfCases[caseName] = dict(zip(perturbedVariables,caseValues))
+        #DictOfCases[caseName] = [float(k) for k in caseData \
+        #                         if k!=caseName and k!=""]
+    fp.close()
+    return perturbedVariables, DictOfCases
+
+
+def write_case_config(caseDict):
+    """Short function that writes a file summarising the values of
+    an input dictionary. 
+    """
+    fout = open('nominal_case.config','w')
+    for k in range(len(caseDict)):
+        fout.write('{0}'.format(caseDict.keys()[k]+':   '))
+        fout.write('{0}'.format(caseDict.values()[k]))
+        fout.write('\n')
+    fout.close()
+
+
+# The below function was originally used in nenzfr_sensitivity. Parts of
+# the code have now been incorporated into the function "set_case_running"
+# (see above). LukeD 05-05-2020
+# 
+#def perturb_CoreRadiusFraction(var, perturbedVariables,\
+#                               DictOfCases, levels):
+#    """
+#    Perturbation of the CoreRadiusFraction may be completed
+#    as a post-processing step using the nominal condition
+#    flow solution. We don't need to compute separate solutions.
+#
+#    This function copies the nominal condition solution to a
+#    new case file and then runs nenzfr with the --just-stats
+#    option and the perturbed CoreRadiusFraction value.
+#    """
+#
+#    for kk in range(levels):
+#        if kk != 0:
+#            caseString = 'case'+"{0:02}".format(\
+#                 perturbedVariables.index(var))+\
+#                 "{0:01}".format(kk)
+#            # Only if required do we proceed with copying the
+#            # nominal case data and calling 'nenzfr.py'
+#            if not os.path.exists(caseString):
+#                print "Perturbing "+var
+#                #command_text = 'rm -r '+caseString
+#                #run_command(command_text)
+#
+#                # Copy nominal case data into a new case
+#                # directory
+#                command_text = 'cp -r case000 '+caseString
+#                run_command(command_text)
+#                # Change into new directory and call nenzfr --just-stats
+#                os.chdir(caseString)
+#                command_text = \
+#                    'nenzfr.py --just-stats --CoreRadiusFraction='+\
+#                    str(DictOfCases[caseString][perturbedVariables.index(var)])
+#                run_command(command_text)
+#                os.chdir('../')
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/app/nenzfr/sensitivity.py	Wed Jun 03 20:53:51 2020 +1000
@@ -0,0 +1,384 @@
+#!/usr/bin/env python
+""" sensitivity.py
+
+This script calculates the sensitivity and uncertainty in each
+output property due to the various input parameters using a set
+of calculations generated by "perturb.py".
+
+The calculation method follows that outlined in the text book:
+    Coleman, HW and Steele, WG (1999). Experimentation and
+    Uncertainty Analysis for Engineers. 2nd edition. New York,
+    John Wiley & Sons, Inc. DOI: 10.2514/2.3773 
+
+Note that this script DOES NOT account for correlated uncertainties;
+that is, it assumes the input variable uncertainties are completely 
+independent of one another. 
+
+The user must supply a function called "read_outfile.py" which reads
+the outfile of the users analysis code and returns a dictionary of 
+results and list of output variables (to be analysed).
+ 
+.. Author: Luke Doherty (luke.doherty@eng.ox.ac.uk)
+           Oxford Thermofluids Institute
+           The University of Oxford
+"""
+
+VERSION_STRING = "05-May-2020"
+
+import shlex, string
+import sys, os, copy
+E3BIN = os.path.expandvars("$HOME/e3bin")
+sys.path.append(E3BIN) # installation directory
+sys.path.append("") # so that we can find user's scripts in current directory
+import optparse
+import numpy as np #import array
+np.seterr(divide='ignore',invalid='ignore')
+from perturb_utils import read_case_summary
+from sensitivity_input_utils import sensitivity_input_checker
+from sensitivity_utils import get_values, write_sensitivity_summary,\
+                              write_uncertainty_summary
+from read_outfile import read_outfile
+
+#from nenzfr_utils import run_command, quote, read_case_summary, \
+#     read_nenzfr_outfile, read_estcj_outfile
+#from nenzfr_input_utils import input_checker, nenzfr_perturbed_input_checker, nenzfr_sensitivity_input_checker
+
+#---------------------------------------------------------------
+def run_sensitivity(cfg):
+    """
+    Function that accepts a config dictionary and runs a 
+    sensitivity analysis
+    """
+    #
+    cfg['bad_input'] = False
+    #
+    # First check our input
+    cfg = sensitivity_input_checker(cfg)
+    # bail out here if there is an issue
+    if cfg['bad_input']:
+        return -2
+    
+    # Read the sensitivity_case_summary file to get the perturbed
+    # variables and their various values
+    perturbedVariables, dictOfCases = read_case_summary()
+    #dontNeed, dictOfCases = read_case_summary()
+    #perturbedVariables = copy.deepcopy(cfg['inputVariables'])
+    
+    #print perturbedVariables
+    #print dictOfCases
+    #print
+    #print cfg
+    
+    # Make a copy of list of variables we want to analyse
+    exitVar = copy.deepcopy(cfg['outputVariables'])
+
+    # Define the name of the nominal case and load the exit plane data
+    nominal = 'case000'
+    nominalData, allExitVar = \
+                        read_outfile('./'+nominal+'/'+cfg['resultsFile']) 
+    
+    #nominalData, exitVar = read_nenzfr_outfile('./'+nominal+'/'+\
+    #                                           cfg['jobName']+'-exit.stats')
+    # Load the nozzle supply data
+    #nominalSupply = read_estcj_outfile('./'+nominal+'/'+cfg['jobName']+'-estcj.dat')
+    # Now add the relevant supply data (T, h) to the nominalData dictionary
+    #nominalData['supply_rho'] = nominalSupply['rho']
+    #nominalData['supply_T'] = nominalSupply['T']
+    #nominalData['supply_h'] = nominalSupply['h']
+    # Add the supply variables to the exitVar list
+    #exitVar.insert(0,'supply_rho')
+    #exitVar.insert(1,'supply_T')
+    #exitVar.insert(2,'supply_h')
+    # Add extra variables of interest (q, Re_u, m_dot)
+    #nominalData, exitVar = add_extra_variables(nominalData, exitVar)
+     
+    #nominalValues = get_values(nominalData, exitVar)
+    #print 'nominalValues: ',nominalValues
+    
+     
+    # Extract the input uncertainties for easier use. We must use the
+    # relative values
+    inputUncertainties = copy.deepcopy(cfg['inputVariablesUncertainty_rel'])
+    
+    
+    # Loop through each of the perturbed variables
+    sensitivity_rel = {}
+    #sensitivity_rel2 = {}
+    sensitivity_abs = {}
+    #sensitivity_abs2 = {}
+    uncertainty = {}
+    #uncertainty2 = {}
+
+    for k in range(len(perturbedVariables)):
+        var = perturbedVariables[k]
+        print var
+        #if var == 'CoreRadiusFraction':
+        #    perturb_CoreRadiusFraction(var, perturbedVariables,\
+        #                               dictOfCases, cfg['levels'])
+        
+        # Define the name of the relevant perturbed cases and load the
+        # associated data
+        high = 'case'+"{0:02}".format(k)+'1'
+        highData, dontNeed = read_outfile('./'+high+'/'+cfg['resultsFile'])
+        
+        #highData, dontNeed  = read_nenzfr_outfile('./'+high+'/'+\
+        #                                          cfg['jobName']+'-exit.stats')
+        #highSupply = read_estcj_outfile('./'+high+'/'+cfg['jobName']+'-estcj.dat')
+        #highData['supply_rho'] = highSupply['rho']
+        #highData['supply_T'] = highSupply['T']
+        #highData['supply_h'] = highSupply['h']
+        #highData, dontNeed = add_extra_variables(highData, [])
+
+        low = 'case'+"{0:02}".format(k)+'2'
+        lowData, dontNeed = read_outfile('./'+low+'/'+cfg['resultsFile'])
+        
+        #lowData, dontNeed = read_nenzfr_outfile('./'+low+'/'+\
+        #                                        cfg['jobName']+'-exit.stats')
+        #lowSupply = read_estcj_outfile('./'+low+'/'+cfg['jobName']+'-estcj.dat')
+        #lowData['supply_rho'] = lowSupply['rho']
+        #lowData['supply_T'] = lowSupply['T']
+        #lowData['supply_h'] = lowSupply['h']
+        #lowData, dontNeed = add_extra_variables(lowData, [])
+
+        print low, nominal, high
+        
+        # Values of the freestream properties at the perturbed conditions
+        #highValues = get_values(highData,exitVar)
+        #lowValues = get_values(lowData,exitVar)
+        #print 'highValues: ',highValues
+        #print 'lowValues: ',lowValues
+        #print 'highData: ',highData
+        #print 'lowData: ',lowData
+        #print 'nominalData: ',nominalData
+        #print dictOfCases[high]
+        #print perturbedVariables
+        #print perturbedVariables.index(var)
+        #print dictOfCases[high][var]
+        
+        # Values of the perturbed input values
+        highX = dictOfCases[high][var]
+        lowX = dictOfCases[low][var]
+        nominalX = dictOfCases[nominal][var]
+        #print 'highX: ',highX
+        #print 'lowX: ',lowX
+        #print 'nominalX: ',nominalX
+        
+        if cfg['levels'] == 3:
+            # As the perturbations may not be centered on the nominal
+            # condition we caculate the gradient by taking a weighted
+            # average of the forward and backward derivatives. The
+            # weightings are such that the truncation error associated
+            # with this gradient estimate is O(3) or higher (i.e. the
+            # weightings are such that the second order terms in the Taylor
+            # series expansion cancel. Thanks to D.Petty for this theory.)
+            highWeighting = (nominalX - lowX)/(highX - lowX)
+            lowWeighting = (highX - nominalX)/(highX - lowX)
+            
+            # Calculate the aboslute sensitivities
+            #sensitivity_abs[var] = ( highWeighting*(np.array(highValues)-\
+            #                                    np.array(nominalValues))/\
+            #                                   (highX - nominalX) + \
+            #                     lowWeighting*(np.array(nominalValues)-\
+            #                                   np.array(lowValues))/\
+            #                                  (nominalX - lowX)     )
+            sensitivity_abs[var] = {k: ( highWeighting*(highData[k]-\
+                                                nominalData[k]) /\
+                                               (highX - nominalX) + \
+                                          lowWeighting*(nominalData[k]-\
+                                                lowData[k]) /\
+                                               (nominalX - lowX) ) \
+                                    for k in exitVar}
+            #print
+            #print 'exitVar: ',exitVar
+            #print
+            #print 'sens_abs[var]: ',sensitivity_abs[var]
+            
+            # Calculate the relative sensitivities
+            #sensitivity_rel[var] = sensitivity_abs[var]*\
+            #                       nominalX/np.array(nominalValues)
+            sensitivity_rel[var] = {k: sensitivity_abs[var][k]*\
+                                        nominalX/nominalData[k] \
+                                     for k in exitVar}
+            
+            #print 'sens_rel[var]: ',sensitivity_rel[var]
+            
+        elif cfg['levels'] == 5:
+            # For 5 levels per variable we have additional cases
+            # that need to be loaded. Again we do not assume that
+            # the levels are equally spaced around the nominal.
+            # The weightings are such that the truncation error
+            # associated with this estimate is O(4) or higher.
+            tooHigh = 'case'+"{0:02}".format(k)+'3'
+            tooHighData, dontNeed = \
+                   read_outfile('./'+tooHigh+'/'+cfg['resultsFile']) 
+            #tooHighData,dontNeed = \
+            #      read_nenzfr_outfile('./'+tooHigh+'/'+cfg['jobName']+'-exit.stats')
+            #tooHighSupply = read_estcj_outfile('./'+tooHigh+'/'+cfg['jobName']+'-estcj.dat')
+            #tooHighData['supply_rho'] = tooHighSupply['rho']
+            #tooHighData['supply_T'] = tooHighSupply['T']
+            #tooHighData['supply_h'] = tooHighSupply['h']
+            #tooHighData, dontNeed = add_extra_variables(tooHighData, [])
+
+            tooLow = 'case'+"{0:02}".format(k)+'4'
+            tooLowData, dontNeed = \
+                  read_outfile('./'+tooLow+'/'+cfg['resultsFile'])
+            
+            #tooLowData,dontNeed = \
+            #      read_nenzfr_outfile('./'+tooLow+'/'+cfg['jobName']+'-exit.stats')
+            #tooLowSupply = read_estcj_outfile('./'+tooLow+'/'+cfg['jobName']+'-estcj.dat')
+            #tooLowData['supply_rho'] = tooLowSupply['rho']
+            #tooLowData['supply_T'] = tooLowSupply['T']
+            #tooLowData['supply_h'] = tooLowSupply['h']
+            #tooLowData, dontNeed = add_extra_variables(tooLowData, [])
+
+            #print tooHigh, high, nominal, low, tooLow
+
+            # Values of the freestream properties at the perturbed conditions
+            #tooHighValues = get_values(tooHighData,exitVar)
+            #tooLowValues = get_values(tooLowData,exitVar)
+            
+            # Values of the perturbed input values
+            tooHighX = dictOfCases[tooHigh][var]
+            tooLowX = dictOfCases[tooLow][var]
+            tooHighDeltaX = tooHighX - nominalX
+            highDeltaX = highX - nominalX
+            tooLowDeltaX = tooLowX - nominalX
+            lowDeltaX = lowX - nominalX
+            weighting = (tooHighX - tooLowX)/(highX - lowX)
+
+            denom = 1/tooHighDeltaX - 1/tooLowDeltaX -\
+                   (1/highDeltaX - 1/lowDeltaX)*weighting
+
+            #numer =  np.array(tooHighValues)/tooHighDeltaX**2 -\
+            #         np.array(tooLowValues)/tooLowDeltaX**2 -\
+            #        ( np.array(highValues)/highDeltaX**2 -\
+            #          np.array(lowValues)/lowDeltaX**2 )*weighting -\
+            #         np.array(nominalValues)*\
+            #           ( 1/tooHighDeltaX**2 - 1/tooLowDeltaX**2 - \
+            #             ( 1/highDeltaX**2 - 1/lowDeltaX**2 )*weighting )
+
+            numer = {k: tooHighData[k]/tooHighDeltaX**2 -\
+                         tooLowData[k]/tooLowDeltaX**2 -\
+                       ( highData[k]/highDeltaX**2 -\
+                         lowData[k]/lowDeltaX**2 )*weighting -\
+                         nominalData[k]*\
+                         ( 1/tooHighDeltaX**2 - 1/tooLowDeltaX**2 - \
+                           ( 1/highDeltaX**2 - 1/lowDeltaX**2 )*weighting ) \
+                      for k in exitVar}   
+            
+            #numer2values = get_values(numer2,exitVar)
+            #numer2check = numer2values - numer 
+            #print
+            #print 'denom: ',denom
+            #print
+            #print 'exitVar: ',exitVar
+            #print
+            #print 'numer: ',numer
+            #print
+            #print 'numer2: ',numer2
+            #print
+            #print 'numer2check: ',numer2check
+            #print
+            #print
+            
+            # Calculate the absolute sensitivities
+            #sensitivity_abs[var] = numer/denom
+            sensitivity_abs[var] = {k: numer[k]/denom for k in exitVar}    
+            
+            #sens_abs2values = get_values(sensitivity_abs2[var],exitVar)
+            #sens_abs2check = sens_abs2values - sensitivity_abs[var]
+            
+            # Calculate relative sensitivities
+            #sensitivity_rel[var] = sensitivity_abs[var]*nominalX/\
+            #                       np.array(nominalValues)
+            sensitivity_rel[var] = {k: sensitivity_abs[var][k]*\
+                                        nominalX/nominalData[k] \
+                                     for k in exitVar}
+
+            #sens_rel2values = get_values(sensitivity_rel2[var],exitVar)
+            #sens_rel2check = sens_rel2values - sensitivity_rel[var]
+            #
+            #print 'sens_abs[var]: ',sensitivity_abs[var]
+            #print
+            #print 'sens_abs2[var]: ',sensitivity_abs2[var]
+            #print
+            #print 'sens_abs2check[var]: ',sens_abs2check
+            #print
+            #print 'sens_rel2check[var]: ',sens_rel2check
+            #print
+            
+        # Now calculate the uncertainty in each exit flow variable
+        # due to the uncertainty in the current (perturbed)
+        # input variable
+        #uncertainty[var] = sensitivity_rel[var]*inputUncertainties[var]
+        
+        uncertainty[var] = {k: sensitivity_rel[var][k]*\
+                               inputUncertainties[var] for k in exitVar}
+        
+        #uncert2values = get_values(uncertainty2[var],exitVar)
+        #uncert2check = uncert2values - uncertainty[var]
+        #print
+        #print 'exitVar: ',exitVar
+        #print  
+        #print 'uncertainty[var]: ',uncertainty[var]
+        #print
+        #print 'uncertainty2[var]: ',uncertainty2[var]
+        #print
+        #print 'inputUncertainties[var]: ',inputUncertainties[var]
+        #print
+        #print 'uncert2check: ',uncert2check
+        #print
+        
+        
+    # Write out a file of the sensitivities
+    write_sensitivity_summary(sensitivity_rel, perturbedVariables,\
+                              exitVar, 'relative')
+    write_sensitivity_summary(sensitivity_abs, perturbedVariables,\
+                              exitVar, 'absolute')
+    
+    # Write out a file of the uncertainties
+    write_uncertainty_summary(uncertainty, perturbedVariables,\
+                              exitVar, inputUncertainties)
+    
+    return
+
+
+
+def main():
+    """
+    Examine the command-line options and then calculate the sensitivties
+    and uncertainties of each property based on a set of results generated 
+    by "perturb.py".
+    """
+    op = optparse.OptionParser(version=VERSION_STRING)
+    op.add_option('-c', '--config_file', dest='config_file',
+                  help=("filename for the config file"))
+    opt, args = op.parse_args()
+    config_file = opt.config_file
+    #   
+    cfg = {}
+    #
+    if not cfg: #if the configuration dictionary has not been filled up already, load it from a file
+        try: #from Rowan's onedval program
+            execfile(config_file, globals(), cfg)
+        except IOError as e:
+	    print "Error {0}".format(str(e)) 
+            print "There was a problem reading the config file: '{0}'".format(config_file)
+            print "Check that it conforms to Python syntax."
+            print "Bailing out!"
+            sys.exit(1)
+    #
+    run_sensitivity(cfg)
+    #         
+    return
+
+#---------------------------------------------------------------
+
+if __name__ == '__main__':
+    if len(sys.argv) <= 1:
+        print "sensitivity:\n Calculate Sensitivity of Shock Tunnel Test Flow Conditions for a varying inputs"
+        print "   Version:", VERSION_STRING
+        print "   To get some useful hints, invoke the program with option --help."
+        sys.exit(0)
+    return_flag = main()
+    sys.exit(return_flag)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/app/nenzfr/sensitivity_input_utils.py	Wed Jun 03 20:53:51 2020 +1000
@@ -0,0 +1,114 @@
+"""
+sensitivity_input_utils.py -- Function to check the sensitivity input 
+    dictionary that is pulled in from a config file. Default values 
+    are added where appropriate.
+
+.. Author: Luke Doherty (luke.doherty@eng.ox.ac.uk)
+           Oxford Thermofluids Institute
+           The University of Oxford
+
+  Version: 04-May-2020
+"""
+
+import sys, os, copy
+E3BIN = os.path.expandvars("$HOME/e3bin")
+sys.path.append(E3BIN) # installation directory
+sys.path.append("") # so that we can find user's scripts in current directory
+from read_outfile import read_outfile
+from perturb_utils import read_case_summary
+
+#--------------------------------------------------------------------------
+def sensitivity_input_checker(cfg):
+    """Takes the config dictionary and checks that all the necessary
+    inputs are in there. Will bail out if important variables are missing,
+    and set default values for those that are not crucial.
+
+    Returns the config dictionary after it has been checked over and will
+    tell sensitivity to bail out if it finds an issue.
+    
+    """
+
+    print "Checking sensitivity inputs."
+
+    # Make empty input uncertainties dictionary that we will populate
+    #cfg['inputUncertainties'] = {}
+
+    if 'levels' not in cfg:
+        cfg['levels'] = 3
+        print "The number of perturbation levels to use is not specified."
+        print "    Setting it to default value of {0}".format(cfg['levels'])
+
+    if 'outputVariables' not in cfg:
+        print "You have not specified a list of which variables are to be analysed."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+
+    if 'inputVariables' not in cfg:
+        print "You have not specified a list of the input variables (that have been perturbed."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+
+    if 'inputVariablesUncertainty' not in cfg:
+        print "You have not specified a dictionary of the input variables uncertainties."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+
+    if 'typeOfUncertainty' not in cfg:
+        print "You have not specified the type of uncertainty (absolute or relative)."
+        print "    Bailing out."
+        cfg['bad_input'] = True
+    
+    # Check that we have uncertainties for each input variable
+    for variable in cfg['inputVariables']:
+        #print variable
+        #print cfg['inputVariablesUncertainty'].keys()
+        
+        if variable not in cfg['inputVariablesUncertainty'].keys():
+            print "You have not specified an uncertainty for {0}".format(variable)
+            print "    Bailing out."
+            cfg['bad_input'] = True
+    
+    # Read the sensitivity_case_summary file to get the perturbed
+    # variables and their various values
+    perturbedVariables, dictOfCases = read_case_summary()
+    
+    # Define the name of the nominal case and load the exit plane data
+    nominal = 'case000'
+    nominalData, dontNeed = read_outfile('./'+nominal+'/'+cfg['resultsFile'])
+    
+    # Check for discrepancies between the perturbedVariables list and the
+    # inputVariables list
+    for var in perturbedVariables:
+        if var not in cfg['inputVariables']:
+            print "'{0}' was perturbed but does not appear in the 'inputVariables' list.".format(var)
+            # Bailing out could be made unnecessary but it would require 
+            # signifiant changes to how the main sensitivity function 
+            # works since the case number is specific to the variable 
+            # position in the inputVariables list. 
+            print "    Bailing out."
+            cfg['bad_input'] = True
+    
+    # Clean up the input uncertainty. Not sure if this is better here
+    # or in the main function
+    if cfg['typeOfUncertainty'] in ['absolute']:
+        cfg['inputVariablesUncertainty_abs'] = copy.deepcopy(cfg['inputVariablesUncertainty'])
+        cfg['inputVariablesUncertainty_rel'] = {}
+        for var in cfg['inputVariables']:
+            cfg['inputVariablesUncertainty_rel'][var] = \
+                cfg['inputVariablesUncertainty_abs'][var]/\
+                dictOfCases['case000'][var]
+    else:
+        cfg['inputVariablesUncertainty_rel'] = copy.deepcopy(cfg['inputVariablesUncertainty'])
+        cfg['inputVariablesUncertainty_abs'] = {}
+        for var in cfg['inputVariables']:
+            cfg['inputVariablesUncertainty_abs'][var] = \
+                cfg['inputVariablesUncertainty_rel'][var]*\
+                dictOfCases['case000'][var]
+    cfg.pop('inputVariablesUncertainty')
+    cfg.pop('typeOfUncertainty')
+    
+    
+    print "Done checking sensitivity inputs."
+
+    return cfg
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/app/nenzfr/sensitivity_utils.py	Wed Jun 03 20:53:51 2020 +1000
@@ -0,0 +1,174 @@
+"""
+sensitivity_utils.py -- Small untility functions needed by the
+               sensitivity  programs.
+
+.. Author: Luke Doherty (luke.doherty@eng.ox.ac.uk)
+           Oxford Thermofluids Institute
+           The University of Oxford
+
+  Version: 04-May-2020
+"""
+#--------------------------------------------------------------------------
+def get_values(dict, propertyList):
+    """
+    Adapted from Peter Jacob's function "get_fractions" which
+    may be found within "cea2_gas.py".
+    """
+    valueList = []
+    #print propertyList
+    for s in propertyList:
+        #print s
+        if s in dict.keys():
+            #print s
+            #print dict[s]
+            valueList.append(dict[s])
+        else:
+            valueList.append(0.0)
+            print "WARNING: "+s+"was not found in the current case dictionary."
+    return valueList
+
+
+def write_sensitivity_summary(sensitivity, perturbedVariables, \
+                              exitVar, type):
+    """
+    Write out a file summarising the sensitivity of each exit flow 
+    parameter to each of the input parameters.
+    """
+     
+    titleFormatDict = {k : '{0:{fill}>'+"{0}".format(max(11,len(k))+2)+"}" \
+                           for k in perturbedVariables}
+    #print 'titleFormatDict: ',titleFormatDict
+    
+    formatDict = {k : '{0:>'+"{0}".format(max(11,len(k))+2)+".5g}" \
+                      for k in perturbedVariables} 
+    #print 'formatDict: ',formatDict
+    
+    # Old nenzfr-specific header formatting dictionaries
+    #titleFormatDict = {'p1':'{0:{fill}>13}', 'T1':'{0:{fill}>13}',
+    #                   'Vs':'{0:{fill}>13}', 'pe':'{0:{fill}>13}',
+    #                   'Tw':'{0:{fill}>13}', 'BLTrans':'{0:{fill}>13}',
+    #                   'TurbVisRatio':'{0:{fill}>14}',
+    #                   'TurbInten':'{0:{fill}>13}',
+    #                   'CoreRadiusFraction':'{0:{fill}>20}'}
+    #formatDict = {'p1':'{0:13.5g}', 'T1':'{0:>13.5g}',
+    #              'Vs':'{0:>13.5g}', 'pe':'{0:>13.5g}',
+    #              'Tw':'{0:>13.5g}', 'BLTrans':'{0:>13.5g}',
+    #              'TurbVisRatio':'{0:>14.5g}',
+    #              'TurbInten':'{0:>13.5g}',
+    #              'CoreRadiusFraction':'{0:>20.5g}'}
+    
+    # Open the file
+    if type in ['relative']:
+        #fout = open('sensitivities_rel.dat','w')
+        fout = open('sensitivities_rel.dat','w')
+    elif type in ['absolute']:
+        fout = open('sensitivities_abs.dat','w')
+        
+
+    # Write header information
+    #fout.write('{0:}\n'.format(type+' sensitivities'))
+    fout.write('{0:}\n'.format(type+' sensitivities'))
+    #fout.write('{0:>13}'.format('variable'))
+    fout.write('{0:>13}'.format('variable'))
+    #
+    for k in perturbedVariables:
+        #fout.write(titleFormatDict[k].format(k,fill=''))
+        fout.write(titleFormatDict[k].format(k,fill=''))
+    fout.write('\n')
+    #
+    for k in perturbedVariables:
+        #fout.write(titleFormatDict[k].format('-',fill='-'))
+        fout.write(titleFormatDict[k].format('-',fill='-'))
+    fout.write('{0:->13}'.format('-'))
+    fout.write('\n')
+    
+    # Now write out all the data
+    for k in exitVar:
+        #fout.write('{0:>13}'.format(k))
+        fout.write('{0:>13}'.format(k))
+        for kk in perturbedVariables:
+            #X_kk = sensitivity[kk][exitVar.index(k)]
+            X_kk = sensitivity[kk][k]
+
+            fout.write(formatDict[kk].format(X_kk))
+        #
+        fout.write('\n')
+    #
+    fout.close()
+    
+    if type in ['relative']:
+       print 'File "sensitivities_rel.dat" written'
+    elif type in ['absolute']:
+       print 'File "sensitivities_abs.dat" written'
+
+    return    
+    
+def write_uncertainty_summary(uncertainty, perturbedVariables, \
+                              exitVar, inputUncertainties):
+    """
+    Write out a file summarising the sensitivity of each exit flow
+    parameter to each of the input parameters.
+    """
+
+    titleFormatDict = {k : '{0:{fill}>'+"{0}".format(max(8,len(k))+2)+"}" \
+                           for k in perturbedVariables}
+    #print 'titleFormatDict: ',titleFormatDict
+    
+    formatDict = {k : '{0:>'+"{0}".format(max(8,len(k))+2)+".2%}" \
+                      for k in perturbedVariables}
+    #print 'formatDict: ',formatDict
+
+    # Old nenzfr-specific header formatting dictionaries
+    #titleFormatDict = {'p1':'{0:{fill}>10}', 'T1':'{0:{fill}>10}',
+    #                   'Vs':'{0:{fill}>10}', 'pe':'{0:{fill}>10}',
+    #                   'Tw':'{0:{fill}>10}', 'BLTrans':'{0:{fill}>10}',
+    #                   'TurbVisRatio':'{0:{fill}>14}',
+    #                   'TurbInten':'{0:{fill}>11}',
+    #                   'CoreRadiusFraction':'{0:{fill}>20}'}
+    #formatDict = {'p1':'{0:10.2%}', 'T1':'{0:>10.2%}',
+    #              'Vs':'{0:>10.2%}', 'pe':'{0:>10.2%}',
+    #              'Tw':'{0:>10.2%}', 'BLTrans':'{0:>10.2%}',
+    #              'TurbVisRatio':'{0:>14.2%}',
+    #              'TurbInten':'{0:>11.2%}',
+    #              'CoreRadiusFraction':'{0:>20.2%}'}
+
+    fout = open('uncertainties.dat','w')
+    
+    # Write header information
+    fout.write('{0:>19}'.format('variable'))
+    for k in perturbedVariables:
+        fout.write(titleFormatDict[k].format(k,fill=''))
+    fout.write('{0:>10}'.format('total'))
+    fout.write('\n')
+    
+    # Write out the uncertainty for each perturbed variable
+    fout.write('{0:>19}'.format('(input uncertainty)'))
+    for k in perturbedVariables:
+        fout.write(formatDict[k].format(inputUncertainties[k]))
+    fout.write('{0:>10}'.format(''))
+    fout.write('\n')
+    for k in perturbedVariables:
+        fout.write(titleFormatDict[k].format('-',fill='-'))
+    fout.write('{0:->29}'.format('-'))
+    fout.write('\n')    
+    
+    # Now write out all the data. We also calculate and write out the 
+    # total uncertainty in each exit flow variable due to the contributions
+    # from all input values.
+    for k in exitVar:
+        fout.write('{0:>19}'.format(k))
+        X_total = 0.0
+        for kk in perturbedVariables:
+            #X_kk = uncertainty[kk][exitVar.index(k)]
+            X_kk = uncertainty[kk][k]
+            fout.write(formatDict[kk].format(X_kk))
+            X_total += X_kk**2 # Calculate total
+        fout.write('{0:10.2%}'.format(X_total**0.5))
+        fout.write('\n')
+
+    fout.close()
+
+    print 'File "uncertainties.dat" written'
+    return
+
+