Project 1 - Parsing NASA thermodynamic data

CODE:

clear all
close all
clc

fname = input('Please enter the filename with extension:','s');
fidr = 1;

while fidr ~= -1
    
    [fidr,ermsg] = fopen(fname,'r'); %Open the file
    fprintf('File ID:%d\n',fidr);
    if (fidr == -1)
        
        disp(ermsg)
        fprintf('Do you want to again enter filename ? (1 - YES, 0 - NO)')
        option = input(':');
        
        if option
            
            fprintf('Please enter a valid file name')
            fname = input(':','s');
            fid = 1;
            
        else
            
            disp('Please try again !')
            break
            
        end
        
    else
        
        break
        
    end
    
end

if fidr ~= -1
    species_cal(fidr);
    fclose(fidr);
end

disp('Shutting down!')
pause(1)

%CODE TO READ THE THERMO30  FILE AND CALCULATE THE THERMODYNAMIC PROPERTIES
%CP, H AND S FOR ALL THE SPECIES PRESENT IN THE FILE ALONG WITH THEIR
%MOLECULAR WEIGHTS-BY ANIL DUBEY

function [] = species_cal(fidr)

%Store and display the global temperature range
disp(fgetl(fidr));
fprintf('\n')
text = fgetl(fidr);% To read a line of text
arr = strsplit(text);%To split a text using space '  ' as delimiter
global_temp = [str2double(arr{2}), str2double(arr{3}), str2double(arr{4})];
fprintf('GLOBAL TEMPERATURE RANGE: ')
fprintf('%.3f  %.3f  %.3f\n',global_temp)
fprintf('\n')
disp('DATA OBTAINED FROM:')
fprintf('\n')
text = fgetl(fidr);
disp(text);
fprintf('\n')
text = fgetl(fidr);
disp(text);
text = fgetl(fidr);

%Taking input from user
fprintf('\n')
disp('OPTIONS:')
fprintf('\n')
disp('1. Type species name and press enter for properties of individual species.')
fprintf('\n')
disp('2. Just press enter for obtaining data for all the species in the THERMO30 file.')
fprintf('\n')
species = input('Waiting for input: ','s');

%Convert the input to uppercase to avoid character mismatch (if any)
species = upper(species);

%Create arrays to store data for a given species
local_temp = zeros(1,3);
highT_coeff = zeros(1,7);
lowT_coeff = zeros(1,7);

if isempty(species) %Executed when user just presses the Enter key
    
%This if block calculates all the properties for all the species and saved
%them in seperate folders
   
    text = fgetl(fidr); %Read the first line of the first species in the file THERMO30.dat

    while ~(text(1)=='E') %The loop runs until the last line of the file which is 'END'
        
        arr = strsplit(text);
        species_weight = 0;
        species = arr{1}; %The strsplit function returns a cell
        element = {};%Cell to store the names of various elements in a species
        number = [];%Array to store the number of corresponding element in a species 
        tag = 1;
        element_no = 1;
        number_no = 1;
        name_index = 1;
        
%We are looping 25 through 44 because the no of elements in each species
%are listed beyond 25th character in the first line of data for each
%species. This can be checked by counting the number of characters
%including the white spaces. We need the exact number of elements to
%calculate the molecular weight.(This was the most challenging task)

        for i = 25:44
            
            
            if isletter(text(i))
                    
                element{element_no}(name_index) = text(i);
                name_index = name_index + 1;
                tag = 1;
                    
            elseif (text(i) == ' ') && tag
                
                tag = 0;
                name_index = 1;
                element_no = element_no + 1;

                
            elseif (~(text(i) == ' ') && ~(isletter(text(i))))
                
                number(number_no) = str2double(text(i));
                number_no = number_no + 1;
                
            end
            
        end

        %Calculate the molecular weight using the info stored in element
        %and number data structures (see lines 56 and 57)
        for i1 = 1:length(number)
            species_weight = add_weight(element{i1},number(i1),species_weight);
        end
        
    disp(' ')
    fprintf('Molecular weight of %s = %.2f\n',species,species_weight)
    disp(' ')
    %Access the local temperature range specific to each species
    local_temp = [str2double(arr{length(arr)-3}), str2double(arr{length(arr)-2})...
            str2double(arr{length(arr)-1})];
    %Access all the 14 co-efficients for the species and store it in the
    %data matrix
    data = fscanf(fidr,'%f %f %f %f %f %f',[6,3]);
    
    %Taking the transpose of data matrix as the row-col order changes
    %while reading through fscanf
    data = data.';
    %Accessing and storing all the required data to compute the properties
    highT_coeff = [data(1,1:5),data(2,1:2)];
    lowT_coeff = [data(2,3:5),data(3,1:4)];
    data_calculator(species,local_temp,highT_coeff,lowT_coeff)
    text = fgetl(fidr);
    
    end
  
