Project 1 - Parsing NASA thermodynamic data

clear all
close all
clc

R = 8.314;
R1 = fopen('THERMO.dat','r');
fgetl(R1);
temp_range = fgetl(R1);
A=strsplit (temp_range,' ');
Global_min_temp=(A{2});
Global_min_temp=str2num(Global_min_temp);
Global_mid_temp=(A{3});
Global_mid_temp=str2num(Global_mid_temp);
Global_max_temp=(A{4});
Global_max_temp=str2num(Global_max_temp);

fgetl(R1);
fgetl(R1);
fgetl(R1);

for i=1:53
Spe=fgetl(R1);
M = strsplit(Spe,' ');
element = M{1};
low_temp = (M{end-3});
low_temp = str2num(low_temp);
mid_temp = (M{end-1});
mid_temp = str2num(mid_temp);
max_temp = (M{end-2});
max_temp = str2num(max_temp);

L1 = fgetl(R1);
L2 = fgetl(R1);
L3 = fgetl(R1);

p1 = findstr(L1,'E');
p2 = findstr(L2,'E');
p3 = findstr(L3,'E');

i1 = L1(1:p1(1)+3);
i1 = str2num(i1);

i2 = L1(p1(1)+4:p1(2)+3);
i2 = str2num(i2);

i3 = L1(p1(2)+4:p1(3)+3);
i3 = str2num(i3);

i4 = L1(p1(3)+4:p1(4)+3);
i4 = str2num(i4);

i5 = L1(p1(4)+4:p1(5)+3);
i5 = str2num(i5);


i6 = L2(1:p1(1)+3);
i6 = str2num(i6);

i7 = L2(p2(1)+4:p2(2)+3);
i7 = str2num(i7);

i8 = L2(p2(2)+4:p2(3)+3);
i8 = str2num(i8);

i9 = L2(p2(3)+4:p2(4)+3);
i9 = str2num(i9);

i10 = L2(p2(4)+4:p2(5)+3);
i10 = str2num(i10);


i11 = L3(1:p3(1)+3);
i11 = str2num(i11);

i12 = L3(p3(1)+4:p3(2)+3);
i12 = str2num(i12);

i13 = L3(p3(2)+4:p3(3)+3);
i13 = str2num(i13);

i14 = L3(p3(3)+4:p3(4)+3);
i14 = str2num(i14);

T = linspace(low_temp,max_temp,500);

% calculating H Cp S
% for high temp
% if(T>mid_temp)
% Cp = (i1 + i2*T + i3*T.^2 + i4*T.^3 + i5*T.^4)*R;
% H = (i1 + i2*(T/2) + i3*(T.^3)/2 + i4*(T.^3)/4 + i5*(T.^5)/4 + i6/T).*R.*T;
% S = (i1*log(T)+ i2*T + i3*(T.^2)/2 + i4*(T.^3)/2 + i5* (T.^4) /4 + i7)*R;
% else
% % for low temp
% Cp = (i8 + i9*T + i10*T.^2 + i11*T.^3 + i12*T.^4)*R;
% H = (i8+i9*(T/2)+i10*(T.^3)./2+i11*(T.^3)./4+i12*(T.^5)./4+i13/T).*R.*T;
% S = (i8*log(T)+i9*T+i10*(T.^2)/2+i11*(T.^3)/2+i12*(T.^4)/4+ i14)*R;
% end

for i=1:length(T)
     if T(i)>mid_temp
         cp(i)= (i1+i2.*T(i)+i3.*T(i).^2+i4.*T(i).^3+i5.*T(i).^4)*R;
         H(i)= R*T(i).*(i1+((i2*T(i)/2))+((i3*T(i).^2)/3)+((i4*T(i).^3)/4)+((i5*T(i).^4)/5)+((i6./T(i))));
         S(i)= (i1.*log(T(i))+i2.*T(i)+i3.*(T(i)^2)/2+i4.*(T(i)^3)/3+i5.*(T(i)^4)/4+i7)*R;
     else
         cp(i)= (i8+i9.*T(i)+i10.*T(i).^2+i11.*T(i).^3+i12.*T(i).^4)*R;
         H(i)= R*T(i).*(i8+((i9*T(i)/2))+((i10*T(i).^2)/3)+((i11*T(i).^3)/4)+((i12*T(i).^4)/5)+((i13./T(i))));
         S(i)= (i8.*log(T(i))+i9.*T(i)+i10.*(T(i)^2)/2+i11.*(T(i)^3)/3+i12.*(T(i)^4)/4+i14)*R;
     end
end

% Saving in a folder
 cd('C:\Users\HARSH ABDAR\Documents\MATLAB\Parsing')
 mkdir(element);
 cd(element);

 % To plot variation of cp with T
 figure(1)
 plot(T,cp, 'linewidth','2','Color','g')
 xlabel('Temperature in K')
 ylabel('Specific heat in KJ/kg.K')
 grid on
 title(sprintf('variation of specific heat with temperature for %s', species))
saveas(figure(1),'cp.png');

% to plot variation of enthalpy with temperature
figure(2)
plot(T,cp,'linewidth','2','colour','r')
xlabel('Temperature in K')
ylabel('Enthalpy in KJ/kg')
grid on
title(sprintf('variation of Enthalpy with temperature for %s', species))
saveas(figure(2),'H.png');

% to plot variation of entropy with temperature 
figure(3)
plot(T,S,'linewidth','2','Color','b')
xlabel('Temperature in K')
ylabel('Entropy in KJ/K')
grid on
title(sprintf('variation of Entropy with temperature for %s', species))
saveas(figure(3),'S.png');

cd('C:\Users\Harsh Abdar\Documents\MATLAB\Parsing')


 %Molecular weight calculation
 elements = {'H', 'C', 'O', 'N', 'A', 'S'};
 atomic_weight = [1,12,16,14,39,32];
 molecular_weight=0;
 for k=1:length(species)
     for l=1:length(elements)
         if strcmp(species(k),elements(1))
             molecular_weight = molecular_weight+atomic_weight(1);
             position=1
         end
     end
     n = str2double(species(k));
     if(n>1)
         molecular_weight = molecular_weight+(atomic_weight(position)*(n-1));
     end
 end
end