Project 1 - Parsing NASA thermodynamic data

AIM

• To parse NASA thermodynamic data from THERM.dat file using MATLAB
• To calculate molecular weight of every species and displaying in command window.
• To plot specific heat, enthalpy and entropy with temperature range distributed from global high temperature to global low temperature.
• To create a separate folder for all the plots and saving them in ‘.png’ format.

THEORY

The official record of thermodynamic values at different temperatures has been stored by NASA in the form of polynomials with specified constants ranging from one to fourteen for various species. In order to quantify the values, the polynomials can be solved with defined temperature values in accordance with our requirements.

The term ‘parsing’ refers to MATLAB reading an m-file (script or function) and translating it into an internal code that it stores in memory. It is this ‘parsed’ file in memory that is actually run when you execute the code.

GOVERNING EQUATIONS

Specific Heat,

`c_p/R = a_1 +a_2T + a_3T^2 + a_4T^3 +a_5T^4`

Enthalpy,

`H/(RT) = a_1 + a_2(T/2) + a_3(T^2/3) + a_4(T^3/4) +a_5(T^4/5) +a_6/T`

Entropy

`S/R =a_1 lnT + a_2T + a_3T^2/2+a_4T^3/3+a_5T^4/4+a_6`

Where,

`c_p`= Specific Heat

H = Enthalpy

S = Entropy

R = Universal Gas Constant (8.314 J/mol.K)

T = Temperature

*The first 7 coefficients are High Temperature Coefficients and the second 7 coefficients are Low Temperature Coefficients

MATLAB CODE

The code is divided into 7 different files. One is for the main program and the other 6 are supporting functions.

Function files:

1. File 1 - plot_data.m

% Function Definition: Plot thermodynamic data for specific species. 
function []=plot_data(f1,species,init_position,plot_choice)

% Initialize while loop status.
status=true;

% Convert User Input to Upper Case.
species = upper(species);

% Seek the cursor position to initial read position.
fseek(f1,init_position,-1);

% Search the species name in THERMO.dat file.
while status==true
    for i=1:4
        x{i}=fgetl(f1);
        if strcmp(x{1},'END')
			% Print error and return to main program if species not found.
            fprintf('n')
            fprintf('Sorry! Thermodynamic data not available for %sn',species)
            status=false;
            return;
        end
    end

% If species is found in THERMO.dat file trim the leading and trailing spaces and get the name.
    x1=strtrim(x{1});
    x1=strsplit(x1);
    s1=x1{1};

% If species is found in THERMO.dat file break out of the while loop.
    if string(s1)==string(species)
        break
    end
end

% Extract the local temperature ranges from the first line.
lt_low=str2num(x1{end-3});
lt_mid=str2num(x1{end-1});
lt_high=str2num(x1{end-2});

% Define Temperature arrays for plotting.
T_l = lt_low:1:lt_mid;
T_h = lt_mid:1:lt_high;
T = [T_l,T_h];

%Extract polynomial coefficients from remaining lines.

% High temperature range coefficients (a1-a5) from Line 2
x2=x{2};
a=findstr(x2,'E');
a1h=str2num(x2(1:a(1)+3));
a2h=str2num(x2(a(1)+4:a(2)+3));
a3h=str2num(x2(a(2)+4:a(3)+3));
a4h=str2num(x2(a(3)+4:a(4)+3));
a5h=str2num(x2(a(4)+4:a(5)+3));

% Coefficients a6-a7 of high range and a1-a3 of low range from line 3
x3=x{3};
a=findstr(x3,'E');
a6h=str2num(x3(1:a(1)+3));
a7h=str2num(x3(a(1)+4:a(2)+3));
a1l=str2num(x3(a(2)+4:a(3)+3));
a2l=str2num(x3(a(3)+4:a(4)+3));
a3l=str2num(x3(a(4)+4:a(5)+3));

% Coefficients a4-a7 low range from line 4.
x4=x{4};
a=findstr(x4,'E');
a4l=str2num(x4(1:a(1)+3));
a5l=str2num(x4(a(1)+4:a(2)+3));
a6l=str2num(x4(a(2)+4:a(3)+3));
a7l=str2num(x4(a(3)+4:a(4)+3));

% Define coefficients array.
coeffs = [a1h a2h a3h a4h a5h a6h a7h a1l a2l a3l a4l a5l a6l a7l];

% Call calculat_properties() function to calculate Cp,H & S arrays.
[Cp,H,S] = calculating_thermodynamic_properties(T_l,T_h,coeffs);

% Get the current folder directory in cf variable using pwd command.
cf = pwd;

% Location of images to be saved
d = strcat('D:\Matlab Challenges',getenv('USERNAME'),'Matlab Programs');

% Define the sub-directory to save plot images.
plot_directory = 'NASA Thermodynamic Data Plots'+"" + species;

% Do not create the folders if they already exist.
if ~exist(plot_directory,'dir')
    mkdir(plot_directory);
end

% Plot directories with file names
s1 = plot_directory + 'Specific Heat Vs Temperature Plot '+species+'.png';
s2 = plot_directory + 'Enthalpy Vs Temperature Plot '+species+'.png';
s3 = plot_directory + 'Entropy Vs Temperature Plot '+species+'.png';

% String for title of plot.
t1 = ("Specific Heat Vs Temperature Plot for "+ species);
t2 = "Enthalpy Vs Temperature Plot for "+ species;
t3 = "Entropy Vs Temperature Plot for "+ species;

% Plot the thermodynamic properties

