Parsing NASA Thermodynamics data - MATLAB

Aim: To write a code in octave to parse the NASA thermodynamics data file and then calculate the thermodynamic properties.

Introduction:

Thermodynamic properties are the characteristics of a system through which it can be specified, there are two types of properties;

1) Extensive Thermodynamic Properties

2) Intensive Thermodynamic Properties

Code:

clear all 
close all
clc

R=8.314;


data=fopen('THERMO.dat','r');


fgetl(data);


global_temp_range=fgetl(data);
a=strsplit(global_temp_range,'.');
global_low_temp=str2num(a{1});
global_mid_temp=str2num(a{2})(2);
global_high_temp=str2num(a{3})(2);

temp=linspace(global_low_temp,global_high_temp,400);

for i = 1:3
    fgetl(data);
    
end      

for i = 1:53
    line1=fgetl(data)
    A=strsplit(line1,'.');
    speciesN= (a{1}(1));
    
    t2=fgetl(data);
    t3=fgetl(data);
    t4=fgetl(data);
    
    j=strfind(t2,'E');
    k=strfind(t3,'E');
    l=strfind(t4,'E');
    
    
    h1=str2num(t2(1:j(1)+3));
    h2=str2num(t2(j(1)+4:j(2)+3));
    h3=str2num(t2(j(2)+4:j(3)+3));
    h4=str2num(t2(j(3)+4:j(4)+3));
    h5=str2num(t2(j(4)+4:j(5)+3));
    h6=str2num(t3(1:k(1)+3));
    h7=str2num(t3(k(1)+4:k(2)+3));
    
    
    n1=str2num(t3(k(2)+4:k(3)+3));
    n2=str2num(t3(k(3)+4:j(4)+3));
    n3=str2num(t3(k(4)+4:j(5)+3));
    n4=str2num(t4(1:l(1)+3));
    n5=str2num(t4(l(1)+4:l(2)+3));
    n6=str2num(t4(l(2)+4:l(3)+3));
    n7=str2num(t4(l(3)+4:l(4)+3));
    
    [c_p]=Specific_heat(h1,h2,h3,h4,h5,n1,n2,n3,n4,n5,R,global_mid_temp);
    [h_rt]=enthalpy(h1,h2,h3,h4,h5,h6,n1,n2,n3,n4,n5,n6,R,global_mid_temp);
    [S]=enthalpy(h1,h2,h3,h4,h5,h6,h7,n1,n2,n3,n4,n5,n6,n7,R,global_mid_temp);
    
  mkdir(speciesN);
  current_dir=pwd;
  cd(speciesN)

  figure(1)

  plot(temp,c_p,'linewidth',2,'color','g');
  title(speciesN)
  xlabel("Temperature")
  ylabel("Specific Heat")
  test=sprintf('frame_cp%d.png',i);
  saveas(gca,test);

  figure(2)

  plot(temp,h_rt,'linewidth',2,'color','b');
  title(speciesN)
  xlabel("Temperature")
  label("Enthalpy")
  test=sprintf('frame_h%d.png',i);
  saveas(gca,test);

  figure(3)

  plot(temp,S,'linewidth',2,'color','b');
  title(speciesN)
  xlabel("Temperature")
  ylabel("Entropy")
  test=sprintf('frame_s%d.png',i);
  saveas(gca,test);

  mw=molecular_weight(speciesN)
  cd(current_dir);

end

Specific heat function:

function[c_p]=Specific_heat(h1,h2,h3,h4,h5,n1,n2,n3,n4,n5,R,global_mid_temp);

  for i=1:400
    if temp(i)<=global_mid_temp

    c_p(i)=R*(n1+(n2*temp(i))+(n3*temp(i)).^2)+(n4*(temp(i)).^3)+(n5*(temp(i)).^4));

    else

    c_p(i)=R*(h1+(h2*temp(i))+(h3*temp(i)).^2)+(h4*(temp(i)).^3)+(h5*(temp(i)).^4));
    end
end  
  

Similarly the programs are written for entropy and enthalpy.

PLOTS:

1) Specific Heat vs Temperature:

2) Entropy vs Temp:

3)Enthalpy vs Temp: