Project 1 - Parsing NASA thermodynamic data
% NASA Thermodynamic Data Processing in MATLAB
clc; clear; close all;
% Open the file
file_path = 'THERMO.dat'; % Ensure this file is in the same directory
fid = fopen(file_path, 'r');
% Check if file opened successfully
if fid == -1
error('File not found or cannot be opened. Check if THERMO.dat is in the same directory.');
end
disp('File opened successfully!');
% Skip the first 5 header lines
for i = 1:5
fgetl(fid);
end
R = 8.314; % Universal gas constant (J/mol·K)
% Process each species
while ~feof(fid)
line = fgetl(fid);
if ~ischar(line) || isempty(line)
break;
end
% Extract species name
species = strtrim(line(1:18)); % First 18 characters
% Extract temperature limits
temp_data = str2num(line(46:end));
if length(temp_data) < 3
continue; % Skip if not enough temperature data
end
T_low = temp_data(1);
T_mid = temp_data(2);
T_high = temp_data(3);
% Read high-temperature coefficients
high_coeffs = [];
for i = 1:3
high_coeffs = [high_coeffs str2num(fgetl(fid))]; %#ok<AGROW>
end
high_coeffs = high_coeffs(1:7); % Ensure only 7 coefficients are used
% Read low-temperature coefficients
low_coeffs = [];
for i = 1:3
low_coeffs = [low_coeffs str2num(fgetl(fid))]; %#ok<AGROW>
end
low_coeffs = low_coeffs(1:7); % Ensure only 7 coefficients are used
% Check if coefficients are valid
if length(high_coeffs) < 7 || length(low_coeffs) < 7
continue; % Skip to next species if data is incomplete
end
% Temperature range for calculations
T = linspace(T_low, T_high, 100);
% Initialize Cp, H, S arrays
Cp = zeros(size(T));
H = zeros(size(T));
S = zeros(size(T));
% Compute Cp, H, S using NASA polynomials
for i = 1:length(T)
if T(i) >= T_mid
coeffs = high_coeffs;
else
coeffs = low_coeffs;
end
% Specific Heat Capacity (Cp)
Cp(i) = R * (coeffs(1) + coeffs(2)*T(i) + coeffs(3)*T(i)^2 + ...
coeffs(4)*T(i)^3 + coeffs(5)*T(i)^4);
% Enthalpy (H)
H(i) = R * T(i) * (coeffs(1) + coeffs(2)*T(i)/2 + coeffs(3)*T(i)^2/3 + ...
coeffs(4)*T(i)^3/4 + coeffs(5)*T(i)^4/5 + coeffs(6)/T(i));
% Entropy (S)
S(i) = R * (coeffs(1)*log(T(i)) + coeffs(2)*T(i) + coeffs(3)*T(i)^2/2 + ...
coeffs(4)*T(i)^3/3 + coeffs(5)*T(i)^4/4 + coeffs(7));
end
% Plot results
figure;
subplot(3,1,1);
plot(T, Cp, 'b', 'LineWidth', 2);
xlabel('Temperature (K)');
ylabel('Cp (J/mol·K)');
title(['Specific Heat Capacity of ', species]);
subplot(3,1,2);
plot(T, H, 'r', 'LineWidth', 2);
xlabel('Temperature (K)');
ylabel('Enthalpy (J/mol)');
title(['Enthalpy of ', species]);
subplot(3,1,3);
plot(T, S, 'g', 'LineWidth', 2);
xlabel('Temperature (K)');
ylabel('Entropy (J/mol·K)');
title(['Entropy of ', species]);
% Save plots in a folder
folder_name = ['Results_', species];
if ~exist(folder_name, 'dir')
mkdir(folder_name);
end
saveas(gcf, fullfile(folder_name, [species, '_plots.png']));
disp(['Plots saved in folder: ', folder_name]);
% Pause to allow visualization
pause(1);
end
% Close the file
fclose(fid);
disp('Processing complete!');