else
    
%This else block is executed when the user enters a name of a specific
%species for which properties are required. The coding logic to calculate 
%the molecular weight and thermodynamic properties is same as that of the 
%if block.
    
    flag = 1;
    
    while flag
        
    text = fgetl(fidr);
    species_weight = 0;
    arr = strsplit(text);
    element = {};
    number = [];
    tag = 1;
    element_no = 1;
    number_no = 1;
    name_index = 1;

        if (arr{length(arr)} == '1' && strcmp(arr{1},species))

            for i2 = 25:44 
            
            
                if isletter(text(i2))

                    element{element_no}(name_index) = text(i2);
                    name_index = name_index + 1;
                    tag = 1;

                elseif (text(i2) == ' ') && tag

                    tag = 0;
                    name_index = 1;
                    element_no = element_no + 1;


                elseif (~(text(i2) == ' ') && ~(isletter(text(i2))))

                    number(number_no) = str2double(text(i2));
                    number_no = number_no + 1;

                end
            
            end
            
            for i3 = 1:length(number)
                species_weight = add_weight(element{i3},number(i3),species_weight);
            end
        
            disp(' ')
            fprintf('Molecular weight of %s = %.2f\n',species,species_weight)
            flag = 0;
            local_temp = [str2double(arr{length(arr)-3}), str2double(arr{length(arr)-2})...
                str2double(arr{length(arr)-1})];
            data = fscanf(fidr,'%f %f %f %f %f %f',[6,3]);
            data = data.';
            highT_coeff = [data(1,1:5),data(2,1:2)];
            lowT_coeff = [data(2,3:5),data(3,1:4)];

        else

            for i4 = 1:3
                fgetl(fidr);
            end

        end
    
    end
    
    data_calculator(species,local_temp,highT_coeff,lowT_coeff)
    
end

disp('All figures saved successfully!')
pause(1)

%Function to add the weight of each element in a given species to calculate
%the total weight

function [species_weight] = add_weight(ele,count,species_weight)

    switch ele
        case 'C'
            species_weight = species_weight + 12*count;
        case 'H'
            species_weight = species_weight + 1*count;
        case 'N'
            species_weight = species_weight + 14*count;
        case 'O'
            species_weight = species_weight + 16*count;
        case 'AR'
            species_weight = species_weight + 40*count;
        otherwise
            error('Element not available!')
    end
    
end

%Function to calculate all the properties Cp, H and S for a given species

