Parsing and Calculating the NASA thermodynamic data

Input Data

• Atomic Weight:

H= 1.007

C= 12.0096

O= 15.99903

N= 14.00643

S= 32.059

A= 40.078

• Universal Gas Constant:

R= 8.314 J/K-mol

Formulae

• `c_p/R = C1+C2*T+C3*T^2+C4*T^3+C5`
• `H/(RT)= C1+C2*T+C3*T^2/3+C4*T^3/4+C5*T^4/5+(C6)/T`
• `S/R= C1*log(T)+C2*T+C3*T^2/2+C4*T^3/3+C5*T^4/4+C7`

Commands

• fileID = fopen(filename) opens the file, filename, for binary read access, and returns an integer file identifier equal to or greater than 3.
• tline = fgetl(fileID) returns the next line of the specified file, removing the newline characters.
• mkdir folderName creates the folder folderName.
• cd displays the current folder.
• token = strtok(str,delimiters) parses str using the characters in delimiters. If delimiters includes more than one character, then strtok treats each character in delimiters as a separate delimiter. Because the delimiters are individual characters, delimiters can be any size, and the characters within delimiters can be in any order.
• X = str2num(chr) converts a character array or string scalar to a numeric matrix. The input can include spaces, commas, and semicolons to indicate separate elements. If str2num cannot parse the input as numeric values, then it returns an empty matrix.
• k = strfind(str,pattern) searches str for occurrences of pattern. The output, k, indicates the starting index of each occurrence of pattern in str. If pattern is not found, then strfind returns an empty array, []. The strfind function executes a case-sensitive search.
• C = strsplit(str,delimiter) splits str at the delimiters specified by delimiter.If str has consecutive delimiters, with no other characters between them, then strsplit treats them as one delimiter. For example, both strsplit('Hello,world',',') and strsplit('Hello,,,world',',') return the same output.

CODE

• Main Code

clear all
clc
close all 

%Pre-processing Data
thermo_data=fopen('thermo.dat','r');
header = fgetl(thermo_data);

%Global Temperatures
l1= fgetl(thermo_data);
l2=strsplit(l1,' ');
global_low_temp=str2num(l2{2});
global_mid_temp=str2num(l2{3});
global_high_temp=str2num(l2{4});

%Omitting Extra Data
for i=1:3
    z=fgetl(thermo_data);
end

%Creating New-Folder
mkdir Compounds
A2=fgetl(thermo_data)

%Parsing the Species
for i=1:53
 
    species_name=strtok(A2) 
    species= char(species_name);
    p=strfind(A2,'.');
    q=strsplit(A2(p-3:end-1),' ');
    
    %Local Temperatures
    T_low=str2num(q{1});
    T_high=str2num(q{2});
    T_mid=str2num(q{3});
    T=linspace(T_low,T_high,100);
   
    %Reading Coefficient Data
    X=fgetl(thermo_data);
    Y=fgetl(thermo_data);
    Z=fgetl(thermo_data);
    
    a=strfind(X,'E');
    b=strfind(Y,'E');
    c=strfind(Z,'E');
    
    %Assigning Variables to the Coefficient Data
    %High temperature Coefficient
    C1=str2num(X(1:a(1)+3));
    C2=str2num(X(a(1)+4:a(2)+3));
    C3=str2num(X(a(2)+4:a(3)+3));
    C4=str2num(X(a(3)+4:a(4)+3));
    C5=str2num(X(a(4)+4:a(5)+3));
    C6=str2num(Y(1:b(1)+3));
    C7=str2num(Y(b(1)+4:b(2)+3));
    
    %Low temperature Coefficient
    C8=str2num(Y(b(2)+4:b(3)+3));
    C9=str2num(Y(b(3)+4:b(4)+3));
    C10=str2num(Y(b(4)+4:b(5)+3));
    C11=str2num(Z(1:c(1)+3));
    C12=str2num(Z(c(1)+4:c(2)+3));
    C13=str2num(Z(c(2)+4:c(3)+3));
    C14=str2num(Z(c(3)+4:c(4)+3));
    
    
    %Calculating the Specific heat,Enthalpy, Entropy,Molecular mass
    [c_p, H, S] = temp_coeff(C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,T_low,T_high,T_mid);
     molecular_mass =molecular_weight(species)
    
     
    %Display the current folder
    cd Compounds
   
    %Creating and Assigning the folder to the Species
    mkdir(species);
    cd (species)
    
    %Plotting and Saving the Cp vs Temp. Graph
    figure(1)
    plot (c_p,T,'linewidth',2)
    xlabel('Specific Heat J/K-mol')
    ylabel('Temperature K')
    title ('Specific Heat vs Temperature')
    saveas(1,'c_p vs temp.fig')
    saveas(1,'c_p vs temp.jpeg')
    
    %Plotting and Saving the H vs Temp. Graph
    figure(2)
    plot (H,T,'linewidth',2)
    xlabel('Enthalpy J/mol')
    ylabel('Temperature K')
    title ('Enthalpy vs Temperature')
    saveas(2,'H vs temp.fig')
    saveas(2,'H vs temp.jpeg')
    
    %Plotting and Saving the S vs Temp. Graph
    figure(3)
    plot (S,T,'linewidth',2)
    xlabel('Entropy J/K-mol')
    ylabel('Temperature K')
    title ('Entropy vs Temperature')
    saveas(3,'S vs temp.fig')
    saveas(3,'S vs temp.jpeg')
    
    %Saving the species properties in to a tex file
    f=fopen('Properties.txt','w');
    fprintf(f, '%16s : %16s \n','Species Name',species);
    fprintf(f, '%16s : %16.4s \n','Molecular Weight',molecular_mass);
    fprintf(f,'Local Temperature Range: %4.1fk to %4.1fk \n\n',T_low,T_high);
    fprintf(f, '%5s %5s %15s %15s \n', 'Temperature', 'Specific Heat J/K-mol ', 'Enthalpy J/mol','Entropy J/K-mol');
   for k= 1: length(T)  
    fprintf(f, '%10.2f %16.2f %18.2f %17.2f \n', T(k), c_p(k),H(k),S(k));
    fprintf(f, '\n');