% Plot Specific Heat Vs Temperature.
if plot_choice==1||plot_choice==4
fig1=figure(1);
fig1.WindowState='maximized';
Plot_Cp = plot(T,Cp,'linewidth',1.5);
grid on
title(t1);
xlabel('Temperature in Kelvin');
ylabel('Specific Heat in kJ/kg-K');
saveas(fig1,s1);
disp("Successfully saved the plot image to:  "+cf+''+string(s1));
end

% Plot Enthalpy Vs Temperature
if plot_choice==2||plot_choice==4
fig2=figure(2);
fig2.WindowState='maximized';
Plot_H = plot(T,H,'linewidth',1.5);
grid on
xlabel('Temperature in Kelvin');
ylabel('Enthalpy in KJ/mol');
title(t2);
saveas(fig2,s2);
disp("Successfully saved the plot image to:"+cf+''+string(s2));
end

% Plot Entropy Vs Temperature.
if plot_choice==3||plot_choice==4
fig3=figure(3);
fig3.WindowState='maximized';
Plot_S = plot(T,S,'linewidth',1.5);
grid on
xlabel('Temperature in Kelvin');
ylabel('Entropy in KJ/mol');
title(t3);
saveas(fig3,s3);
disp("Successfully saved the plot image to:"+cf+''+string(s3));
end

% Calculate Molar Mass using calculating_molar_mass() function
MM = calculating_molar_mass(species,1);
fprintf("nMolar Mass of %s is %fn",species,MM);
fprintf('n');

%Save all the Therodynamic properties to a file
% File name with plot directory
dat_d = strcat(plot_directory,'Species_Data.dat');

%Open the file for writing data
f2 = fopen(dat_d,'w');

% Write the data to the file
fprintf(f2,'%15s %15s %15s %15sn','Temperature','Specific-Heat','Enthalpy','Entropy');

