Project 1 - Parsing NASA thermodynamic data

%Parsing NASA thermodynamic data

clear all
close all
clc

%Universal Gas Constant in KJ/K-mol-K
R = 8.314;

%Importing the external file
K = fopen('THERMO.dat','r');
line1 = fgetl(K);
line2 = fgetl(K);

%Parcing the global temp values
global_values = strsplit(line2,' ');
global_low_temp = str2num(global_values{2})
global_mid_temp = str2num(global_values{3})
global_high_temp = str2num(global_values{4})

line3 = fgetl(K);
line4 = fgetl(K);
line5 = fgetl(K);

for i = 1:53
    T = fgetl(K);
    species_name = T(1:9);
    species_name = deblank(species_name);
   
    species_temp = T(49:75);
    local_values = strsplit(species_temp,' ');
    local_low_temp = str2num(local_values{1});
    local_mid_temp = str2num(local_values{3});
    local_high_temp = str2num(local_values{2});
    
    temp = [local_low_temp local_mid_temp local_high_temp];
    temp_range = linspace(local_low_temp,local_high_temp,500);
    
    coeff = fgetl(K);
    coeff1 = findstr(coeff,'E');
    
    a1 = coeff(1:coeff1(1)+3);
    a1 = str2num(a1);
    a2 = coeff(coeff1(1)+4:coeff1(2)+3);
    a2 = str2num(a2);
    a3 = coeff(coeff1(2)+4:coeff1(3)+3);
    a3 = str2num(a3);
    a4 = coeff(coeff1(3)+4:coeff1(4)+3);
    a4 = str2num(a4);
    a5 = coeff(coeff1(4)+4:coeff1(5)+3);
    a5 = str2num(a5);

    new_line = fgetl(K);
    coeff2 = findstr(new_line,'E');
    
    a6 = new_line(1:coeff2(1)+3);
    a6 = str2num(a6);
    a7 = new_line(coeff2(1)+4:coeff2(2)+3);
    a7 = str2num(a7);
    a8 = new_line(coeff2(2)+4:coeff2(3)+3);
    a8 = str2num(a8);
    a9 = new_line(coeff2(3)+4:coeff2(4)+3);
    a9 = str2num(a9);
    a10 = new_line(coeff2(4)+4:coeff2(5)+3);
    a10 = str2num(a10);
    
    new_line1 = fgetl(K);
    coeff3 = findstr(new_line,'E');
    
    a11 = new_line(1:coeff3(1)+3);
    a11 = str2num(a11);
    a12 = new_line(coeff3(1)+4:coeff3(2)+3);
    a12 = str2num(a12);
    a13 = new_line(coeff3(2)+4:coeff3(3)+3);
    a13 = str2num(a13);
    a14 = new_line(coeff3(3)+4:coeff3(4)+3);
    a14 = str2num(a14);
    
    %Specific heat
    cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,temp_range,R);
    
    %enthalpy 
    H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,temp_range,R);
    
    %entropy
    S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R);
    
    %molecular weight
    m_weight = weight(species_name);
    
    %save molecular weight
    species_w = fopen('species_mass.txt','W');
    fprintf(species_w,'molecular mass of %s is %f',species_name,m_weight);
    fclose(species_w);
    
    %plot & save
    cur_dir = pwd;
    
    %create a seperate folder for each species
    mkdir(species_name);
    cd(species_name);
    
    figure(1)
    plot(temp_range,cp,'linewidth',2,'color',[0.2345 0.4567 0.5678])
    xlabel('Temperature(k)');
    ylabel('Specific heat(kJ/kg-k)');
    title(species_name);
    
    figure(2)
    plot(temp_range,H,'linewidth',2,'color',[0.5678 0.3423 0.7874])
    xlabel('Temperature(k)');
    ylabel('Enthalpy(kJ/kg)');
    title(species_name);    
    
    figure(3)
    plot(temp_range,S,'linewidth',2,'color',[0.4532 0.1234 0.3324])
    xlabel('Temperature(k)');
    ylabel('Entropy(kJ/k)');
    title(species_name);   
    
    %save the plots in seperate folder for each species
    saveas(1,'specific heat.jpg');
    saveas(2,'enthalpy.jpg');
    saveas(3,'entropy.jpg');
    
    cd(cur_dir);
    
    %printing the molecular weight of the species
    sprintf('molecular weight of %s = %f',species_name,m_weight)
    
end

%close the file
fclose(K)

function[Cp]= specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,temp,temp_range,R)
 %a1 a2 a3 a4 a5 are high temperature coefficients
 %a8 a9 a10 a11 a12 are low temperature coefficients
 for i= 1:length(temp_range)
 if temp_range(i)<=temp(2)
 Cp(i)= R*(a8 + (a9*temp_range(i)) + (a10*(temp_range(i).^2)) + (a11*(temp_range(i).^3)) + (a12*(temp_range(i).^4)));
 else
 Cp(i)= R*(a1 + (a2*temp_range(i)) + (a3*(temp_range(i).^2)) + (a4*(temp_range(i).^3)) + (a5*(temp_range(i).^4))); 
 
 end
 end

function[H]= enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,temp,temp_range,R)
 %a1 a2 a3 a4 a5 are high temperature coefficients
 %a8 a9 a10 a11 a12 are low temperature coefficients
 
 for i= 1:length(temp_range)
 if temp_range(i)<=temp(2)
 H(i)= R*((a8*temp_range(i)) + (a9*((temp_range(i).^2)/2)) + (a10*((temp_range(i).^3)/3)) + (a11*((temp_range(i).^4)/4)) + (a12*((temp_range(i).^5)/5)) + a13);
 
 else
 H(i)= R*((a1*temp_range(i)) + (a2*((temp_range(i).^2)/2)) + (a3*((temp_range(i).^3)/3)) + (a4*((temp_range(i).^4)/4)) + (a5*((temp_range(i).^5)/5)) + a6); 
 end
end

function[S]= entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,temp,temp_range,R)
 
 %a1 a2 a3 a4 a5 are high temperature coefficients
 %a8 a9 a10 a11 a12 are low temperature coefficients
 
 for i= 1:length(temp_range)
 if temp_range(i)<=temp(2)
 S(i)= R.*((a8*log(temp_range(i))) + (a9*temp_range(i)) + (a10*((temp_range(i).^2)/2)) + (a11*((temp_range(i).^3)/3)) + (a12*((temp_range(i).^4)/4)) + a14);
 
 else
 S(i)= R.*((a1*log(temp_range(i))) + (a2*temp_range(i)) + (a3*((temp_range(i).^2)/2)) + (a4*((temp_range(i).^3)/3)) + (a5*((temp_range(i).^4)/4)) + a7);
 
 end
end

function[m_weight]= weight(species_name)
 
 element= {'H','O','N','Ar','C','S'};
 atomic_mass= [1.008, 15.999, 14.007, 39.948, 12.011, 32.065];
 m_weight=0;
 for i= 1:length(species_name)
 for j= 1:length(element)
 if strcmp(species_name(i),element(j))
 m_weight= m_weight + atomic_mass(j);
 a= j;
 end
 
 end
 
 b= str2num(species_name(i));
 if b>1
 m_weight= m_weight + atomic_mass(a)*(b-1);
 end
 end
end