Parsing NASA thermodynamic data

1.AIM- To parse of thermodynamic data which is collected by NASA. Based on data calculate the enthalpy, entropy and specific heat in terms of temperature by calculating coefficients from data.

2.Introduction- Here we have to write a code which takes our data file “THERMO.dat” as an input.

So firstly we read this file line by line and get information what we want.

For open this file use fopen() command. Now for line by line reading use command fgetl (). In this file first, we have 3 temperature values 300K,1000K and 5000K. For this information use strsplit command to get temperature values.

In the data file every data is saved in the same format like 1^st species name then 14 variables 7 for low-temperature range and 7 for high-temperature range. Values below 3000K are classified as low-temperature range and above this high-temperature range. From string value to num value use command “str2double”.

Specific heat, enthalpy and entropy are calculated according to a given formula.

_{Cp/R = a1 + a2 T + 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 ln(T) + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7}

_{Main code -}

clear all
close all
clc
%universal gas constant
R=8.314;
%opening and reading data file 
f1= fopen('THERMO.dat','r');
%HEADER information
line1=fgetl(f1);
%global temp
temp=fgetl(f1);
A=strsplit(temp,' ');
global_low_temp=str2double(A{2});    %global low temperature
global_mid_temp=str2double(A{3});    %global mid temperature
global_high_temp=str2double(A{4});   %global high temperature
%skipping comments line
des_1=fgetl(f1);
des_2=fgetl(f1);
des_3=fgetl(f1);
%species information
for j=1:53
    
    species=fgetl(f1);
    species_name=species(1:6)
   
   %reading 1st line data 
    line_1=fgetl(f1);
    C=strfind(line_1,'E');
    a1=line_1(1:C(1)+3);
    a1=str2double(a1);   %1st coefficient
    a2=line_1(C(1)+4:C(2)+3);
    a2=str2double(a2);   %2nd coefficient
    a3=line_1(C(2)+4:C(3)+3);
    a3=str2double(a3);    %3rd coefficient
    a4=line_1(C(3)+4:C(4)+3);
    a4=str2double(a4);    %4th coefficient
    a5=line_1(C(4)+4:C(5)+3);
    a5=str2double(a5);    %5th coefficient
    
    %reading 2nd line data
    line_2=fgetl(f1);
    D=strfind(line_2,'E');
    a6=line_2(1:D(1)+3);
    a6=str2double(a6);     %6th coefficient
    a7=line_2(D(1)+4:D(2)+3);
    a7=str2double(a7);     %7th coefficient
    a8=line_2(D(2)+4:D(3)+3);
    a8=str2double(a8);     %8th coefficient
    a9=line_2(D(3)+4:D(4)+3);
    a9=str2double(a9);     %9th coefficient
    a10=line_2(D(4)+4:D(5)+3);
    a10=str2double(a10);     %10th coefficient
    
    %3rd line data reading
    line_3=fgetl(f1);
    D=strfind(line_3,'E');
    a11=line_3(1:D(1)+3);
    a11=str2double(a11);    %11th coefficient
    a12=line_3(D(1)+4:D(2)+3);
    a12=str2double(a12);    %12th coefficient
    a13=line_3(D(2)+4:D(3)+3);
    a13=str2double(a13);    %13th coefficient
    a14=line_3(D(3)+4:D(4)+3);
    a14=str2double(a14);    %14th coefficient
   
   
    %Create temp range for calculation
    T=linspace(300,5000,100);
    %calculation of properties based on low and high temperature range
    for j=1:length(T)
        if T(j)>3000
            cp=(a1+a2.*T+a3.*T.^2+a4.*T.^3+a5.*T.^4).*R;
            H=(a1+(a2.*T)/2+(a3.*T.^2)/3+(a4.*T.^3)/4+(a5.*T.^4)/5+a6./T).*R.*T;
            S=(a1.*log(T)+a2.*T+(a3.*T.^2)/2+(a4.*T.^3)/3+(a5.*T.^4)/4+a7).*R;
        else
            cp=(a8+a9.*T+a10.*T.^2+a11.*T.^3+a12.*T.^4).*R;
            H=(a8+(a9.*T)/2+(a10.*T.^2)/3+(a11.*T.^3)/4+(a12.*T.^4)/5+a13./T).*R.*T;
            S=(a8.*log(T)+a9.*T+(a10.*T.^2)/2+(a11.*T.^3)/3+(a12.*T.^4)/4+a14).*R;
        end
    end
    %saving as a folder
    spec=char(species_name);
    fd=sprintf('%s',spec);
    cur_dir=pwd;
    mkdir(spec);
    cd(spec)
    %plotting 
    %Temperature vs enthalpy
    figure(1)
    plot(T,H,'linewidth',4,'color','b')
    title('species_name')
    xlabel('Temperature')
    ylabel('Enthalpy')
    
    %Temperature vs entropy
    figure(2)
    plot(T,S,'linewidth',4,'color','r')
    title('species_name')
    xlabel('Temperature')
    ylabel('Entropy')
    
    %Temperature vs specific heat
    figure(3)
    plot(T,cp,'linewidth',4,'color','y')
    title('species_name')
    xlabel('Temperature')
    ylabel('specific heat')
   
    %save the plot
    saveas(1,'enthalpy.jpg');
    saveas(2,'entropy.jpg');
    saveas(3,'specificheat.jpg');
   cd(cur_dir);
    %molecular mass
    MM=molecular_mass(species_name)
   
end
fclose(f1);
    

_{For calculation of Molecular mass of species, we write a separate code.}

_{supplement code-}

%function for molecular mass
function MM= molecular_mass(species_name)
species=upper(species_name);
atoms={'H','C','N','O','AR'};
atomicweight=[1 12 14 16 40];
MM=0;
for i=1:length(species)
    for j=1:length(atoms)
        if strcmp(species(i),atoms(j))
            MM=MM+atomicweight(j);
            m=j;
        end
    end
    n= str2double(species(i));
    if n>1
       MM= MM+atomicweight(m)*(n-1);
    end
end
end
        

_{Plot for O2-}

_{Plot for N2}

_{Plot for CO2-}

_{Command Window-MM value is printed for every species.}