clear all
close all
clc
%step 1
f1 = fopen('THERMO.txt','r'); % f1 is used to open file using fopen command
header = fgetl(f1); % Header is used for getting header of the file
temp_range = fgetl(f1); % temp_range is to get temperature range
a = strsplit(temp_range,' '); % strsplit is used to split the temperature range in a seperate values
global_min_temp = str2double(a{2}); % str2double command is used to get the values and assaign to a charater 'a'
global_mid_temp = str2double(a{3}); % global_min_temp,global_mid_temp & global_max_temp are said to be temperature values
global_max_temp = str2double(a{4}); % global_temp is used to determine the temperature of the species
% Range of temperature
T = linspace(global_min_temp,global_max_temp,400); % T is specified as temperature and linspace command is used determine temperature range in graph
R = 8.314; % universal gas constant
% Skipping the next three comments in the file
unwanted1 = fgetl(f1); % unwanted is specified that the line not to be considered
unwanted2 = fgetl(f1); % so the next three line in the data file is not considered
unwanted3 = fgetl(f1); % so fgetl is use to get the lines and not considered
for i = 1:53 % for loop is used, i is said to be increament value
% Extracting line 1 & the name of the species
line1 = fgetl(f1); % line1 is used to get first line of the species
name = strtok(line1); % name is used to specifie the species, strtok command is used to determine as token
species_data(i).name = name; % species_data is used to save name of the species
% Extracting higher & lower temperature coeffs
line2 = fgetl(f1); % line2, 3 & 4 are used to determine line of the species data line
line3 = fgetl(f1); % fgetl is to get the line of species data
line4 = fgetl(f1); % f1 is the data file to be read
A = strfind(line2,'E'); % A, B & C are the values of the line2 ,3 & 4 are to be save
B = strfind(line3,'E'); % strfind is said to be finding the values or letter
C = strfind(line4,'E'); % E is said to be the finding the values of the command
% Extracting higher temperature coeffs
a1 = str2double(line2(1:A(1)+3)); % a1,2,3,4,5 is set of minmum temp calulation formula
a2 = str2double(line2(A(1)+4:A(2)+3)); % a6 & 7 are the middel temp calculation formula
a3 = str2double(line2(A(2)+4:A(3)+3)); % line2 is said to be the value of the species at minmum temp
a4 = str2double(line2(A(3)+4:A(4)+3)); % A is used to give the values from lines of species
a5 = str2double(line2(A(4)+4:A(5)+3));
a6 = str2double(line3(1:B(1)+3)); % line3 is said to be the value of the species at middel temp
a7 = str2double(line3(B(1)+4:B(2)+3)); % B is used to give the values from lines of species
% Extracting low temperature coeffs
b1 = str2double(line3(B(2)+4:B(3)+3)); % b1,2,3 is a set of middel temp calculation formula
b2 = str2double(line3(B(3)+4:B(4)+3)); % b4,5,6,7 are the maximum temp calculation formula
b3 = str2double(line3(B(4)+4:B(5)+3)); % line3 are to be the middel temp values of the species
b4 = str2double(line4(1:C(1)+3)); % B is mentioned above
b5 = str2double(line4(C(1)+4:C(2)+3)); % line4 is said to be the maximum temp values of the species
b6 = str2double(line4(C(2)+4:C(3)+3)); % C is to give the values from lines of species
b7 = str2double(line4(C(3)+4:C(4)+3));
% Specific heat calculation
cp = specificheat (a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_mid_temp) ;
% Enthalpy calculation
H = enthalpy (a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,global_mid_temp) ;
% Entropy calculation
S = entropy_calc(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_mid_temp) ;
% Creating individal folder for species & plotting
Plotting(cp,H,S,T,name) ;
% Molecular mass of individual species
Molecular_weight = molecular_weight(species_data(i).name);
end % for loop end here
fclose(f1); % closing f1
% Asking user whether to show plots of any particular element or not
element_plot_showcase(species_data) ;
%% function programs startes
% Calculation of specific heat
function cp = specificheat (a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_mid_temp)
% Specific heat function is started
% Comparing with global temp
% If less than global temp it will take min coeffs
if T > global_mid_temp
% T is grater than global mid temperature and if loop starts
cp = R.*(a1 + a2.*T + a3.*T.^2 + a4.*T.^3 + a5.*T.^4);
% R is executed in this formula to get the value for cp
% If temp is more it will take max coeffs
else
% If the value of the T is greater than global mid temp then else is executed
cp = R.*(b1 + b2.*T + b3.*T.^2 + b4.*T.^3 + b5.*T.^4);
% Then R is executed in this formula to get the value for cp
end
% If loop is ended and the specific heat values is taken
end
% Specific heat function is ended
% Calculation of enthalpy
function H = enthalpy (a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,global_mid_temp)
% Enthalpy function is starting
if T > global_mid_temp
% To calculate, T must be higher than global mid temp and if loop is starting
H = R.*(a1.*T + (a2.*T.^2)./2 + (a3.*T.^3)./3 + (a4.*T.^4)./4 + (a5.*T.^5)./5 +a6);
% R is executed in this formula to get the value for H
else
% Else is executed if the T is higher than global mid temp
H = R.*(b1.*T + (b2.*T.^2)./2 + (b3.*T.^3)./3 + (b4.*T.^4)./4 + (b5.*T.^5)./5 +b6);
% then R is executed in this formula to get the value for H
end
% If loop is ended and the enthalpy values is taken
end
% Enthalpy function is ended
% Calculation of entropy
function S = entropy_calc(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_mid_temp)
% Entropy function is starting
if T <= global_mid_temp
% If loop is starting and T is higher than global mid temp for calculation
S = R.*(a1.*log(T) + a2.*T + (a3.*T.^2)./2 + (a4.*T.^3)./3 + (a5.*T.^4)./4 +a7);
% R is executed in the formula to get the value for S
else
% Else is executed if the value is higher than global mid temp
S = R.*(b1.*log(T) + b2.*T + (b3.*T.^2)./2 + (b4.*T.^3)./3 + (b5.*T.^4)./4 +b7);
% than R is executed in this formula to get the value for S
end
% If loop is ended and the values of the entropy is taken
end
% Entropy function is ended
% Calculation molecular weight of species
function MW = molecular_weight(name) % MW is representing molecular weight and the function is starting
atoms = ['H','C','O','N','R']; % Atoms is used to representing atom molecules
fprintf('molecular weight of '); % Fprint is used to print a statement
fprintf(name); % Name of the species is printed
fprintf(' : '); % Printing : for the values to be printed after it
atomic_weight = [1 12 16 14 40]; % Atomic weight of the molecules is taken
MW = 0; % Executing when MW is equel to 0
for i = 1:length(name) % For loop is starting and i is used to increament
for j = 1:length(atoms) % Again one more for loop is starting and j is used to increament
if strcmp(name(i),atoms(j)) % if loop is started and strcmp is used to show name and atoms are exected
MW = MW + atomic_weight(j); % MW formula is used
position = j; % Position vector is used
end % End of if loop
end % End of for loop of j
n = str2double(name(i)); % Name of i is show in string n
if n > 1 % If loop is started and n is grater than 1
MW = MW + atomic_weight(position)*(n-1); % MW is executed and the
end % End of if loop of n
end % End of for loop of i
fprintf('%d',MW); % Printing MW values
disp(' '); % Disp is used to display
end % End of function of molecular weight
% Creating individal folder for each species & plotting
function Plotting(cp,H,S,T,name) % Plotting function starting
cd (pwd) % Cd is to creat folder and pwd is to returns the path of the folder
dir_name = name; % Dir name is to list the folder or file
mkdir(dir_name); % Mkdir is to vreify the file or folder is present
figure(1) % Figure 1 is created
plot (T,cp,'LineWidth',2,'Color','b') % Plot command is to create a graph
title (sprintf('Specific heat Vs Temperature for %s',name)) % Title is to give the graph a name
xlabel ('TEMPERATURE') % Xlabel is to label x axis of the graph
ylabel ('SPECIFIC HEAT') % Ylabel is to label y axis of the graph
figure(2) % Figure 2 is created
plot (T,H,'LineWidth',2,'Color','r') % Plot command is used to plot a graph
title (sprintf('Enthalpy Vs Temperature for %s',name)) % Title is to give the graph a name
xlabel ('TEMPERATURE') % Xlabel is to label x axis of the graph
ylabel ('ENTHALPY') % Ylabel is to label y axis of the graph
figure(3) % Figure 3 is created
plot (T,S,'LineWidth',2,'Color','g') % Plot command is used to plot a graph
title(sprintf('Entropy Vs Temperature for %s',name)) % Title is to give the graph a name
xlabel('TEMPERATURE') % Xlabel is to label x axis of the graph
ylabel('ENTROPY') % Ylabel is to label y axis of the graph
cd(dir_name) % Cd is to craete a folder and dir name is to list the path of folder
saveas(1,'cp_vs_t.png') % Saveas is to save the graph 1 in name of png extension
saveas(2,'h_vs_t.png') % Saveas is to save the graph 2 in name of png extension
saveas(3,'s_vs_t.png') % Saveas is to save the graph 3 in name of png extension
cd .. % Cd is to create a folder
end % End of the plotting function
% Element plot showcase
function [] = element_plot_showcase(species_data)
choice = questdlg('Do you want to view H-T,S-T,C-T plot of any element?','user respouse','YES','NO','NO');
ctrl = 0;
while ctrl == 0
switch choice
case 'YES'
h = msgbox({'The elements present in the database are as follows :' ...
'O, O2, H, H2, OH, H2O, HO2, H2O2, C, CH, CH2, CH2(S)' ...
'CH3, CH4, CO, CO2, HCO, CH2O, CH2OH, H3O, CH3OH, C2H' ...
'C2H2, C2H3, C2H4, C2H5, C2H6, CH2CO, HCCO, HCCOH, H2CN' ...
'HCN, HNO, N, NNH, N2O, NH, NH2, NH3, NO, NO2, HCNO, HOCN' ...
'HNCO, NCO, CN, HCNN, N2, AR, C3H8, C3H7, CH3CHO, CH2CHO'});
user = inputdlg('Enter the element name : ');
delete(h);
pos = 0;
for i = 1:53
temp_name = species_data(i).name;
if strcmpi(temp_name,user)==1
pos = i;
break
end
end
if pos ~= 0
path = '/Users/dhilipansivasamy/Desktop/MatLAB program/NASA thermodynamic';
filename1 = fullfile(path,temp_name,'cp_vs_t.png');
image1 = imread(filename1);
figure(1)
imshow(image1);
filename2 = fullfile(path,temp_name,'h_vs_t.png');
image2 = imread(filename2);
figure(2)
imshow(image2);
filename3 = fullfile(path,temp_name,'s_vs_t.png');
image3 = imread(filename3);
figure(3)
imshow(image3);
choice = questdlg('Do you want to view H-T,S-T,C-T plot of any other element?','user responce','YES','NO','NO');
else
uiwait(msgbox('please enter a valid element name . . . !' ))
end
case 'NO'
ctrl = 1;
end
end
end
Folders of the data stored
THERMO.txt