function [] = data_calculator(species,local_temp,highT_coeff,lowT_coeff)

    R = 8.314; %Universal gas constant in J/mol K 
    T_low = local_temp(1):1:local_temp(3);
    T_high = local_temp(3)+1:1:local_temp(2);
    T = [T_low,T_high];
    Cp = zeros(1,length(T));
    H = zeros(1,length(T));
    S = zeros(1,length(T));

    for i5  = 1:length(T)

        if T(i5) > local_temp(3)

            a1 = highT_coeff(1);
            a2 = highT_coeff(2);
            a3 = highT_coeff(3);
            a4 = highT_coeff(4);
            a5 = highT_coeff(5);
            a6 = highT_coeff(6);
            a7 = highT_coeff(7);

        else

            a1 = lowT_coeff(1);
            a2 = lowT_coeff(2);
            a3 = lowT_coeff(3);
            a4 = lowT_coeff(4);
            a5 = lowT_coeff(5);
            a6 = lowT_coeff(6);
            a7 = lowT_coeff(7);
        end

        Cp(i5) = (a1 + a2*T(i5) + a3*T(i5)^2 + a4*T(i5)^3 + a5*T(i5)^4)*R;
        H(i5) = (a1 + (a2*T(i5))/2 + (a3*T(i5)^2)/3 + (a4*T(i5)^3)/4 + (a5*T(i5)^4)/5 + a6/T(i5))*R*T(i5);
        S(i5) = (a1*log(T(i5)) + a2*T(i5) + (a3*T(i5)^2)/2 + (a4*T(i5)^3)/3 + (a5*T(i5)^4)/4 + a7)*R;

    end
    
    folder_name = strcat('fig_',species);%Create name of the folder for a given species name
    mkdir(folder_name); %Make a new folder to store the figures of a given species
    path = pwd; %Set path to the present working directory
    
    %Cp vs T plot
    figure(1)
    plot(T,Cp,'-k','LineWidth',2)
    xlabel('$Temperature (K)$','Interpreter','latex','FontSize',14)
    ylabel('$Cp (J/mol-K)$','Interpreter','latex','FontSize',14)
    title(strcat('Cp for ',' ',species),'Interpreter','latex','FontSize',16)
    %Create Plot name
    name = strcat(species,'_Cp','.jpg');
    %Create full name of the path for saving the figure
    full_name = fullfile(path,folder_name,name);
    print(gcf, full_name, '-djpeg', '-r600');
    
    %H vs T plot
    figure(2)
    plot(T,H,'-k','LineWidth',2)
    xlabel('$Temperature (K)$','Interpreter','latex','FontSize',14)
    ylabel('$H (J/mol-K)$','Interpreter','latex','FontSize',14)
    title(strcat('Enthalpy change for: ',' ',species),'Interpreter','latex','FontSize',16)
    %Create Plot name
    name = strcat(species,'_H','.jpg');
    %Create full name of the path for saving the figure
    full_name = fullfile(path,folder_name,name);
    print(gcf, full_name, '-djpeg', '-r600');

    %S vs T plot
    figure(3)
    plot(T,S,'-k','LineWidth',2)
    xlabel('$Temperature (K)$','Interpreter','latex','FontSize',14)
    ylabel('$S (J/mol-K)$','Interpreter','latex','FontSize',14)
    title(strcat('Entropy change for: ',' ',species),'Interpreter','latex','FontSize',16)
    %Create Plot name
    name = strcat(species,'_S','.jpg');
    %Create full name of the path for saving the figure
    full_name = fullfile(path,folder_name,name);
    print(gcf, full_name, '-djpeg', '-r600');
    
end

end

OUTPUT SCREEN SCREENSHOT:

FOLDER SCREENSHOT:

LOGIC:

1. The file has to be opened through fopen() function in read mode. If the file does not exist file ID returned by fopen is -1 which is used in a while loop to decide if the file name is incorrect.

2. The user is again given a chance to give the correct file name in case the file could not be opened. If the error persists then the file does not exist in the present working directory.

3. If the file exists then the file is opened in the read mode and the file ID is passed to species_cal() function to carry out the computation.

4. First five lines of the THERMO30.dat file were displayed on the screen for the user.

5. The user is then given option to enter a species name to generate its plot or to press enter only to generate the plot for all the species.

6. Rest of the code is self explanatory. One may look at the comments within the code for details.

7. The most challenging part was to determine the number of elements in each species to calculate the molecular weight.

The first piece of code developed for this purpose is as follows:

The molar weights of each element was stored in a structure named mol_weight as mol_weight.C = 12; mol_weight.H = 1; mol_weight.N = 14 and so on. The variable species_weight was initially assigned a value of 0.

