Project 1 - Parsing NASA thermodynamic data

Main Program:

clear
close
clc
% universal gas constant
R = 8.314;               
%using fopen to open the file
f = fopen('THERMO.dat','r');

for i=1:5
    fgetl(f);
end

for j = 1:53
    %extraction of local temperatures
    a = fgetl(f);
    b = strsplit(a,' ');
    c=length(b);
    
    % to extract species name
    name=b{1};
    
    local_low = str2num(b{c-3}); 
    local_high = str2num(b{c-2});
    local_mid = str2num(b{c-1});
    
    %calculation of cp, h, s
    T = linspace(local_low,local_high,100);
    
    %extraction of co-efficients
    coeff_line1 = fgetl(f);
    pos_e = strfind(coeff_line1, 'E');
    a1 = str2num(coeff_line1(1:pos_e(1)+3));
    a2 = str2num(coeff_line1((pos_e(1)+4):(pos_e(2)+3)));
    a3 = str2num(coeff_line1((pos_e(2)+4):(pos_e(3)+3)));
    a4 = str2num(coeff_line1((pos_e(3)+4):(pos_e(4)+3)));
    a5 = str2num(coeff_line1((pos_e(4)+4):(pos_e(5)+3)));
    
    coeff_line2 = fgetl(f);
    pos_e1 = strfind(coeff_line2, 'E');
    a6 = str2num(coeff_line2(1:pos_e1(1)+3));
    a7 = str2num(coeff_line2((pos_e1(1)+4):(pos_e1(2)+3)));
    a8 = str2num(coeff_line2((pos_e1(2)+4):(pos_e1(3)+3)));
    a9 = str2num(coeff_line2((pos_e1(3)+4):(pos_e1(4)+3)));
    a10 = str2num(coeff_line2((pos_e1(4)+4):(pos_e1(5)+3)));
    
    coeff_line3 = fgetl(f);
    pos_e2 = strfind(coeff_line3, 'E');
    a11 = str2num(coeff_line3(1:pos_e2(1)+3));
    a12 = str2num(coeff_line3((pos_e2(1)+4):(pos_e2(2)+3)));
    a13 = str2num(coeff_line3((pos_e2(2)+4):(pos_e2(3)+3)));
    a14 = str2num(coeff_line3((pos_e2(3)+4):(pos_e2(4)+3)));
    
    %calculation of cp,h,s
    cp = NASA_specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,local_mid);
    h = NASA_enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,local_mid);
    s = NASA_entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,local_mid);
    
    %calculating molecular weight of species
    MW = molecular_weight(name);
    k=sprintf('%s mol_wt=%f',name,MW); 
    disp(k)
    
    % plotting
    cd('E:\MatLab\Nasa project')
    
    mkdir(name)
    cd(name)
    
    figure (1)
    plot(T,cp,'linewidth',2,'color','g');
    xlabel('Local Temperature [K]');
    ylabel('Specific heat [J/kg-k]');
    grid on
    title('Local Temperature vs Specific heat');
    
    figure (2)
    plot(T,h,'linewidth',2,'color','b');
    xlabel('Local Temperature [K]');
    ylabel('Enthalpy [J/kg]');
    grid on
    title('Local Temperature vs Enthalpy');
    
    figure (3)
    plot(T,s,'linewidth',2,'color','r');
    xlabel('Local Temperature [K]');
    ylabel('Entropy [j/k]');
    grid on
    title('Local Temperature vs Entropy');
    
    saveas(figure(1),'specific heat.jpg')
    saveas(figure(2),'Enthalpy.jpg')
    saveas(figure(3),'Entropy.jpg')

end

Specific Heat (Cp):

%calculation of specific heat
function cp = NASA_specificheat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,local_mid)
%if the temperature is more it will take max coefficients
    if T <= local_mid
       cp = ((a8+(a9*T)+(a10*T.^2)+(a11*T.^3)+(a12*T.^4))*R);
%if the temperature is more it will take min coefficients
    else    
         cp = ((a1+(a2*T)+(a3*T.^2)+(a4*T.^3)+(a5*T.^4))*R);
    end
end

Entropy(s):