for i=1:length(T)
    fprintf(f2,'%15.5f %15.5f %15.5f %15.5fn',T(i),Cp(i),H(i),S(i)');
end

% Close the file
fclose(f2);

% End of function.
end

Code Explanation:

Line 2 is the function definition.

Input variables

f1 – is the file ID of THERM.dat

species – Name of species to plot data for in string.

initial_position – position of read cursor from which the name of first species begins.

plot_choice – choice of which properties to plot.

No output variables for this function.

Line 5 -8 – Initialize while loop status to true and capitalize the string characters using upper function.

Line 11 sets the cursor read position to initial_position using fseek() function.

Line 14 to 35 searches the name of species in the THERMO.dat file.

Line 14 to 16: While and for loop for searching the species name in data file. For each iteration of while loop the for loop iterates 4 times and each time it reads a line from the THERMO.dat file and stores it in an array x in line 16.

Line 17 to 22: Checks if end of the file is reached and it means that the requested species name is not found. Line 19, 20 prints an error if the species name is not found. Line 22 returns the control to the main program.

Line 23 end of if statement & Line 24 is end of for loop.

Line 27 to 29: The name of species is in the first element of array x. Line 27 trims the leading and trailing spaces in the first element of x using strtrim() function and assign it to variable x1. Line 28 splits string x1 into words using strsplit. Line 29 first word of x1 is the species name and is stored in s1.

Line 32 s1 is compared with input species, if a match is found line 33 breaks out of the while loop.

Lines 38 to 40: Gets the local temperature from the end of x1 array using indexing from the end and converts it to numbers and assign it to lt_low, lt_mid & lt_high variables.

Line 43 & 44: Converts temperature range values to arrays with intervals of 1 degree celcius.

Line 45: Assembles low range and high range temperature arrays into one single array.

Lines 50 to 56: Extracts polynomial coefficients from strings.

Line 50: extracts second line from string array from x and assigns it to x2.

Line 51: Searches for position of letter ‘E’ in the second line using findstr(x2, E) method. First input is x2 which is the string in which search operation has to be done and second is ‘E’ the character to be searched. Variable a is the array which contains the index position of ‘E’ in x2.

Line 52 to 56: uses string slicing with index values of ‘E’ (which is stored in a) to get the coefficients. The coefficients are converted to number using str2num() function.

Similarly coefficients from the 3^rd & 4^th lines are extracted in lines 59 to 65 and 68 to 73.

Line 76 converts all the coefficients into an array coeffs.

Line 79 calls the calculate properties method to calculate the Cp, H and S values for the corresponding temperature values in the form of an array.

Line 82 gets the address of current folder and saves it in cf variable using pwd command.

Note: This program will work on windows platform only because the address to desktop might be different on other platforms.

Line 91 defines the plot directory in a string variable plot_directory which contains the address of sub-folders to be created for each species.

Line 94 checks if the directory exists if not then line 95 creates required directory.

Line 96 is end of if statement.

Lines 99 to 101 defines the file names for images along with the directory in which these should be saved for 3 thermodynamic properties.

Lines 104 to 106 defines the plot title lines for the plots.

Line 111 is the if condition for plotting the specific heat. If user choses 1 or 4 option statements under this if condition (lines 112-120) are executed.

Line 112 opens a figure for plotting. Line 113 changes the WindowState variable of figure object to ‘maximized’ which maximizes the window.

Line 114 plots specific heat vs temperature using plot command.

Line 115 switches on the grid line 116 sets the title, Lines 117 & 118 sets the x & y label.

Line 119 saves the figure as an image where fig1 is the variable that holds the address of figure object and s1 is the sub-directory with image file name.

Line 120 displays a message to user that the plot was successfully saved to the directory. Line 121 is end of if statement.

Similarly lines 123 to 134 & lines 137 to 147 plot and save images for enthalpy and entropy respectively.

Line 150 changes the current directory back to the old one where the function calculating_molar_mass() can be found.

Line 151 calls the calculating_molar_mass() function to calculate the molar mass of the species.

Line 151 & 152 prints the calculated molar mass.

Lines 157 to 175 saves the thermodynamic data into a DAT file.

Line 157 changes the directory back tto the plot directory.

Line 159 creates the file name with the directory.

Line 162 opens the file for writing using fopen command.

Line 165 prints the header information into the DAT file.

Line 167 to 169 is the for loop for writing the data into the file line by line. Each iteration writes one line and the for loop variable i is used for indexing the temperature, Cp, H & S arrays.

Line 172 closes the file using fclose() command.

Line 175 changes the directory back to the old one.

Line 177 is the end of function.

2. File 2 - plot_all_data.m

% Function definition: Plot thermodynamic data for all species.
function []=plot_all_data(f1,init_position)

% Initialize while loop status to true.
status=true;

% Set cursor read position to initial position.
fseek(f1,init_position,-1);

% While loop: 1 execution for each species.
while status==true

% For loop: 4 iterations to extract coefficients of 1 species.
 
for i=1:4
	% Extract data using fgetl(fileID).
    x{i}=fgetl(f1);
	
	% Returns the control to the calling function if end of file is reached.
    if strcmp(x{1},'END')
        close(fig1);
        status=false;
        return;
    end
end

% Extract the name of species: first word first line.
x1=strtrim(x{1});
x1=strsplit(x1);
species=x1{1};

% Local Temperature range
lt_low=str2num(x1{end-3}); % Local Low Temperature.
lt_mid=str2num(x1{end-1}); % Local Mid Temperature.
lt_high=str2num(x1{end-2}); % Local High Temperature.
lt=[lt_low lt_mid lt_high];

T_l = lt_low:1:lt_mid; % Local Low Temperature Range.
T_h = lt_mid:1:lt_high; % Local High Temperature Range.
T = [T_l,T_h]; % Local Temperature Range.

x2=x{2}; % 2nd line Extraction from cell array to string.
a=findstr(x2,'E');% Search for character 'E' in x2.
a1h=str2num(x2(1:a(1)+3)); %Coefficient a1 high temperature range.
a2h=str2num(x2(a(1)+4:a(2)+3));%Coefficient a2 high temperature range.
a3h=str2num(x2(a(2)+4:a(3)+3));%Coefficient a3 high temperature range.
a4h=str2num(x2(a(3)+4:a(4)+3));%Coefficient a4 high temperature range.
a5h=str2num(x2(a(4)+4:a(5)+3));%Coefficient a5 high temperature range.
x3=x{3}; % Extract 3rd line to x3
a=findstr(x3,'E');% Search for character 'E' in x3
a6h=str2num(x3(1:a(1)+3));%Coefficient a6 high temperature range.
a7h=str2num(x3(a(1)+4:a(2)+3));%Coefficient a7 high temperature range.
a1l=str2num(x3(a(2)+4:a(3)+3));%Coefficient a1 low temperature range.
a2l=str2num(x3(a(3)+4:a(4)+3));%Coefficient a2 low temperature range.
a3l=str2num(x3(a(4)+4:a(5)+3));%Coefficient a3 low temperature range.
x4=x{4};% Extract 4th line from cell array x to x3.
a=findstr(x4,'E');% Search for character 'E' in x3.
a4l=str2num(x4(1:a(1)+3));%Coefficient a4 low temperature range.
a5l=str2num(x4(a(1)+4:a(2)+3));%Coefficient a5 low temperature range.
a6l=str2num(x4(a(2)+4:a(3)+3));%%Coefficient a6 low temperature range.
a7l=str2num(x4(a(3)+4:a(4)+3));%Coefficient a7 low temperature range.

% Assemble all coefficients in one array named coeffs.
coeffs = [a1h a2h a3h a4h a5h a6h a7h a1l a2l a3l a4l a5l a6l a7l];

% Calculate Cp, H, & S using calculate_properties() function.
[Cp,H,S] = calculating_thermodynamic_properties(T_l,T_h,coeffs);

% Get the current folder directory in cf variable using pwd command.
cf=pwd;

% Location of images to be saved
d=strcat('D:Users',getenv('USERNAME'),'Karan');

% Define the string path for plot directory
plot_directory = 'NASA Thermodynamic Data Plots'+"" + species;

% Create the folders if they do not exist.
if ~exist(plot_directory,'dir')
    mkdir(plot_directory);
end

% Open figure window and maximize it
fig1=figure(1);
fig1.WindowState='maximized';

% Image file names with directories.
s1 = plot_directory + 'Specific Heat Vs Temperature Plot '+species+'.png';
s2 = plot_directory + 'Enthalpy Vs Temperature Plot '+species+'.png';
s3 = plot_directory + 'Entropy Vs Temperature Plot '+species+'.png';

% Title string for plots.
t1 = ("Specific Heat Vs Temperature Plot for "+ species);
t2 = "Enthalpy Vs Temperature Plot for "+ species;
t3 = "Entropy Vs Temperature Plot for "+ species;

% Specific Heat Vs Temperature plot.
Plot_Cp = plot(T,Cp,'linewidth',1.5);
grid on
title(t1);
xlabel('Temperature in Kelvin');
ylabel('Specific Heat in kJ/kg-K');
saveas(fig1,s1); % Save plot to specified directory.

% Enthalpy Vs Temperature Plot
Plot_H = plot(T,H,'linewidth',1.5);
grid on
xlabel('Temperature in Kelvin');
ylabel('Enthalpy in KJ/mol');
title(t2);
saveas(fig1,s2);% Save plot to specified directory.

% Entropy Vs Temperature Plot. 
Plot_S = plot(T,S,'linewidth',1.5);
grid on
xlabel('Temperature in Kelvin');
ylabel('Entropy in KJ/mol');
title(t3);
saveas(fig1,s3);% Save plot to specified directory.

% Display the directory of saved plots to command window.
Message="Successfully saved the plot images for "+string(species)+" to:  "+cf+""+plot_directory;
disp(Message);
% Calculate the molar mass using Molar_Mass() function.
MM = calculating_molar_mass(species,1);

% Print the molar mass to command window.
fprintf("Molar Mass of %s is %fnn",species,MM);

%Save all the Therodynamic properties to a file
% File name with plot directory
dat_d = strcat(plot_directory,'Species_Data.dat');

%Open the file for writing data
f2 = fopen(dat_d,'w');

% Write the data to the file
fprintf(f2,'%15s %15s %15s %15sn','Temperature','Specific-Heat','Enthalpy','Entropy');

for i=1:length(T)
    fprintf(f2,'%15.2f %15.5f %15.5f %15.5fn',T(i),Cp(i),H(i),S(i)');
end

% Close the file
fclose(f2);

% Change back the directory
cd(cf);
% End of loop
end
% End of function.
end

Code Explanation:

Line 2 is the function definition.

Input variables

f1 – is the file ID of THERM.dat

initial_position – position of read cursor from which the name of first species begins.

No output variables for this function.

Line 5 initializes while loop status to true.

Line 8 sets the cursor read position to initial_position using fseek() function.

Line 11 defines the while loop, while loop iterates once for each species.

Line 15 defines the for loop. The for loop iterates 4 times for each species.

Line 17 reads the line from the THERMO.dat file and assigns it to x.

Line 20 compares x with the string ‘END’ to check whether it is the end of the file.

Line 21 closes the figure. Line 22 sets the while loop status to false. Line 23 returns the function control.

Line 24 is end of if statement. Line 25 is end of for loop.

Line 28 trims the leading and trailing spaces from the first line read in every first iteration of the for loop. Line 29 splits the line into words array & Line 30 extracts the first word and assigns it to species variable.

Line 33 extracts the last 3^rd word (which is the local low temperature) from line 1 by indexing the string array and converts it to number using str2num() function. Similarly, lines 34 & 35 extract local mid and high temperatures respectively.

Line 38 & 39 Converts the temperature values in array of low temperature and high temperature range with intervals of 1 degree celcius.

Line 40 Assembles low range and high range temperature arrays into one single array.

Lines 42 to 48: Extracts polynomial coefficients from strings.

Line 42: extracts second line from string array from x and assigns it to x2.

Line 43: Searches for position of letter ‘E’ in the second line using findstr(x2, E) method. First input is x2 which is the string in which search operation has to be done and second is ‘E’ the character to be searched. Variable a is the array which contains the index position of ‘E’ in x2.

Line 44 to 48: uses string slicing with index values of ‘E’ (which is stored in a) to get the coefficients. The coefficients are converted to number using str2num() function.

Similarly coefficients from the 3^rd & 4^th lines are extracted in lines 49 to 55 and 56 to 61.

Line 64 converts all the coefficients into an array coeffs.

Line 67 calls the calculate properties method to calculate the Cp, H and S values for the corresponding temperature values in the form of an array.

Line 70 Gets the current folder directory and stores in variable cf using pwd command.

Line 73 creates a string variable d which contains the address of Desktop folder of the current user using getenv(‘USERNAME’) command.

Line 76 changes the current directory to desktop using cd command.

Line 79 defines the plot directory in a string variable plot_directory.

Line 82 checks if the directory exists if not then line 83 creates required directory.

Line 84 is end of if statement.

Line 87 opens a figure window and Line 88 maximizes the window by changing the WindowState variable to ‘maximized’.

Lines 91 to 93 creates string variables for holding the image file names with their respective directories.

Lines 96 to 98 creates string variables for title of plots.

Line 101 plots the specific heat vs temperature plot using the plot command.

Line 102 switches on the grid, line 103 sets the plot title, lines 104 & 105 sets the x & y labels respectively, line 106 saves the plot as image using saveas command where fig1 is the variable for figure and s1 is the image name with path.

Similarly, enthalpy and entropy vs temperature are plotted in lines 109 to 114 & lines 117 to 122.

Line 150 creates a message variable to display in command line that the plot was successfully saved.

Line 151 prints the message variable.

Line 128 changes the current folder to the old one using cf variable so that calculating_molar_mass() function can be found by MATLAB.

Line 130 calls the user defined function calculating_molar_mass() to compute the molar mass of species.

Line 133 prints the molar mass to the command window.

Line 136 is end of while loop.

Lines 137 to 155 writes the thermodynamic data to a DAT file and saves it in the plot directory.

Line 137 changes the directory back to plot directory.

Line 139 creates the file name with the directory for storing the data in a DAT file.

Line 142 opens a file for writing using the fopen command.

Line 145 prints the header information into the DAT file using fprintf() function.

Lines 147 to 149 is the for loop for writing the data into the file line by line. Each iteration writes one line and the for loop variable i is used for indexing the temperature, Cp, H & S arrays.

Line 152 closes the file using the fclose command.

Line 155 changes back the directory to the old one.

Line 157 is the end of while loop.

Line 159 is the end of function.

3. File 3 – calculating_molar_mass.m

% Function definition: Calculates the molar mass of species from the species name.
% This function can be used even if more complex species ( Heavier carbon compounds with long carbon chains where multiple brackets are required to write the chemical formula) in the THERMO.dat file.
function [results] = calculating_molar_mass(species,multiplier_out)

% Convert species to uppercase character array and replace argon element symbol 'AR' to just 'A'
species=char(strrep(upper(species),'AR','A'));
z=species;

% Return zero for empty string.
if strcmp(species,'');
    result=0;
    return
end

% Define Element Symbol characters and their atomic weights
Species_Array = ['O' 'H' 'C' 'S' 'N' 'A'];
Atomic_Weight = [15.999,1.00784,12.0107,32.065,14.0067,39.948];

%Initialize Molar Mass to 0.
M=0;

% If the species contains round braces, extract the contents outside the braces
% Using function extract_contents_outside_brackets().
if contains(species,'(')
    Outside_Brackets = extract_contents_outside_brackets(species);
    z=Outside_Brackets;
end

% Calculate the Molar Mass for contents outside brackets.
for i=1:length(z)
	%Check if the character is number, if yes assign it to n.
    [n,status]=str2num(z(i));
	
	% Check if the character is letter, if yes check which element it is and add its molar mass
    if isletter(z(i))
        ch = z(i);
        index = strfind(Species_Array,ch);
        M = M + Atomic_Weight(index);
    
	% If the character is number then increment the molar mass by number-1 times atomic weight of previous element 
	elseif status
        M = M + Atomic_Weight(index)*(n-1);       
    end
% End of loop
end

% Multiply the molar mass with the supplied multiplier.
M=M*multiplier_out;% This molar mass is only the molar mass of those elements which are outside the brackets

% Calculate the molar mass of contents inside brackets with multiplier.
if contains(species,'(')
    % Extract contents inside brackets
	[Inside_Brackets,multiplier_in] = extract_contents_inside_brackets(species);
	
	% Add the molar mass of inside brackets elements. 
    M = M + calculating_molar_mass(Inside_Brackets,(multiplier_in*multiplier_out));
end

% Assign M to results.
results=M;% Function definition: Calculates the molar mass of species from the species name.
% This function can be used even if more complex species ( Heavier carbon compounds with long carbon chains where multiple brackets are required to write the chemical formula) in the THERMO.dat file.
function [results] = calculating_molar_mass(species,multiplier_out)

% Convert species to uppercase character array and replace argon element symbol 'AR' to just 'A'
species=char(strrep(upper(species),'AR','A'));
z=species;

% Return zero for empty string.
if strcmp(species,'');
    result=0;
    return
end

% Define Element Symbol characters and their atomic weights
Species_Array = ['O' 'H' 'C' 'S' 'N' 'A'];
Atomic_Weight = [15.999,1.00784,12.0107,32.065,14.0067,39.948];

%Initialize Molar Mass to 0.
M=0;

% If the species contains round braces, extract the contents outside the braces
% Using function extract_contents_outside_brackets().
if contains(species,'(')
    Outside_Brackets = extract_contents_outside_brackets(species);
    z=Outside_Brackets;
end

% Calculate the Molar Mass for contents outside brackets.
for i=1:length(z)
	%Check if the character is number, if yes assign it to n.
    [n,status]=str2num(z(i));
	
	% Check if the character is letter, if yes check which element it is and add its molar mass
    if isletter(z(i))
        ch = z(i);
        index = strfind(Species_Array,ch);
        M = M + Atomic_Weight(index);
    
	% If the character is number then increment the molar mass by number-1 times atomic weight of previous element 
	elseif status
        M = M + Atomic_Weight(index)*(n-1);       
    end
% End of loop
end

% Multiply the molar mass with the supplied multiplier.
M=M*multiplier_out;% This molar mass is only the molar mass of those elements which are outside the brackets

% Calculate the molar mass of contents inside brackets with multiplier.
if contains(species,'(')
    % Extract contents inside brackets
	[Inside_Brackets,multiplier_in] = extract_contents_inside_brackets(species);
	
	% Add the molar mass of inside brackets elements. 
    M = M + calculating_molar_mass(Inside_Brackets,(multiplier_in*multiplier_out));
end

% Assign M to results.
results=M;

Code Explanation:

Line 3 is the function definition. Calculates the molar mass of species from the species name. This function can be used even if more complex species ( Heavier carbon compounds with long carbon chains where multiple brackets are required to write the chemical formula) in the THERMO.dat file.

Inputs:

species – Chemical formula of species as string. Multiple brackets can be used if there are multiple group of compounds like sulphate, carbonate, nitrate. Currently only 6 elements are defined in this function.

multiplier_out – This function is a recursive function i.e. a function which calls itself until it finds the end result. The function uses a different multiplier value based on the numbers in the chemical formula when it calls itself, but when an outside program call this function multiplier_out should be 1.

Outputs:

results – molar mass of the species.

Line 6 reconditions the species variable using 3 functions. One is upper() which capitalizes the characters in it, second is strrep() finds and replaces ‘AR’ with ‘A’ to designate element argon with a single letter & third is char() which converts string array into character array in case someone passes a string array variable.

Line 7 assigns value of string to a new variable z.

Line 10 to 13 checks if the string is empty, if yes it returns the control of function.

Lines 16 & 17 Defines Element Symbol characters and their atomic weights in two different arrays.

Line 20 initializes the molar mass variable to 0.

Lines 24 to 27 checks if the chemical formula contains any brackets. If yes, the user defined function extract_contents_outside_brackets() extracts the contents outside the first open bracket and last close bracket and assigns it to z.

Lines 30 to 45 calculates the molar mass of elements in variable z.

Line 30 is for loop, iterates for every character in z.

Line 32 checks if the character is a number or not and changes the status variable to true if it is a number. And it converts from string to number.

Line 35 checks if the character is a letter or not. If it is a letter then lines 36 to 38 are executed.

Line 37 finds the character in the Species_Array to check which element it is using strfind() method. If found it returns an index of the character position in the Species_Array.

Line 38 adds the atomic weight of that species to molar mass variable.

Line 41 is the elseif statement, if status is true the character is a number then the atomic weight of the previous character is multiplied with the number minus one and added it to molar mass variable. One is subtracted because in the previous iteration the atomic mass of that element was already added.

Line 43 is the end of ifelse statement.

Line 45 is the end of for loop.

Line 48 multiplies the multiplier_out variable to the molar mass. This makes no difference during the first call of calculating_molar_mass() function. During the second call, which the function itself makes, this statement multiplies the number immediately after the bracket to the molar mass of contents inside the bracket.

Line 51 check if the species contains brackets. If yes it executes lines 53 and 56.

Line 53 extracts the contents of species variable inside the first open bracket and last close bracket and assigns it to inside brackets variable using the function extract_contents_inside_bracket(). It also assigns any numbers immediately after the last close bracket to multiplier_in variable.

Line 56 calls the calculating_molar_mass() function itself to calculate the molar mass of the contents inside the first open bracket & last close bracket and the process repeats until the last bracket pair.

Line 57 ends the if statement.

Line 60 assigns the molar mass to result variable.

4. File 4 - extracting_contents_outside_brackets.m

% Function definition: This function extracts all the contents outside the first open bracket and last close bracket, except for any numbers if any, immediately after the last close bracket.
function [results] = extract_contents_outside_brackets(species)

%If the input does not contain any brackets then return ''.
if ~contains(species,'(')
    results='';
    return
end
% Assign species to x.
x=species;
% Search for start brackets.
a=strfind(x,'(');
% Get the index of first start bracket.
a=a(1);
% Search for end brackets.
b=strfind(x,')');
% Get the index of last close bracket.
b=b(end);

% Assign all characters before first start bracket to x1 and all characters after last close bracket to x2.
x1 = x(1:a-1);
x2 = x(b+1:end);
% Initialize x3.
x3='';

% Check if x2 starts with numerical characters, If yes, remove them all and assign the resulting string to x3.
for i=1:length(x2);
    [m,status]=str2num(x2(i));
    
    if status;
        x3=strip(x2,'left',x2(i));
    else
        break
    end
end

% Concantinate x1 & x3 and return the resulting string.
results=strcat(x1,x3);

Code Explanation:

Line 2 is the function definition

Inputs

species – takes string input, chemical formula of species

output – returns the chemical formula of all elements outside the first open bracket and last close bracket.

Line 5 checks if species contains any brackets. If no, the line 6 sets the result to ‘’ and line 7 returns the control of the function. Line 8 is end of if statement.

Line 10 assigns contents of species to variable x.

Line 12 searches for open bracket in x using strfind() function and assigns the indexes to a.

Line 14 assigns the first element (i.e. index of first open bracket) of a to a itself.

Line 16 searches for close bracket in x using strfind() function and assigns it to b.

Line 18 assigns the last element (index of last close bracket) of b to b itself.

Line 21 slices variable x and assigns x1 the characters before the first open bracket.

Line 22 slices variable x and assigns x2 the characters after the last close bracket

Line 24 initializes x3 to ‘’.

Line 27 to 35 removes any integers in the beginning of x2 string.

Line 27 is the for loop which iterates for every character in x2.

Line 28 checks if the character is a number. It changes the status variable to true and converts string to number and assigns it to variable m.

Line 30 is the if statement, if the status is true (i.e. character is a number) then it removes the character from x2 and assigns it to x2.

Line 32 & 33 otherwise if it is a character, it breaks out of the for loop.

Line 35 is end of the for loop.

Line 38 assigns the results variable the concatenation of x1 & x3.

5. File 5 – extract_contents_inside_brackets.m

% Function definition: this function extracts string inside first open bracket and last close bracket.
% and the resulting string is stored in results and any integer after the close bracket is stored in multiplier variable.
function [results n] = extract_contents_inside_brackets(species)

% If the input does not contain any brackets then set results='' & n=1 return the function
if ~contains(species,'(')
    results='';
    n=1;
    return
end

% Assign string to x
x=species;
% Search for open bracket.
a=strfind(x,'(');
% Get the index of first open bracket.
a=a(1);

%Search for close bracket.
b=strfind(x,')');

% Get the index for last close bracket.
b=b(end);

% Slice the contents inside brackets and store it in results variable.
results=x(a+1:b-1);

% Get the remaining characters after the last close bracket and store it in y.
y = species((b+1):end);

%Initialize n to ''.
n='';

% Search for any numbers immediately after the last close bracket.
for i=1:length(y)
	%Check for numbers.
    [m,status]=str2num(y(i));
    
	%If first one itself is a alphabet then break the loop.
    if isletter(y(i))
        n='1';
        break;
	% If it is a number then convert it to string and concantinate it to n.
    elseif status
        m=num2str(m);
        n=strcat(n,m);
    end
end
% If there are no numbers set n to '1'.
if strcmp(n,'')
    n='1';
end

% Convert n to number.
n=str2num(n);

Code Explanation:

Line 3 is the function definition. This function extracts string inside first open bracket and last close bracket and the resulting string is stored in results and any integer after the close bracket is stored in multiplier variable.

Line 6 checks if the species contains any brackets, if not then line 7 & 8 sets the results variable to ‘’ and n to 1 and line 9 returns the control to the calling function.

Line 10 ends the if statement.

Line 13 assigns species to x.

Line 15 searches for open brackets in x and returns the array of bracket position (or index) to variable a.

Line 17 assigns the first value (position of first open bracket) in array to a itself.

Line 20 searches for close bracket in x and assigns the array of bracket position to b.

Line 23 assigns the last element in b (index of last close bracket) to b itself.

Line 26 slices the variable x using indexes a & b to get the contents inside the first open bracket and last close bracket.

Line 29 assigns all the characters after the last close bracket to variable y using string slicing methods.

Line 32 initializes n variable to ‘’.

Line 35 to 48 Searches for any numbers immediately after the last close bracket.

Line 35 is the for loop which iterates once for every character in y.

Line 37 checks if the character is numeric. If yes it sets status variable to true and converts it to number and assigns it to m.

Line 40 checks if the character is an alphabet. On finding the first alphabet it breaks out of the for loop in line 42 and assigns n = ‘1’.

Line 44 if the character is numeric statements under elseif (line 45 & 46) executes.

Line 45 & 46: if character is number m is converted to number and concatenated to n.

Line 47 is end of elseif statement.

Line 48 is end of for loop.

Line 50 to 52: If there are no numbers after the last close bracket then the multiplier n is assigned 1.

Line 55 converts n to number.

6. File 6 – calculating_thermodynamic_properties.m

% Function Definition: Calculates Thermodynamic Properties from polynomial coefficients and Temperature values.
function [Cp H S] = calculating_thermodynamic_properties(T_l,T_h,coeffs)

% Universal Gas Constant.
R = 8.314;

% Low Temperature Range polynomial coefficients.
a1l = coeffs(8);
a2l = coeffs(9);
a3l = coeffs(10);
a4l = coeffs(11);
a5l = coeffs(12);
a6l = coeffs(13);
a7l = coeffs(14);

% High Temperature Range polynomial coefficients.
a1h = coeffs(1);
a2h = coeffs(2);
a3h = coeffs(3);
a4h = coeffs(4);
a5h = coeffs(5);
a6h = coeffs(6);
a7h = coeffs(7);

% Calculate Cp or specific heat.
Cp_low = R*(a1l+a2l*T_l+a3l*T_l.^2+a4l*T_l.^3+a5l*T_l.^4);
Cp_high = R*(a1h + a2h*T_h + a3h*T_h.^2 + a4h*T_h.^3+a5h*T_h.^4);
Cp = [Cp_low, Cp_high];

% Calculate H or Enthalpy.
H_low = R*T_l.*(a1l + (a2l*T_l)/2 + (a3l*T_l.^2)/3 + (a4l*T_l.^3)/4 + (a5l*T_l.^4)/5 + a6l./T_l);
H_high = R*T_h.*(a1h + (a2h*T_h)/2 + (a3h*T_h.^2)/3 + (a4h*T_h.^3)/4 + (a5h*T_h.^4)/5 + a6h./T_h);
H = [H_low, H_high];

% Calculate S or Entropy
S_low = R*(a1l*log(T_l) + a2l*T_l + (a3l*T_l.^2)/2 + (a4l*T_l.^3)/3 + (a5l*T_l.^4)/4 + a7l);
S_high = R*(a1h*log(T_h) + a2h*T_h + (a3h*T_h.^2)/2 + (a4h*T_h.^3)/3 + (a5h*T_h.^4)/4 + a7h);
S = [S_low, S_high];

Code Explanation:

Line 2 is the function definition. This function calculates the thermodynamic properties Cp, H & S from the polynomial coefficients and temperature.

Outputs

Cp – Specific heat values array at specified temperature points.

H – Enthalpy values array

S – Entropy values array

Inputs

T_l – Low temperature range array (local)

T_h – High Temperature range array (local)

coeffs – the 14 polynomial coefficients for high and low temperature range.

Line 5 defines the universal gas constant R = 8.314

Line 8 to 14 creates variables for polynomial coefficients for the low temperature range from coeffs array.

Line 17 to 23 creates variables for polynomial coefficients for the high temperature range from coeffs array.

Line 26 & 27 calculates array of Cp values for low and high temperature range using the polynomial coefficients and temperature range values. Line 28 combines both the range of cp values into one array. Similarly, lines 31 to 33 & 36 to 38 calculate enthalpy and entropy respectively.

Code for main program:

clear all
close all
clc

%Open filw in read mode
f1 = fopen('THERMO.dat','r');

%Skip the header using fgetl
fgetl(f1);

%get the first line
gt = fgetl(f1);
%trim the leading and trailing spaces using srtrim()
gt=strtrim(gt);

%Get the global temperature range values
%Split the values with strsplit( ,' ')
gt=strsplit(gt,' ');

%Convert string to numbers
for n=1:3
    gt_num(n)=str2num(gt{n});
end
gt=gt_num;
gt_low=gt(1);
gt_mid=gt(2);
gt_high=gt(3);

%Skip the comments starting with !
tline='!';
position=0;
while tline(1)=='!'
    %Use ftell(fileID) to get the current cursor position
    before=ftell(f1);
    tline=fgetl(f1);
    if tline(1)~='!'
        position=fseek(f1,before,-1);
    end
end

%Store the initial position usinf ftell(fieID) to a variable
init_position=ftell(f1);
%Command Line Interface for getting the thermodynamic data
while true
    %Main menu
    disp('Main Menu');
    disp('----')
    disp('Enter 1 to save the thermodynamic plots as images for all the species');
    disp('Enter 2 to plot thermodynamic data for specific species');
    disp('Or just Press Enter to Exit');
    fprintf('nn');

    %Get user's choice
    species_choice = input('Enter your choice: ','s');
    
    %For choice 2 - Command line interfaceto get data for specific species
        if species_choice=="2"
            while true
                %Get the species name from the user
                fprintf('nn');
                fprintf('Enter the species name for which you want the thermodynamic data plot?');
                species=input('Species: ','s')
                fprintf('nn');
                fprintf('Thermodynamic plots for %sn', upper(species));
                disp('Enter 1 for Specific Heat Vs Temperature Plot');
                disp('Enter 2 for Enthalpy Vs Temperature Plot');
                disp('Enter 3 for Entropy Vs Temperature Plot');
                disp('Enter 4 for all of the above plots');
                % Get user's choice for thermodynamic properties to plot.
                plot_choice=input('Enter your choice:  ');
                plot_data(f1,species,init_position,plot_choice);
                fprintf('n');        
                choice=input('Enter Y to plot data for another species or nEnter M for Main Menu  or nJust Press Enter to exitn','s');
                if upper(choice)=='M'
                    clc
                    break;
                elseif choice==""
                    return;
                end
            end
    % For choice 1 - Get data for all species
        elseif species_choice=="1"
        plot_all_data(f1,init_position);
        fclose(f1);
        return;
		
    %For choice 3: Exit the program.
        elseif species_choice==""
        fclose(f1);
        return;
        end
end

Code explanation:

Lines 2 to 4: Clear command window & workspace and close all plots.

Line 7: Open the THERMO.dat file using fopen() command. Where first input ‘THERMO.dat’ is the name of the file and second input is mode which is ‘r’ (means open the file in read mode).

Line 10: Skip the first line (header) using fgetl() command. fgetl reads the first line and moves the pointer to the start of second line.

Line 13: Get the second line using fgetl in variable gt (short for global temperature).

Line 15: Trim the leading and trailing spaces in gt using strtrim() ffunction.

Line 19: Split one single string into an array by using strsplit(). First input is name of variable and second input is the character (‘ ’ in our case) to be used for splitting.

Line 22 to 24: Convert string to number using str2num() function.

Line 25 to 28: Create variables and assign values for global low, global mid and global high temperatures.

Lines 31 to 40: Skip comments starting with ‘!’. If the file THERMO.dat is revised and number of comment lines are increased in the revised version than this code can be used without any changes.

Line 31 & 32: Defines a variable tline = ‘!’ and position=0.

Line 33 defines the while loop. The loop goes on till the first character of the next line read from the file is not ‘!’.

Line 35 uses the ftell(fileID) command to get the current position of read cursor.

Line 36 gets the line as string from the THERMO.dat file using fget() function.

Line 37 checks if the first character is ! or not. If it is not ! then the cursor read position is changed in line 38 to before using fseek(fileID,offset,origin) function. Where f1 is the file ID, offset is the position to be offset(before variable in our case) from the origin (-1 means origin is start of the file, 0 is current position & 1 is end of file).

Line 39 is end of if statement and line 40 is end of while loop.

Line 43 stores the position of read cursor which is at the beginning of first species in THERMO.dat file to initial_position variable.

Line 46 define the while loop for the command line interface for plotting the thermodynamic data.

Lines 48 to 53: Prints instructions for user to decide whether to save the thermodynamic data for all the species as images or plot data in MATLAB only for some species and then save it as images.

Line 56: Prompts user to enter his choice.

Line 59, 86 & 92: The if and elseif statements based on the users’s choice, decide whether to plot data for all the species or only plot data for specific species or to exit the program.

Line 60: The while loop is for user to plot data for multiple species but not all the species. After plotting the properties for one species if the user wants to plot data for another species he can specify his choice.

Lines 62 to 64: Prompts the user to enter the name of species.

Lines 65 to 73: Prompts the user to enter his choice of thermodynamic properties to plot.

Line 74: Calls the user defined function plot_data() to plot the required data. Explanation for plot_data function is provided under section File 2 plot_data.m.

Line 75 to 76: Asks user whether he wants to plot data for another species or exit to main menu or exit the program. Based on his choice, if he wants exit the main menu, lines 77 to 79 breaks out of the inner while loop and if he wants to exit the program lines 80 to 82 returns the control of main program. And if he wants to plot data for another species, the while loop once again iterates.

Lines 86 to 89: In the main menu line 73 if user enters 1 then these statements are executed. Line 87 calls user defined function plot_all_data() to plot thermodynamic data for all the species. Line 88 closes the file THERMO.dat. Line 89 returns the control to MATLAB.

Line 91 to 94: If user just presses enter without any inputs, line 93 closes the file THERMO.dat and returns the control to MATLAB.

Line 95 is end of if-elseif statements and line 96 is end of while loop.

OBSERVATIONS

Below are the outputs for CO and the plots for other species can be found similarly.

Here are the all the folders showing all the files in them:

ERRORS FACED

1. Named the function for calculating molar mass as calculating_molar_mass() and called it by Molar_Mass() in the plot_data.m and plot_all_data.m functions. Corrected it by removing the making amends in the code.
2. A couple of other typos which were corrected as I read the code again.
3. Tried to change directory to D: drive when it was already the current directory. Corrected it by removing unnecessary statements.