Project 1 - Parsing NASA thermodynamic data

clear all
close all
clc
% opening file in matlab
f1=fopen('THERMO.dat','r')
% skipping 1st five lines
for i=1:5
    fgetl(f1)
end
 
for i=1:53
% reading & assigning 6th line 
header=fgetl(f1)
% spilting line based on space
a=strsplit(header,' ')
% assigning 1st species
species=a{1}
% finding low,mid & high temperatures
llt=str2double(a{end-3})
lmt=str2double(a{end-1})
lht=str2double(a{end-2})
 % reading 7th line from file
line1=fgetl(f1);
% finding the positions of E from line 7
e1=strfind(line1,'E');
% finding 1st coefficient
a1=str2double(line1(1:e1(1)+3));
% finding 2nd coefficient
a2=str2double(line1(e1(1)+4:e1(2)+3));
% finding 3rd coefficient
a3=str2double(line1(e1(2)+4:e1(3)+3));
% finding 4th coefficient
a4=str2double(line1(e1(3)+4:e1(4)+3));
% finding 5th coefficient
a5=str2double(line1(e1(4)+4:e1(5)+3));

% reading 8th line
line2=fgetl(f1);
% finding the positions of E from 8th line
e2=strfind(line2,'E');
% finding 6th coefficient
a6=str2double(line2(1:e2(1)+3));
% finding 7th coefficient
a7=str2double(line2(e2(1)+4:e2(2)+3));
% finding 8th coefficient
a8=str2double(line2(e2(2)+4:e2(3)+3));
% finding 9th coefficient
a9=str2double(line2(e2(3)+4:e2(4)+3));
% finding 10th coefficient
a10=str2double(line2(e2(4)+4:e2(5)+3));

% reading 9th line
 line3=fgetl(f1);
 %finding the positions of E from 9th line
 e3=strfind(line3,'E');
 % finding 11th coefficient
 a11=str2double(line3(1:e3(1)+3));
 % finding 12th coefficient
  a12=str2double(line3(e3(1)+4:e3(2)+3));
% finding 13th coefficient
a13=str2double(line3(e3(2)+4:e3(3)+3));
% finding 14th coefficient
a14=str2double(line3(e3(3)+4:e3(4)+3));
% giving temperature array
temp=linspace(llt,lht,500);
% universal gas constant
r=8.314;

% specific heat
cp=specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,lmt,temp,r);
% enthalphy
h=enthalphy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,temp,r,lmt);
% entropy
[s]=entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,lmt,r,temp);

% to plot the variation of cp with temp
figure(1)
plot(temp,cp,'linewidth',2,'color','g')(plotting graph for temperature vs cp)
xlabel('temperature in K')  (giving name for xaxis)
ylabel('specific heat in KJ/Kg.K') (giving name for y axis)
grid on
title(sprintf('variation of specific heat with temperature for %s',species))
(giving title for the graph)

% to plot the variation of h with temp
figure(2)
plot(temp,h,'linewidth',2,'color','r')(plotting graph for temperature vs h)
xlabel('temperature in K')
ylabel('enthalpy in KJ/Kg')
grid on
title(sprintf('variation of enthalpy with temperature for %s',species))

% to plot the variation of entropy with temp
figure(3)
plot(temp,s,'linewidth',2,'color','y')
xlabel('temperature in K')
ylabel('entropy in KJ/Kg')
grid on
title(sprintf('variation of entropy with temperature for %s',species))
% creating folder
mkdir(species)
% changing directory
cd(species)
% saving figures in png format
saveas(figure(1),'cp.png'); 
saveas(figure(2),'h.png');
saveas(figure(3),'s.png');
% changing to previous folder
cd ..
% to calculate the molecular weight
  [molecular_weight]=mw(species) 
end (ending for loop)


I wrote seperate functions for finding cp,h,s and molecular weight and that functions I am giving below:

Function for finding cp:

function[cp]=specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,lmt,temp,r)
(function initializing and giving input arguments)

if(temp<=lmt) (opening if loop with given condition)
    cp=r*(a8+(a9*temp)+(a10*(temp).^2)+(a11*(temp).^3)+(a12*(temp).^4))
(equation for finding cp for lower temperatures)
else
    cp=r*(a1+(a2*temp)+(a3*(temp).^2)+(a4*(temp).^3)+(a5*(temp).^4))
(equation for finding cp for higher temperatures)
end (ending if loop)
end  (ending for loop)

Function for finding h:

function [h]=enthalphy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,temp,r,lmt)
if(temp<=lmt)
    h=(a8+(a9*temp/2)+(a10*(temp).^2/3)+(a11*(temp).^3/4)+(a12*(temp).^4/5)+(a13./temp))*r.*temp;
(equation for finding h for lower temperatures)
else
    h=(a1+(a2*temp/2)+(a3*(temp).^2/3)+(a4*(temp).^3/4)+(a5*(temp).^4/5)+(a6./temp))*r.*temp;
(equation for finding h for higher temperatures)
end
end

Function for finding s:

function[s]=entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,lmt,r,temp)
if(temp<=lmt)
    s=r*(a8*(log(temp))+(a9*(temp))+(a10*(temp).^2/2)+( a11*(temp).^3/3)+(a12*(temp).^4/4)+a14);
(equation for finding s for lower temperatures)
else
    s=r*(a1*(log(temp))+(a2*(temp))+(a3*(temp).^2/2)+( a4*(temp).^3/3)+(a5*(temp).^4/4)+a7);
(equation for finding s for higher temperatures)
end
end

 Function for finding molecular_weight:

elements=['H','C','O','N','A','S']; (giving basic elements)
 atomic_weight=[1,12,16,14,39,32];    (giving atomic weight of each element)
 molecular_weight=0;                   (initializing moleculat wt as 0)  
 for i=1:length(species)
     for j=1:length(elements)
         if strcmp(species(i),elements(j))  (comparing strings)
             molecular_weight=molecular_weight+atomic_weight(j); 
(equation for calculating molecular weight)
             position=j;
         end
     end
     n=str2double(species(i)); (converting string into number)
     if(n>1)
         molecular_weight=molecular_weight+atomic_weight(position)*(n-1);
(calculating molecular_weight for n>1)
     end
 end
 disp(molecular_weight)(displaying molecular_weight)


I got some errors while doing this programe and I am attaching the screenshots of that errors along with this.I corrected that errors by taking proper steps.

I attaching the screenshots for saving the plots into different folders.
Also I am attaching that 3 plots for required species along with this.