%calculation of entropy
function s = NASA_entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,local_mid)
    if T <= local_mid
        s =( R*(a8*log(T)+(a9*T)+a10*(T.^2/2)+a11*(T.^3/3)+a12*(T.^4/4)+a14));
    else    
        s = (R*(a1*log(T)+(a2*T)+a3*(T.^2/2)+a4*(T.^3/3)+a5*(T.^4/4)+a7));
    end
end

Enthalpy(h):

%calculation of enthalpy
function h= NASA_enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,T,R,local_mid)
    if T <= local_mid  
        h =( R*(a8*T+a9*(T.^2/2)+a10*(T.^3/3)+a11*(T.^4/4)+a12*(T.^5/5)+a13));
    else    
        h = (R*(a1*T+a2*(T.^2/2)+a3*(T.^3/3)+a4*(T.^4/4)+a5*(T.^5/5)+a6));
    end
end

Molecular Weight:

%Function to calculate molecular weight
function MW =molecular_weight(name)
    atoms = ['H','C','O','N','R'];
    atomic_weight= [1 12 16 14 40];
    MW = 0;
    for i = 1:length(name)
        for j = 1:length(atoms)
            if strcmp(name(i),atoms(j))
                MW = MW +atomic_weight(j);
                position = j;
            end
        end
            n = str2num(name(i));
        if n>1
            MW = MW + atomic_weight(position)*(n-1);
        end
    end
end

Command window:

O mol_wt=16.000000
O2 mol_wt=32.000000
H mol_wt=1.000000
H2 mol_wt=2.000000
OH mol_wt=17.000000
H2O mol_wt=18.000000
HO2 mol_wt=33.000000
H2O2 mol_wt=34.000000
C mol_wt=12.000000
CH mol_wt=13.000000
CH2 mol_wt=14.000000
CH2(S) mol_wt=14.000000
CH3 mol_wt=15.000000
CH4 mol_wt=16.000000
CO mol_wt=28.000000
CO2 mol_wt=44.000000
HCO mol_wt=29.000000
CH2O mol_wt=30.000000
CH2OH mol_wt=31.000000
CH3O mol_wt=31.000000
CH3OH mol_wt=32.000000
C2H mol_wt=25.000000
C2H2 mol_wt=26.000000
C2H3 mol_wt=27.000000
C2H4 mol_wt=28.000000
C2H5 mol_wt=29.000000
C2H6 mol_wt=30.000000
CH2CO mol_wt=42.000000
HCCO mol_wt=41.000000
HCCOH mol_wt=42.000000
H2CN mol_wt=28.000000
HCN mol_wt=27.000000
HNO mol_wt=31.000000
N mol_wt=14.000000
NNH mol_wt=29.000000
N2O mol_wt=44.000000
NH mol_wt=15.000000
NH2 mol_wt=16.000000
NH3 mol_wt=17.000000
NO mol_wt=30.000000
NO2 mol_wt=46.000000
HCNO mol_wt=43.000000
HOCN mol_wt=43.000000
HNCO mol_wt=43.000000
NCO mol_wt=42.000000
CN mol_wt=26.000000
HCNN mol_wt=41.000000
N2 mol_wt=28.000000
AR mol_wt=40.000000
C3H8 mol_wt=44.000000
C3H7 mol_wt=43.000000
CH3CHO mol_wt=44.000000
CH2CHO mol_wt=43.000000
>> 

Code logic:

fgetl, fopen are used to open and read the file, as we can say fgetl will consider one line for each time it runs. strsplit is used to split the line based on the space between each element(word). str2num is used to consider the numeric values in a line based on their position and spaces, if it is not a number it will show as an error. coeff_line(1,2,3) is used to define each line which has numeric data(14 coefficients) and strfind is used to find the position of E and with the help of it, by using str2num coefficients were divided as these doesn't have spaces between them.

for the calculation of Cp, s, h formulas for each were written as a function and it was added into the main program. Molecular weight calculation was also written as a function by defining the atomic weight (approx) and giving the formula.

to draw the plots for this data, I have created separate folder for each species(using mkdir,cd) and it was defined in this program to save all of them. all the plots were labeled with x,y axis, title and thtey were saves as jpg images. to create a new folder for each species,, I have used name = b{1}, as the name od species is the first word of the line.

for loop is used to run the program for 52 more iterations to calculate for all 53 species.

I have attached plots from species (CO2, O2, N2) as asked and to show that I have given title, x and y axis and I have used the local temperature not the global temperature. I am attaching the command window as well to show the calculated molecular weights.