Project 1 - Parsing NASA thermodynamic data

%Project 1 nasa thermodynamic
close all
clear all
clc
R = 8.314; % J/mol/K (Universal gas constant)
data = fopen('THERMO.dat','r');  %read the file
for i = 1:2
    global_tem_range = fgetl(data); % Read global temperature
end
global_tem = str2num(global_tem_range); %converting character to num
low_tem = global_tem(1);
mid_tem = global_tem(2);
high_tem = global_tem(3);

% By pass comment line
for i = 1:3
    useless_comment = fgetl(data);
end

% For loop to find value of coefficient
for i = 1:53
    
    spe = fgetl(data); % Read the line
    name_str = findstr(spe,' ');
    species_name = spe(1:name_str(1)-1); % Read the species name
    t = findstr(spe,'.'); % Find the position of '.' in the line, so that we can calculate local temperature
    t1 = str2double(spe(t(1)-3:t(1)+3)); % local low temperature
    t2 = str2double(spe(t(2)-4:t(2)+3)); % local high temperature
    t3 = str2double(spe(t(3)-4:t(3)+3)); % local mid temperature
    T = linspace(t1,t2,500); %(local low temperature,local high temperature,no. of values)
    
    l1 = fgetl(data); % Read line one of coefficient
    l1e = findstr(l1,'E'); % find E in the line
    
    l2 = fgetl(data); % Read line two of coefficient
    l2e = findstr(l2,'E');
    
    l3 = fgetl(data); % Read line three of coefficient
    l3e = findstr(l3,'E');
    
    a1 = str2double( l1(1:l1e(1)+3)); % defining the coefficient
    a2 = str2double(l1(l1e(1)+4:l1e(2)+3));
    a3 = str2double(l1(l1e(2)+4:l1e(3)+3));
    a4 = str2double(l1(l1e(3)+4:l1e(4)+3));
    a5 = str2double(l1(l1e(4)+4:l1e(5)+3));
    
    a6 = str2double(l2(1:l2e(1)+3));
    a7 = str2double(l2(l2e(1)+4:l2e(2)+3));
    a8 = str2double(l2(l2e(2)+4:l2e(3)+3));
    a9 = str2double(l2(l2e(3)+4:l2e(4)+3));
    a10 = str2double(l2(l2e(4)+4:l2e(5)+3));
    
    a11 = str2double(l3(1:l3e(1)+3));
    a12 = str2double(l3(l3e(1)+4:l3e(2)+3));
    a13 = str2double(l3(l3e(2)+4:l3e(3)+3));
    a14 = str2double(l3(l3e(3)+4:l3e(4)+3));
    
    [c_p,H,S] = thermo_property(T,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,mid_tem); % call the function and find specific heat, enthalpy and entropy
    [m] = molecular_mass(species_name); % call the function and calculate molecular mass and display into command window
    
    %create directory
    mkdir(species_name)
    cd(sprintf(species_name))
    
    figure(1) %ploting the graph
    plot(T,c_p,'linewidth',2,'color','g')
    grid on
    xlabel('Temperature')
    ylabel('Specific_heat')
    title('Species name=',species_name)
    
    figure(2)
    plot(T,H,'linewidth',2,'color','r')
    grid on
    xlabel('Temperature')
    ylabel('Enthalpy')
    title('Species_name',species_name)
    
    figure(3)
    plot(T,S,'linewidth',2,'color','b')
    grid on
    xlabel('Temperature')
    ylabel('Entropy')
    title('Species_name',species_name)
    %Save the plot image
    saveas(figure(1),'Specific_heat.jpg')
    saveas(figure(2),'Enthalpy.jpg')
    saveas(figure(3),'Entropy.jpg')
    
    change_folder = cd('D:\Matlab\Exam\PROJECT_ONE_THERMODYNAMIC') %change the folder where function is kept
    
end

%Function to find molecular weight
function [m] = molecular_mass(species_name)
atom = ['H';'C';'N';'O';'A';'S']; %Atom name
atomic_weight = [1 12 14 16 40 32]; %atomic weight
m = 0;
%loop to find molecular weight of each species
for i = 1:length(species_name)
    for j = 1:length(atom)
        if strcmp(species_name(i),atom(j))
            m = m + atomic_weight(j);
            position = j;
        end
    end
    n = str2double(species_name(1));
    if n>1
        m = m + atomic_weight(position)*(n-1); %Molecular weight
    end
end
%Display molecular weight into the command window
sprintf('Molecular_weight_of %s : ',species_name)
sprintf('%f',m)
disp(' ') 
end

%Find specific heat, Enthalpy,Entropy
function [c_p,H,S] = thermo_property(T,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,mid_tem)

if T > mid_tem % according to problem first 7 coefficient are high temperature coefficient and second 7 are low temperature coefficient
    c_p = (a1 + (a2.*T) + (a3*T.^2) + (a4*T.^3) + (a5*T.^4)).*R;
    H = (a1 + (a2*T./2) + (a3*(T.^2)/3) + (a4*(T.^3)/4) + (a5*(T.^4)/5) + (a6./T)).*(R.*T);
    S = (a1.*log(T) + a2.*T +a3.*(T.^2)/3 + a4.*(T.^3)/3 + a5.*(T.^4)/4 +a7).*R;
else 
    c_p = (a8 + (a9*T) + (a10*T.^2) + (a11*T.^3) + (a12*T.^4)).*R;
    H = (a8 + (a9*T./2) + (a10*(T.^2)/3) + (a11*(T.^3)/4) + (a12*(T.^4)/5) + (a13./T)).*(R.*T);
    S = (a8.*log(T) + a9.*T +a10.*(T.^2)/3 + a11.*(T.^3)/3 + a12.*(T.^4)/4 +a14).*R;
end
end