Parsing NASA Thermodynamic data using MatLab
Objective:
The main objective of this project is to write a code in matlab to prase the NASA Thermodynamic data file and to calculate the thermodynamic properties of various gas species.
NASA Thermodynamic data:
NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.
To create Function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file by reading this file NASA thermodynamic data.Also to calculate the molecular weight of each species.
File Parsing: (Parse = split)
It is the method of splitting the file or other input into a piece of data that can be easily stored or manipulated.By doing this we are going to split a string into parts then recognizing the parts to convert it into something simpler than a string.
Here in this project we are going to split the NASA thermodynamic data which is having 53 gas species values in 213 lines. To perform the calculation for thermodynamic properties, we need some sort of coefficient values and temperature ranges which is avaliable in the Nasa file.So from that we are going to parse the values which are all neede for our calculation.
Matlab program steps:
Sample NASA data:
Step:1 Open the file and read the data
Command used (fopen,fgetl,strsplit,str2double)
Here we are going to open the THEMO.dat file for reading, by using the command fopen() and assingnig to the variable called f1.
Then, the fgetl(f1) is the command acts like a line cursor, which will retrive the next time while running this command for second time and so on.
After reading the second line we have splitted the string with strsplit() command by giving the string name as GlobalTemp and the delimiter used for splitting is space(' '). Then it is assinged to the varible called temp. This will split the whole line into number of cells.
Then the string values are converted into number by using the command called str2double or str2num and assinged to the varible t. In this line t(2) indicates the index position. We can observe that in the second line we are having some unwanted space infront of 300.00, which is considered as the first value in that line. So 300.00 takes the 2nd position. This 300.00 string is converted as the number 300 and assinged to the variable called Global_low_Temp. Same command is followed for the next values.
The line with ! says that they are all comment lines which is not needed for calculation. so we can skip those by execting in the for loop command.
f1 = fopen('THERMO.dat','r'); % opening the file for reading
fgetl(f1); % Reading file Header
GlobalTemp = fgetl(f1); % Reading second line
temp = strsplit(GlobalTemp,' '); % Splitting string from the line
t = str2double(temp); % converting string into number
Global_low_Temp = t(2); % Assinging indices of 't' to a variable
Global_mid_Temp = t(3);
Global_high_Temp = t(4);
% Reading the comment lines with for loop
for i = 3:5
fgetl(f1);
endStep:2 Extracting the all species data
Command used(fgetl,strtok,strsplit)
In step 2 we are going to extract the species data, with the help of for loop and stored in the varibale called 'A'. This is will extract all the data in the single row. so we are converitn into colomn data with transpose operation.
After transposing, we are going to extract the row which is having only species names. For that species=A(1:4:end) is used, this will fectch the row with the differnce of 4, that means it will fetch 1st row and then 4th row and then 8th up to the end and stored in the varibale called species.
The extracted species data contains 53 rows, from that we are going to extract only the species names which is located in the first in all the rows. For this operation strtok command is used. This function parses str from left to right, using whitespace characters as delimiters, and returns part or all of the text in row.
Then from species data, local temperature is extracted with the strsplit command and stored as species_temp. Again this species_temp data is transposed into colomn and stored as tem_range.
These operations are done inside the for loop, so that we can extract all 53 species data.
% Extracting the all species data
for j = 1:213
A{j} = fgetl(f1);
end
A = A'; % Transpose operation from row to colomn
species = A(1:4:end); % Extracting line which is having only species name
% Extracting the species names from the cell
for k = 1:53 % strtok parses str from left to right with delimiters as whitespace
species_name = strtok(species);
species_temp{k} = strsplit(species{k},' ');
temp_range = species_temp';
endStep:3 Getting Input from the user
Here we are getting the input from the user end by prompting "Enter the valid species name ". This input is getting from the user in the form of string. Then that is converted into upper case.
After getting user input, then inside the for loop under if statement string is compared with the stored species_name value with the command called strcmp.
After comaprison the linecursor postion is saved as X.
% Input from the user
user_input = input('Enter valid species name: ','s');
species_input = upper(user_input);
% string comparison
for k = 1:53
if strcmp(species_input,species_name{k})
X = k;
end
end Step:4 Extracting the coefficient for all the gas species
Here we are going to extract the 14 coeffecient of gas species. The coefficients locations are like 5in 2ndrow , next 5 in 3rd row and rest 4 in the 4th row, which can be reffered in the above sample NASA data.
Now for extracting this we are going to locate the string 'E' first, which is present in each coefficient value. For doing that strfind command is used. 4*X-3+1 says the excat location of line from the user input to species location. From that line it fetch the string 'E' location value.
After finding the location of E, for the first coefficient we are saying that A{4*X-3+1}(1:line_1(1)+3). This means that form A,(contains 213 rows) {4*X-3+1} - gives the string E location cell value and (1:line_1(1)+3) gives the value from 1 to 12, (i.e upto string E) by adding +3 it will give the whole coefficient value.(line_1 = 12 27 42 57 72) these are all the location value of E.
After getting the first coefficient , the second coefficient value starts friom the end of first coefficient value. line_1(1) says the end location of string E position in the first coefficient. So line_1(1)+4 says the starting point of second coefficient. Then line_1(2) says the location of string E in the second coefficient and by adding +3 we will reach to the end point of second coefficient.
By doing this we can extract the 14 coefficient values. But remember the extracted coefficient values are still in the form of string only. For converting this to the number strdouble command is used inside the for loop for 14 coefficient.
After converting this into the number , then we are dividing the first 7(!:14) coefficient as high temperature coefficient and the next 7(8:14) is low temperature coefficient and stored in the varibale respectively.
% Extracting co-efficient of gas species
line_1 = strfind(A{4*X-3+1},'E');
a{1,1} = A{4*X-3+1}(1:line_1(1)+3);
a{2,1} = A{4*X-3+1}(line_1(1)+4:line_1(2)+3);
a{3,1} = A{4*X-3+1}(line_1(2)+4:line_1(3)+3);
a{4,1} = A{4*X-3+1}(line_1(3)+4:line_1(4)+3);
a{5,1} = A{4*X-3+1}(line_1(4)+4:line_1(5)+3);
line_2 = strfind(A{4*X-3+2},'E');
a{6,1} = A{4*X-3+2}(1:line_2(1)+3);
a{7,1} = A{4*X-3+2}(line_2(1)+4:line_2(2)+3);
a{8,1} = A{4*X-3+2}(line_2(2)+4:line_2(3)+3);
a{9,1} = A{4*X-3+2}(line_2(3)+4:line_2(4)+3);
a{10,1} = A{4*X-3+2}(line_2(4)+4:line_2(5)+3);
line_3 = strfind(A{4*X-3+3},'E');
a{11,1} = A{4*X-3+3}(1:line_3(1)+3);
a{12,1} = A{4*X-3+3}(line_3(1)+4:line_3(2)+3);
a{13,1} = A{4*X-3+3}(line_3(2)+4:line_3(3)+3);
a{14,1} = A{4*X-3+3}(line_3(3)+4:line_3(4)+3);
% conversion of stringvalue into numbers
for i = 1:14
coeff{i} = str2double(a{i});
end
% First 7 are High temperature co-efficients from (1 to 7)
% second 7 are Low temperature co-efficients form (8 to 14)
high_temp_coeff_range = coeff(1:7);
low_temp_coeff_range = coeff(8:14);Step:5 Extracting the local Temperature values
Each species is having their local temperature values. These temperature are different from species to species. For getting this first we are going to see the cell size for each species. It also differs from species to species, but it ranges from 1x8 to 1x12. So ifelse condition statement is used for extracting local temp value. The temp_range size is 53x1, by saying Temp= size(temp_range{X,1} gives the species row value from the user input. And for that, if the species is having the cell size with 8,9,10,11,or 12 as the colum value, then we are executing the condition statement if Temp(2) denotes species column size is equal to 8, then fetch the data for (1,5) which is having local temp value in the particuar species row. Then it stored as temp1 and next (1,7) and (1,6) is also having the local temp value.
After extracting all the local temperature values, then it is converted from string to number and sored as low_temp, mid_temp and high_temp.
Then linspace command is used from low_temp to mid_temp and stores as low_temp_range. The same is followed for high_temp_range but with mid_temp to high_temp as the range.
% Extracting local Temp ranges for each species
Temp = size(temp_range{X,1});
if Temp(2) == 8
temp1 = temp_range{X,1}(1,5);
temp2 = temp_range{X,1}(1,7);
temp3 = temp_range{X,1}(1,6);
else
if Temp(2) == 9
temp1 = temp_range{X,1}(1,6);
temp2 = temp_range{X,1}(1,8);
temp3 = temp_range{X,1}(1,7);
else
if Temp(2) == 10
temp1 = temp_range{X,1}(1,7);
temp2 = temp_range{X,1}(1,9);
temp3 = temp_range{X,1}(1,8);
else
if Temp(2) == 11
temp1 = temp_range{X,1}(1,8);
temp2 = temp_range{X,1}(1,10);
temp3 = temp_range{X,1}(1,9);
else
if Temp(2) == 12
temp1 = temp_range{X,1}(1,9);
temp2 = temp_range{X,1}(1,11);
temp3 = temp_range{X,1}(1,10);
else
end
end
end
end
end
% conversion of Local temp stringvalue into numbers
low_temp = str2double(temp1); % Low temperature
mid_temp = str2double(temp2); % Mid temperature
high_temp = str2double(temp3); % High temperature
% creating the range of low & high values of local temperatures
low_temp_range = linspace(low_temp,mid_temp,20);
high_temp_range = linspace(mid_temp,high_temp,20);Step:6 Function creation for calculation
For calculating the Specific heat, enthalpy and entropy function is used, by declaring the input paramters outside the function.Then the calculation is done for both high and low temperature coefficient values.
% Function Input for calculations
LT = low_temp_range;
HT = high_temp_range;
R = 8.314; % universal gas constant in j/mol.k
a = low_temp_coeff_range; % low Temp co-eff range from 8 to 14
b = high_temp_coeff_range; % High Temp co-eff range from 1 to 7
% Specific Heat Calculation
Cp_low = func_specific_heat_low(a,LT,R);
Cp_high = func_specific_heat_high(b,HT,R);
% Enthalpy Calculation
H_low = func_enthalpy_low(a,LT,R);
H_high = func_enthalpy_high(b,HT,R);
% Entropy Calculation
S_low = func_entropy_low(a,LT,R);
S_high = func_entropy_high(b,HT,R);Specific Heat calculation function:
function Cp_low = func_specific_heat_low(a,LT,R)
% Given Equ Cp/R = a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4
specific_heat_low = (a{1} + a{2}.*LT + a{3}.*LT.^2 + a{4}.*LT.^3 + a{5}.*LT.^4)*R;
Cp_low = specific_heat_low;
end
function Cp_high = func_specific_heat_high(b,HT,R)
specific_heat_high = (b{1} + b{2}.*HT + b{3}.*HT.^2 + b{4}.*HT.^3 + b{5}.*HT.^4)*R;
Cp_high = specific_heat_high;
endEnthalpy Function:
function H_low = func_enthalpy_low(a,LT,R)
% Given Equ H/RT = a1 + (a2*T)/2 + (a3*T^2)/3 + (a4*T^3)/4 + (a5*T^4)/5 + (a6)/T
enthalpy_low = (a{1} + (a{2}.*LT)/2 + (a{3}.*LT.^2)/3 + (a{4}.*LT.^3)/4 + (a{5}.*LT.^4)/5 + a{6}./LT).*R.*LT;
H_low = enthalpy_low;
end
function H_high = func_enthalpy_high(b,HT,R)
enthalpy_high = (b{1} + (b{2}.*HT)/2 + (b{3}.*HT.^2)/3 + (b{4}.*HT.^3)/4 + (b{5}.*HT.^4)/5 + b{6}./HT).*R.*HT;
H_high = enthalpy_high;
endEntropy Function:
function S_low = func_entropy_low(a,LT,R)
% GivenEqu S/R = a1*logT + a2*T + (a3*T^2)/2 + (a4*T^3)/3 + (a5*T^4)/4 + a7
entropy_low = (a{1}.*log(LT) + (a{2}.*LT) + (a{3}.*LT.^2)/2 + (a{4}.*LT.^3)/3 + (a{5}.*LT.^4)/4 + a{7})*R;
S_low = entropy_low;
end
function H_high = func_entropy_high(b,HT,R)
entropy_high = (b{1}.*log(HT) + (b{2}.*HT) + (b{3}.*HT.^2)/2 + (b{4}.*HT.^3)/3 + (b{5}.*HT.^4)/4 + b{7})*R;
H_high = entropy_high;
endStep:7
Molecular weight calculation:
The speceis elements are declared in the form of cell array. Their atomic mass is decalred in the form of array. Then initialy moleuclat weight is declared as zero. Then in the for loop for the species_input which is the user input string comparison is done for species_input and for the species_elements.
Then in the formula molecular weight and the atomic mass is added for that particular species.Then that position is stored as y.
AFter the nested for loop the species_input which is in the form of string is converted as number with strdouble command. Then the species_input which is now as number is greater than 1, molecular weight is again calculated by having the atomic mass for the current species * n-1.
% Molecular weight calculation
species_elements = {'O','H','C','N','S','A'}; % Declaring species elements
atomic_mass = [16,1,12,14,32,40]; % Declaring the atomic mass value
mol_wt = 0; % declaring molecular weight is '0' in initial stage
for i = 1:length(species_input)
for u = 1:length(species_elements)
if strcmp(species_input(i),species_elements(u))
mol_wt = mol_wt + atomic_mass(u);
y = u;% position is stored
end
end
n = str2double(species_input(i));
if n>1
mol_wt = mol_wt + atomic_mass(y)*(n-1);
end
endStep:8
Plotting code: Plotting is done for Specific Heat and Temperature , Enthalpy vs Temperature and Entropy vs Tempertaure for both low and high range.
% Plotting specific heat vs temp
figure(1)
plot(low_temp_range,Cp_low,'b','linewidth',2)
hold on
plot(high_temp_range,Cp_high,'r','linewidth',2)
xlabel('Temperature (K)');
ylabel('Specific Heat (J/K Kg)');
legend('Low temp','High Temp');
title({'Specific Heat Vs Temperature',['Species name = ',species_input,' & Molecular weight = ',num2str(mol_wt)]});
% plotting Enthalpy vs Temp
figure(2)
plot(low_temp_range,H_low,'g','linewidth',2)
hold on
plot(high_temp_range,H_high,'m','linewidth',2)
xlabel('Temperature (K)');
ylabel('Enthalpy (J/K)');
legend('Low temp','High Temp');
title({'Enthalpy Vs Temperature',['Species name = ',species_input,' & Molecular weight = ',num2str(mol_wt)]});
% plotting Entropy vs Temp
figure(3)
plot(low_temp_range,S_low,'r','linewidth',2)
hold on
plot(high_temp_range,S_high,'m','linewidth',2)
xlabel('Temperature (K) rightarrow');
ylabel('Entropy (J/K) rightarrow');
legend('Low temp','High Temp');
title({'Entropy Vs Temperature',['Species name = ',species_input,' & Molecular weight = ',num2str(mol_wt)]});Step:9
Saving the plot and Molecular number: The plots are saved in the seperate directory folder by using matlab command. The firstline mkdir says create a new directory or folder with a name called Fileparsing and the specific folder name inside fileparsing is taken from the species_input.
cd says current directory, with this plots are saved inside the folder name called fileparsing with the specific species input as the folder name. so each species folder saves 3 plots inside it.
Then in the same directory a new .txt file is open in the name of molecular weight and the inputs which is given from the user end i.e species input will apppend into the txt file in the form of string.
fwrite is used to write the input and sprintf is used for writing the value as string. Then format specifier %s is used to specifiy the format and rn is used as enter which write the input in the next row.
% Saving the plotted file in the system directory with the folder name
mkdir('Fileparsing',species_input)
cd(['D:MatlabFileparsing',species_input])
saveas(figure(1),'Cp vs Temp.png')
saveas(figure(2),'Enthalpy vs Temp.png')
saveas(figure(3),'Entropy vs Temp.png')
% Writing the molecular weight in a .txt file
cd('D:MatlabFileparsing');
fileID = fopen('molecular weight.txt','a');
fwrite(fileID,sprintf('%s %6srn','speceies_name',' mol_wt'));
fwrite(fileID,sprintf('%s %20.2frn',species_input,mol_wt));
fclose(fileID);Full program:
close all
clear all
clc
f1 = fopen('THERMO.dat','r'); % opening the file for reading
fgetl(f1); % Reading file Header
GlobalTemp = fgetl(f1); % Reading second line
temp = strsplit(GlobalTemp,' '); % Splitting string from the line
t = str2double(temp); % converting string into number
Global_low_Temp = t(2); % Assinging indices of 't' to a variable
Global_mid_Temp = t(3);
Global_high_Temp = t(4);
% Reading the comment lines with for loop
for i = 3:5
fgetl(f1);
end
% Extracting the all species data
for j = 1:213
A{j} = fgetl(f1);
end
A = A'; % Transpose operation from row to colomn
species = A(1:4:end); % Extracting line which is having only species name
% Extracting the species names from the cell
for k = 1:53 % strtok parses str from left to right with delimiters as whitespace
species_name = strtok(species);
species_temp{k} = strsplit(species{k},' ');
temp_range = species_temp';
end
% Input from the user
user_input = input('Enter valid species name: ','s');
species_input = upper(user_input);
% string comparison
for k = 1:53
if strcmp(species_input,species_name{k})
X = k;
end
end
% Extracting co-efficient of gas species
line_1 = strfind(A{4*X-3+1},'E');
a{1,1} = A{4*X-3+1}(1:line_1(1)+3);
a{2,1} = A{4*X-3+1}(line_1(1)+4:line_1(2)+3);
a{3,1} = A{4*X-3+1}(line_1(2)+4:line_1(3)+3);
a{4,1} = A{4*X-3+1}(line_1(3)+4:line_1(4)+3);
a{5,1} = A{4*X-3+1}(line_1(4)+4:line_1(5)+3);
line_2 = strfind(A{4*X-3+2},'E');
a{6,1} = A{4*X-3+2}(1:line_2(1)+3);
a{7,1} = A{4*X-3+2}(line_2(1)+4:line_2(2)+3);
a{8,1} = A{4*X-3+2}(line_2(2)+4:line_2(3)+3);
a{9,1} = A{4*X-3+2}(line_2(3)+4:line_2(4)+3);
a{10,1} = A{4*X-3+2}(line_2(4)+4:line_2(5)+3);
line_3 = strfind(A{4*X-3+3},'E');
a{11,1} = A{4*X-3+3}(1:line_3(1)+3);
a{12,1} = A{4*X-3+3}(line_3(1)+4:line_3(2)+3);
a{13,1} = A{4*X-3+3}(line_3(2)+4:line_3(3)+3);
a{14,1} = A{4*X-3+3}(line_3(3)+4:line_3(4)+3);
% conversion of stringvalue into numbers
for i = 1:14
coeff{i} = str2double(a{i});
end
% First 7 are High temperature co-efficients from (1 to 7)
% second 7 are Low temperature co-efficients form (8 to 14)
high_temp_coeff_range = coeff(1:7);
low_temp_coeff_range = coeff(8:14);
% Extracting local Temp ranges for each species
Temp = size(temp_range{X,1});
if Temp(2) == 8
temp1 = temp_range{X,1}(1,5);
temp2 = temp_range{X,1}(1,7);
temp3 = temp_range{X,1}(1,6);
else
if Temp(2) == 9
temp1 = temp_range{X,1}(1,6);
temp2 = temp_range{X,1}(1,8);
temp3 = temp_range{X,1}(1,7);
else
if Temp(2) == 10
temp1 = temp_range{X,1}(1,7);
temp2 = temp_range{X,1}(1,9);
temp3 = temp_range{X,1}(1,8);
else
if Temp(2) == 11
temp1 = temp_range{X,1}(1,8);
temp2 = temp_range{X,1}(1,10);
temp3 = temp_range{X,1}(1,9);
else
if Temp(2) == 12
temp1 = temp_range{X,1}(1,9);
temp2 = temp_range{X,1}(1,11);
temp3 = temp_range{X,1}(1,10);
else
end
end
end
end
end
% conversion of Local temp stringvalue into numbers
low_temp = str2double(temp1); % Low temperature
mid_temp = str2double(temp2); % Mid temperature
high_temp = str2double(temp3); % High temperature
% creating the range of low & high values of local temperatures
low_temp_range = linspace(low_temp,mid_temp,20);
high_temp_range = linspace(mid_temp,high_temp,20);
% Function Input for calculations
LT = low_temp_range;
HT = high_temp_range;
R = 8.314; % universal gas constant in j/mol.k
a = low_temp_coeff_range; % low Temp co-eff range from 8 to 14
b = high_temp_coeff_range; % High Temp co-eff range from 1 to 7
% Specific Heat Calculation
Cp_low = func_specific_heat_low(a,LT,R);
Cp_high = func_specific_heat_high(b,HT,R);
% Enthalpy Calculation
H_low = func_enthalpy_low(a,LT,R);
H_high = func_enthalpy_high(b,HT,R);
% Entropy Calculation
S_low = func_entropy_low(a,LT,R);
S_high = func_entropy_high(b,HT,R);
% Molecular weight calculation
species_elements = {'O','H','C','N','S','A'}; % Declaring species elements
atomic_mass = [16,1,12,14,32,40]; % Declaring the atomic mass value
mol_wt = 0; % declaring molecular weight is '0' in initial stage
for i = 1:length(species_input)
for u = 1:length(species_elements)
if strcmp(species_input(i),species_elements(u))
mol_wt = mol_wt + atomic_mass(u);
y = u;% position is stored
end
end
n = str2double(species_input(i));
if n>1
mol_wt = mol_wt + atomic_mass(y)*(n-1);
end
end
% Plotting specific heat vs temp
figure(1)
plot(low_temp_range,Cp_low,'b','linewidth',2)
hold on
plot(high_temp_range,Cp_high,'r','linewidth',2)
xlabel('Temperature (K)');
ylabel('Specific Heat (J/K Kg)');
legend('Low temp','High Temp');
title({'Specific Heat Vs Temperature',['Species name = ',species_input,' & Molecular weight = ',num2str(mol_wt)]});
% plotting Enthalpy vs Temp
figure(2)
plot(low_temp_range,H_low,'g','linewidth',2)
hold on
plot(high_temp_range,H_high,'m','linewidth',2)
xlabel('Temperature (K)');
ylabel('Enthalpy (J/K)');
legend('Low temp','High Temp');
title({'Enthalpy Vs Temperature',['Species name = ',species_input,' & Molecular weight = ',num2str(mol_wt)]});
% plotting Entropy vs Temp
figure(3)
plot(low_temp_range,S_low,'r','linewidth',2)
hold on
plot(high_temp_range,S_high,'m','linewidth',2)
xlabel('Temperature (K) rightarrow');
ylabel('Entropy (J/K) rightarrow');
legend('Low temp','High Temp');
title({'Entropy Vs Temperature',['Species name = ',species_input,' & Molecular weight = ',num2str(mol_wt)]});
% Saving the plotted file in the system directory with the folder name
mkdir('Fileparsing',species_input)
cd(['D:MatlabFileparsing',species_input])
saveas(figure(1),'Cp vs Temp.png')
saveas(figure(2),'Enthalpy vs Temp.png')
saveas(figure(3),'Entropy vs Temp.png')
% Writing the molecular weight in a .txt file
cd('D:MatlabFileparsing');
fileID = fopen('molecular weight.txt','a');
fwrite(fileID,sprintf('%s %6srn','speceies_name',' mol_wt'));
fwrite(fileID,sprintf('%s %20.2frn',species_input,mol_wt));
fclose(fileID);Plotting:
O2 plot:
Plottting is done for cp_low vs low_temp_range and cp_high vs high_temp_range and ploted as single figure. The legend shows the differnce between high and low tempaertaure.
Species name and Molecular weight for that particular species is also printted below the title of the plot. The same procedure is followed for the below plots with enthalpy and entropy.
N2 Plot:
Co2 Plot:
Folder saving the plots:
Molecular weight txt file:
Conclusion:
Thus File parsing is done for NASA Thermodynamic Data and the thermodynamic properties like specific heat , enthalpy and entropy are calculated. Also plotting is done for the same and stored in the seperate directory by using matlab code.
Refrences:-
https://en.wikipedia.org/wiki/Parsing
https://www.mathworks.com/matlabcentral/answers/317006-create-new-folders-by-mkdir
https://www.mathworks.com/matlabcentral/answers/17580-how-to-define-a-path-in-saveas-command
https://www.mathworks.com/matlabcentral/answers/106320-make-a-loop-to-calculate-molecular-weight