for i = 1:length(species)

        if species(i) == 'C'

            try
                if isnumeric(species(i+1))
                    species_weight = species_weight + mol_weight.C*str2double(species(i+1));
                else
                    species_weight = species_weight + mol_weight.C;
                end
            catch ME
                if strcmp(ME.identifier,'MATLAB:badsubscript')
                    species_weight = species_weight + mol_weight.C;
                else
                    rethrow(ME)
                end
            end

        elseif species(i) == 'H'

            try
                if isnumeric(species(i+1))
                    species_weight = species_weight + mol_weight.H*str2double(species(i+1));
                else
                    species_weight = species_weight + mol_weight.H;
                end
            catch ME
                if strcmp(ME.identifier,'MATLAB:badsubscript')
                    species_weight = species_weight + mol_weight.H;
                else
                    rethrow(ME)
                end
            end

        elseif species(i) == 'N'

            try
                if isnumeric(species(i+1))
                    species_weight = species_weight + mol_weight.N*str2double(species(i+1));
                else
                    species_weight = species_weight + mol_weight.N;
                end
            catch ME
                if strcmp(ME.identifier,'MATLAB:badsubscript')
                    species_weight = species_weight + mol_weight.N;
                else
                    rethrow(ME)
                end
            end

        elseif species(i) == 'O'

            try
                if isnumeric(species(i+1))
                    species_weight = species_weight + mol_weight.O*str2double(species(i+1));
                else
                    species_weight = species_weight + mol_weight.O;
                end
            catch ME
                if strcmp(ME.identifier,'MATLAB:badsubscript')
                    species_weight = species_weight + mol_weight.O;
                else
                    rethrow(ME)
                end
            end

        elseif species(i) == 'AR'

            try
                if isnumeric(species(i+1))
                    species_weight = species_weight + mol_weight.AR*str2double(species(i+1));
                else
                    species_weight = species_weight + mol_weight.AR;
                end
            catch ME
                if strcmp(ME.identifier,'MATLAB:badsubscript')
                    species_weight = species_weight + mol_weight.AR;
                else
                    rethrow(ME)
                end
            end
        end
    end

The drawbacks of the above code piece are as follows:

i. It was unable to successfully calculate the weight of species whose name was in double characters such as 'AR'.

ii. Adding a new element required adding 13 new lines in the code.

iii. As it was using the species name itself (i.e. HCNN, H2O etc.) to identify the number of elements erroneous results were produced with species name such as 'AR' were encountered.

The second code piece developed for the same purpose is as follows:

This used the first line of every species to find out the number of molecules and used the add_weight function given in the final code to add the weight of each element to the species_weight variable. This significantly reduced the size of the code and made it easy to add a new element to the code for the calculation of weight. But this code also failed to successfully calculate the weights of elements with names in double characters such as AR. Indices 24 to 44 are used in the for loop because the no of elements in each species are found within these indices only in the first line corresponding to details of each species in the thermo30.dat file.

for i = 25:44
            if (isletter(text(i)) && isletter(text(i+1)))
                num = 0;
                j = i+2;
                while ~num
                    if text(j) ~= ' '
                        num = str2double(text(j));
                    end
                    j = j+1;
                end

                species_weight = add_weight(text(i:i+1),num,species_weight);

            elseif (text(i+1) == ' ') && isletter(text(i))
                num = 0;
                j = i+1;
                while ~num

                    if text(j) ~= ' '
                        num = str2double(text(j));
                    end
                    
                    j = j+1;

                end

                species_weight = add_weight(text(i),num,species_weight);
                
            end
    end

Finally, the final code which ran successfully for all the species including AR, was smaller in size, and easier to modify to include any new element for calculation of weight, was arrived at in the final code attached in the beginning. In this code instead of directly using the first line of every species' data to calculate the weight, first the names of each species and their corresponding number of atoms in each molecule were saved in a cell names elements and array named number. The pair was finally fed to the add_weight() function to calculate the molecular weight successfully.