end

    %Switching back to initial folder
    cd('D:\matlab\bin\win64\CODES')
    A2=fgetl(thermo_data);
    
end

• Molecular Weight

function [m] =molecular_weight(X)
Z = ['H' 'C' 'O' 'N' 'S' 'A'];

A=[1.007 12.0096 15.99903 14.00643 32.059 40.078 ];
m=0;

for i= 1:length(X)
    for j= 1:length(Z)
        if X(i)==Z(j)
            m = m + A(j);
            y=j;
            
        end
    end
    n= str2num(X(i));
           
    if isempty(n)==1
        ;
    else
            m = m+ A(y)*(n-1);
    end


end

• Thermodynamic Properties

function [c_p H S] = temp_coeff(C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,T_low,T_high,T_mid)
R=8.3145;
T=linspace(T_low,T_high,100)
T1= T_mid;
for i=1:length(T)
if  (T(i) <= T1)
    
  c_p(i)  = ((C8) + (C9.*T(i))+(C10.*(T(i).^2))+(C11.*(T(i).^3))+(C12.*(T(i).^4))).*R
   H(i)   = ((C8) +( C9.*(T(i)./2))+(C10.*((T(i).^2)./3))+(C11.*((T(i).^3)./4))+(C12.*((T(i).^4)./5))+(C13./T(i)))*R.*T(i)
   S(i)   = ((C8.*log(T(i)))+(C9.*T(i))+(C10.*((T(i).^2)./2))+(C11.*((T(i).^3)./3))+(C12.*((T(i).^4)./4)+C14)).*R

else   
 
  c_p(i)  = ((C1)+(C2.*T(i))+(C3.*(T(i).^2))+(C4.*(T(i).^3))+(C5.*(T(i).^4))).*R
   H (i)  = ((C1)+( C2.*(T(i)./2))+(C3.*((T(i).^2)/3))+(C4.*((T(i).^3)./4))+(C5.*((T(i).^4)./5))+(C6./T(i)))*R.*T(i)
   S (i)  = ((C1.*log(T(i)))+(C2.*T(i))+(C3.*((T(i).^2)./2))+(C4.*((T(i).^3)./3))+(C5.*((T(i).^4)./4)+C7)).*R
end
end
end

Explanation

1. As each species has its properties specified in the 4 lines wherein, 1st line gives the local temperatures and the rest 3 gives us the temperature coefficient of that species.
2. To separate these values we use commands like fgetl, strfind and strsplit and store them in the suitable variables from each of those lines.
3. Since there are 53 such species we use 'for loop'.
4. Functions are called in to calculate the properties and the molecular weight.
5. Text file is created to save the calculated values.
6. mkdir and cd commands are used to create the dump folders and save the plots and the text file.

Plots and Snapshots

• O2

• N2

• CO2

• Directory

• Text File

• Molecular mass