BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Project 1 - Parsing NASA thermodynamic data

RRaghupathruni Venkatesh6 Dec 2021
MATLAB
% Matlab program to parse the NASA data to find Specific heat(Cp), ENtropy(s), %Enthalpy(h) and Molecular weight(MW) for each species

clear all
close all
clc

f1 = fopen('THERMO.dat','r');
first_line = fgetl(f1); %reading the first line from data

% Obtaining global maximum, medium and minimum temperatures 
global_temp = fgetl(f1)  %reading 2nd line from data
temp_range = strsplit(global_temp,' ') %splitting second line using blank space as a  
                                       %delimiter

% Obtaining temperature values from splitting data 
global_min_temp = str2num(temp_range{2})
global_med_temp = str2num(temp_range{3})
global_max_temp = str2num(temp_range{4})

% Creating a vector of 400 evenly spaced temperature points

R = 8.314 %univarsal gas constant
comments1 = fgetl(f1); %reading third line from data
comments2 = fgetl(f1); %reading fourth line frtom data
comments3 = fgetl(f1); %reading fifth line from data

for i = 1:53
    line1 = fgetl(f1);
    line_1 = strsplit(line1,' ')
    T_med = str2double(line_1(length(line_1)-1));
    T_max = str2double(line_1(length(line_1)-2));
    T_min = str2double(line_1(length(line_1)-3));
    T = linspace(T_min,T_max,400);
    element_name = strtok(line1);
    
    line2 = fgetl(f1);
    A = strfind(line2,'E');
    
    line3 = fgetl(f1);
    B = strfind(line3,'E');
    
    line4 = fgetl(f1);
    C = strfind(line4,'E');
    
    a1 = str2num(line2(1:A(1)+3));
    a2 = str2num(line2(A(1)+4:A(2)+3));
    a3 = str2num(line2(A(2)+4:A(3)+3));
    a4 = str2num(line2(A(3)+4:A(4)+3));
    a5 = str2num(line2(A(4)+4:A(5)+3));
    a6 = str2num(line3(1:B(1)+3));
    a7 = str2num(line3(B(1)+4:B(2)+3));
    
    b1 = str2num(line3(B(2)+4:B(3)+3));
    b2 = str2num(line3(B(3)+4:B(4)+3));
    b3 = str2num(line3(B(4)+4:B(5)+3));
    b4 = str2num(line4(1:C(1)+3));
    b5 = str2num(line4(C(1)+4:C(2)+3));
    b6 = str2num(line4(C(2)+4:C(3)+3));
    b7 = str2num(line4(C(3)+4:C(4)+3));
   for u = 1:length(T) 
    % Specific heat
    cp(u) = specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T(u),R,global_med_temp);
    
    % Enthalpy
    h(u) = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T(u),R,global_med_temp);

    % Entropy
    s(u) = entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T(u),R,global_med_temp);
   end

    % Calculating molecular weight of species
    Molecular_wt(i) = molecular_weight(element_name);
    
    % Saving the data
    
    cd('C:\Users\venka\Onedrive\Documents\MATLAB\nasa_plots')
    dir_name = element_name;
    mkdir(dir_name);

    % Plotting the data

    figure(1)
    plot(T,cp,'linewidth',2,'color','b')
    title(sprintf('Specific heat Vs Temperature for %s',element_name))
    xlabel('TEMPERATURE')
    ylabel('SPECIFIC HEAT')
    grid on
    
    figure(2)
    plot(T,h,'linewidth',2,'color','r')
    title(sprintf('Enthalpy Vs Temperature for %s',element_name))
    xlabel('TEMPERATURE')
    ylabel('ENTHALPY')
    grid on
    
    figure(3)
    plot(T,s,'linewidth',2,'color','g')
    title(sprintf('Entropy Vs Temperature for %s',element_name))
    xlabel('TEMPERATURE')
    ylabel('ENTROPY')
    grid on
    
    %pause the loop so we can see each species plot in real time
    pause(0.0001)
    cd(dir_name);
    saveas(1,'Specific heat.png');
    saveas(2,'Enthalpy.png');
    saveas(3,'Entropy.png');

end

fclose(f1)

