Project 1 - Parsing NASA thermodynamic data

AIM-NASA came up the polynoials that can be used to evaulate thermodynamic properties such as Cp,H and S.

formulae-:

Cp/R=a1+a2+a3*T^2+a4*T^3+a5*T^4

H/RT=a1+a2*(T/2)+a3*(T^2/3)+a4*(T^3/4)+a5*(T^4/5)+a6/T

S/R=a1 lnt +a2*T+a3*(T^2/2)+a4*(T^3/3)+a5*(T^4/4)+a7

R= universal gas constant

T= Temeprature

clear all
close all
 clc

 f1=fopen('THERMO.dat','r');
 
 for i=1:5
     A=fgetl(f1);
 end
 for j=1:53
 header=fgetl(f1);
 b=strsplit(header,'  ');
 species_name=b{1};
 local_low=str2double(b{end-3});
 local_high=str2double(b{end-2});
 local_mid=str2double(b{end-1});

  l1=fgetl(f1);

 e1=findstr(l1,'E');
 a1=str2double(l1(1:e1(1)+3));
 a2=str2double(l1(e1(1)+4:e1(2)+3));
 a3=str2double(l1(e1(2)+4:e1(3)+3));
 a4=str2double(l1(e1(3)+4:e1(4)+3));
 a5=str2double(l1(e1(4)+4:e1(5)+3));
 
   l2=fgetl(f1);
 e2=findstr(l2,'E');
 a6=str2double(l2(1:e2(1)+3));
 a7=str2double(l2(e2(1)+4:e2(2)+3));
 a8=str2double(l2(e2(2)+4:e2(3)+3));
 a9=str2double(l2(e2(3)+4:e2(4)+3));
 a10=str2double(l2(e2(4)+4:e2(5)+3));

 l3=fgetl(f1);
 e3=findstr(l3,'E');
 a11=str2double(l3(1:e3(1)+3));
 a12=str2double(l3(e3(1)+4:e3(2)+3));
 a13=str2double(l3(e3(2)+4:e3(3)+3));
 a14=str2double(l3(e3(3)+4:e3(4)+3));
 T=linspace(local_low,local_high,200);
 R=8.314;
for K=1:200
    if T(K)>=local_mid 
        cp(K)=R*(a1+a2*T(K)+a3*T(K)^2+a4*T(K)^3+a5*T(K)^4);
        H(K)=R*T(K)*(a1+a2*T(K)/2+(a3*T(K)^2)/3+(a4*T(K)^3)/4+(a5*T(K)^4)/5+a6/T(K));
        S(K)=R*(a1*log(T(K))+a2*T(K)+(a3*T(K)^2)/2+(a4*T(K)^3)/3+(a5*T(K)^4)/4+a7);
    else

        cp(K)=R*(a8+a9*T(K)+a10*T(K)^2+a11*T(K)^3+a12*T(K)^4);
        H(K)=R*T(K)*(a8+a9*T(K)/2+(a10*T(K)^2)/3+(a11*T(K)^3)/4+(a12*T(K)^4)/5+a13/T(K));
        S(K)=R*(a8*log(T(K))+a9*T(K)+(a10*T(K)^2)/2+(a11*T(K)^3)/3+(a12*T(K)^4)/4+a14);
    end
end

 figure(1)
 plot(T,cp)
 xlabel('temperature')
 ylabel('specific heat')
 title(species_name)

 figure(2)
 plot(T,H)
 xlabel('temperature')
 ylabel('enthalpy')
 title(species_name)
        
figure(3)
 plot(T,S)
 xlabel('temperature')
 ylabel('entropy')
 title(species_name) 
mkdir('C:\Users\Radhe\Documents\MATLAB',species_name)
cd(species_name)
saveas(figure(1),'specfic heat.png')
saveas(figure(2),'enthalpy.png')
saveas(figure(3),'entropy.png')
cd ..
M=molecular_weight(species_name)
 end



Molecular_Weight

function W= molecular_weight(species_name)
atomic_name=['H','C','O','N','Ar'];
atomic_weight=[1.00784,12.0107,15.999,14.0067,39.948];
W=0;
for i=1:length(species_name)
    for j=1:length(atomic_name)
        if species_name(i)==atomic_name(j)
            W=W+atomic_weight(j);
            pos=j;
        end
    end
    n=str2num(species_name(i));
    if n>=1
        W=W+atomic_weight(pos)*(n-1);
    end
end
fprintf('molecular weight of %s = %f',species_name,W)
disp(' ')
end