* The below functions are used in the main code to find out the values of Specific heat(Cp), Enthalpy(h) and Entropy(s) for each specie

% Function to calculate Specific heat(cp)

function cp = specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,global_med_temp)

% Comparing with global temp
if T>global_med_temp
    cp = R.*(a1+a2.*T+a3.*T^2+a4.*T.^3+a5.*T.^4);
else
    cp = R.*(b1+b2.*T+b3.*T.^2+b4.*T.^3+b5.*T.^4);
end

---------------------------------------------------------------------------------------

% Function to claculate Enthalpy(h)

function h = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,globa_med_temp)

if T>globa_med_temp
    h = R.*(a1.*T+(a2.*T.^2)./2+(a3.*T.^3)./3+(a4.*T.^4)./4+(a5.*T.^5)./5+a6);
else 
    h = R.*(b1.*T+(b2.*T.^2)./2+(b3.*T.^3)./3+(b4.*T.^4)./4+(b5.*T.^5)./5+b6);
end

---------------------------------------------------------------------------------------

% Function to calculate Entropy(s)

function s = entropy(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,global_med_temp)
if T<=global_med_temp
    s =  R*(a1.*log(T)+a2.*T+a3.*(T.^2)./2+a4.*(T.^3)./3+a5.*(T.^4)./4+a7);
else
    s = R*(b1.*log(T)+b2.*T+b3.*(T.^2)./2+b4.*(T.^3)./3+b5.*(T.^4)./4+b7);
end

* Function for Molecular Weight(MW)

- For molecular weight calculation, created a function called 'molecular_weight'
- This function takes the weight of hydrogen, carbon, nitrogen and argon as inputs. All the elements consists of all these atoms.
- I have used 'strcmp' command to compare elements with species. If they both matches then the molecular weight will be calculated.
- Take the string of species name and position to calculate their weight. Here 'n' is number of molecules of elements present in the species. - The function commands are given below.

%Function to calculate molecular weight

function MW =molecular_weight(element_name)
atoms = ['H','C','N','O','A'];
     fprintf('Molecular Weight of ');
    fprintf(element_name);
    fprintf(' : ');
    atomic_weight = [1 12 14 15 39];
    MW = 0;
    for i = 1:length(element_name)
        for j = 1:length(atoms)
            if strcmp(element_name(i),atoms(j))
                MW = MW +atomic_weight(j);
                position = j;
            end
        end
        n = str2num(element_name(i));
        if n>1
            MW = MW + atomic_weight(position)*(n-1);
        end
    end
    fprintf('%d',MW);
    disp(' ');
    end

>> Molecular Weights:

Molecular Weight of O : 15
Molecular Weight of O2 : 30
Molecular Weight of H : 1
Molecular Weight of H2 : 2
Molecular Weight of OH : 16
Molecular Weight of H2O : 17
Molecular Weight of HO2 : 31
Molecular Weight of H2O2 : 32
Molecular Weight of C : 12
Molecular Weight of CH : 13
Molecular Weight of CH2 : 14
Molecular Weight of CH2(S) : 14
Molecular Weight of CH3 : 15
Molecular Weight of CH4 : 16
Molecular Weight of CO : 27
Molecular Weight of CO2 : 42
Molecular Weight of HCO : 28
Molecular Weight of CH2O : 29
Molecular Weight of CH2OH : 30
Molecular Weight of CH3O : 30
Molecular Weight of CH3OH : 31
Molecular Weight of C2H : 25
Molecular Weight of C2H2 : 26
Molecular Weight of C2H3 : 27
Molecular Weight of C2H4 : 28
Molecular Weight of C2H5 : 29
Molecular Weight of C2H6 : 30
Molecular Weight of CH2CO : 41
Molecular Weight of HCCO : 40
Molecular Weight of HCCOH : 41
Molecular Weight of H2CN : 28
Molecular Weight of HCN : 27
Molecular Weight of HNO : 30
Molecular Weight of N : 14
Molecular Weight of NNH : 29
Molecular Weight of N2O : 43
Molecular Weight of NH : 15
Molecular Weight of NH2 : 16
Molecular Weight of NH3 : 17
Molecular Weight of NO : 29
Molecular Weight of NO2 : 44
Molecular Weight of HCNO : 42
Molecular Weight of HOCN : 42
Molecular Weight of HNCO : 42
Molecular Weight of NCO : 41
Molecular Weight of CN : 26
Molecular Weight of HCNN : 41
Molecular Weight of N2 : 28
Molecular Weight of AR : 39
Molecular Weight of C3H8 : 44
Molecular Weight of C3H7 : 43
Molecular Weight of CH3CHO : 43
Molecular Weight of CH2CHO : 42

>> References:

1. https://www.mathworks.com/help/matlab/matlab_prog/parse-function-inputs.html

2. https://www.mathworks.com/help/matlab/ref/fscanf.html

3. https://www.mathworks.com/help/instrument/scanstr.html

4. https://www.mathworks.com/help/matlab/ref/mkdir.html

>> Step wise explanation of code:

The main objective of the project is to write a function to extract 14 coefficients (The first seven coefficients are high-temperature coefficients and the second seven coefficients are low-temperature coefficients) with the input information of given species and calculate the Molecular weight (MW), Specific heat (Cp), Entropy (s) and Enthalpy (h) for each species with help of NASA thermodynamic polynomial data ‘THERMO’.

  1. I have used ‘fopen’ command to open NASA’s experimental data (THERMO).
  2. The first line consists text (THERMO), ‘fgetl’ command is used to read the first line in the data.
  3. The second line consist of three temperature values, here also ‘fgetl’ command is used to read the second line and stored the values in the variable called global_temp.
  4. The ‘strsplit’ command is used to split the data using given delimiter. The ‘blank space’ is a delimiter, the temperature range consist of three values in a vector.
  5. Obtained the global minimum, medium and maximum temperatures from the splitted data, here ‘str2num’ command is used to convert character array of temperature into numeric array of the temperature.
  6. Now created a vector of 400 evenly spaced temperature points by using ‘linspace’ command and defined universal gas constant value used in the calculations.
  7. Again ‘fgetl’ command is used to read the third, fourth and fifth lines which are text comments in the data file THERMO, to be read in text form.
  8. From line 6 to 218, there is a set of 4 lines, first line consists of the name of the species and second, third and fourth lines will have experimental data.
  9. Line 1, I have used ‘fgetl’ command to read line 6 and stored the line 6 data in the “line­_1” variable.
  10. ‘strtok’ command is used to get the selected part of the string, used to select the specie which have ‘O’ and store the respective string (Specie) in line_1 variable.
  11. Line 2, I have used ‘fgetl’ command to read line 7 and stored the data in variable line_2
  12. To find the places of ‘E’ in line 7, ‘strfind’ command is used and stored respective species in an array named ‘A’.
  13. Similarly, I have repeated the same procedure for next two lines (line three and line four) and saved them in the arrays called ‘B’ and ‘C’.
  14. With the help of ‘str2num’ command, extracted all the coefficients (a1, a2, a3, a4, a5, a6, a7, b1, b2, b3, b4, b5, b6 and b7).
  15. This entire process needs to repeat from line 6 to line 218 for all 53 species, a ‘for loop’ is used for above process.
  16. By calling the functions [specific heat (cp), enthalpy (h), entropy (s) and molecular weight (MW)] into the main code to find out the values of cp, h, s and MW for each species.
  17. ‘plot’ command is used to plot the

Temperature and specific heat

Temperature and entropy

Temperature and enthalpy for all the species

  1. Used ‘cd’ command to create a new folder for every iteration of loop and it creates a folder with the species name and pasted its directory in cd
  2. ‘mkadir’ command is used to provide name for the folder which is created in the earlier step
  • The below commands are used to save all the plots in a separate folder

cd('C:\Users\venka\Onedrive\Documents\MATLAB\nasa_plots')

dir_name = element_name;

mkdir(dir_name)

  1. ‘saveas’ command is used to save the plots in ‘png’ format, the commands are given below

cd(dir_name);

saveas(1,'Specific heat.png');

saveas(2,'Enthalpy.png');

saveas(3,'Entropy.png');

  1. At the end, closed the NASA’s thermodynamic experimental data file (THERMO) with the help of ‘fclose’ command.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.