Project 1 - Parsing NASA thermodynamic data

Aim : To develop a program to parse data from Nasa Thermodynamic Data file and calculate enthaply, entropy, specific heat and molecular weight.

Objective :
1. To get co-efficients from Nasa thermodynamic data file and calculate the enthalpy, entropy and specific heat with the given formula.
2. Calulate the molecular weight for every species and display the eoutput.
3.Plot the graphs for enthalpy, entropy,specific heat vs temperature.

Formula Used :
• Cp = R( a₁ + a₂T + a₃T² +a₄T³ + a₅T⁴)
• H = RT( a₁ + ½a₂T + ⅓a₃T² + ¼a₄T³ + ⅕a₅T⁴ +a₆T )
• S = R( a₁ln(T) + ½a₂T + ⅓a₃T² + ¼a₄T³ + ¼a₅T⁴ + a₇ )
where,
Cp - Specific heat capacity
H - Enthalpy
S – Entropy
a₁, a₂, a₃, a₄, a₅, a₆, a₇ are the co-efficients from the thermo.dat file.
T – Temperature
R – Universal Gas Constant

Algorithm :
1. To get the species, temperature and co-efficients from the thermo.dat file.
2. Calculate enthaply with co-efficients,temperature and R as argument.
3. Calculate entropy with co-efficients,temperature and R as argument.
4. Calculate specific heat capacity with co-efficients,temperature and R as argument.
5. Create a function to calculate molecular weight of the species.
6. Display the output of molecular weights for different species.
7. Plot the graph for specific heat, entropy, enthaply vs temperature and save all the plot in the respective species files.

Syntax of the code used in the program :
File handle = fopen('filename','r'); - The function fopen opens the file ‘filename’ in the current directory with read mode and stores it in the file handle.
Line =fgetl(file handle); - The function fgetl reads the line one by one from the file handle
line=strsplit(line); - strsplit splits the function with the given identifier.
GLT=str2num(GLT); - str2num converts a string into a number.
E = strfind(line2,'E'); - strfind finds the given argument in a particular line.
cur_dir=pwd; - this sets the present working directory as current directory.
mkdir(mol); - this command makes a directory with given argument name in current directory.
saveas(figure(1),'Specific Heat.jpg');- this command saves the current figure as per the name given in the argument.

clear all
close all
clc

%getting file thermo.dat in read mode
f1 = fopen('THERMO.dat','r'); 
%getting line
l1=fgetl(f1);
l2=fgetl(f1);

%global temeprature
l2=strsplit(l2);    
GLT=str2num(l2{2});   %global low temperature
GMT=str2num(l2{3});   %global mid temperature
GHT=str2num(l2{4});   %global high temperature

%getting comment lines 
l3=fgetl(f1);
l4=fgetl(f1);
l5=fgetl(f1);

for i=1:53
    
  line1 = fgetl(f1);
  MOL = strsplit(line1,' ');
  mol=(MOL{1});
  
  %local temperature 
  llt=str2num(MOL{end-3});  %local low temperature
  hlt=str2num(MOL{end-2});  %local high temperature 
  mlt=str2num(MOL{end-1});  %local mid temperature
  
  %extracting co-eeficients
  line2=fgetl(f1);            
  E = strfind(line2,'E');
  a1 = line2(1:E(1)+3);
  a1=str2num(a1);
  a2 = line2(E(1)+4:E(2)+3);
  a2=str2num(a2);
  a3 = line2(E(2)+4:E(3)+3);
  a3=str2num(a3);
  a4 = line2(E(3)+4:E(4)+3);
  a4=str2num(a4);
  a5 = line2(E(4)+4:E(5)+3);
  a5=str2num(a5);
  
  line3=fgetl(f1);
  E2 = strfind(line3,'E');
  a6 = line3(1:E2(1)+3);
  a6=str2num(a6);
  a7 = line3(E2(1)+4:E2(2)+3);
  a7=str2num(a7);
  a8 = line3(E2(2)+4:E2(3)+3);
  a8=str2num(a8);
  a9 = line3(E2(3)+4:E2(4)+3);
  a9=str2num(a9);
  a10 = line3(E2(4)+4:E2(5)+3);
  a10=str2num(a10);
   
  line4=fgetl(f1);
  E3 = strfind(line4,'E');
  a11 = line4(1:E3(1)+3);        
  a11=str2num(a11);
  a12 = line4(E3(1)+4:E3(2)+3);
  a12=str2num(a12);
  a13 = line4(E3(2)+4:E3(3)+3);
  a13=str2num(a13);
  a14 = line4(E3(3)+4:E3(4)+3);
  a14=str2num(a14);
  
  R=8.314; %universal gas constant
  
  %temperature range 
  Tl=linspace(llt,mlt,500);
  Th=linspace(mlt,hlt,500);
  
  %properties for higher temperature range
  cph = (a1+ (a2*Th) + (a3*Th.^2) + (a4*Th.^3) + (a5*Th.^4))*R;
  hh = (a1+(a2*Th/2)+(a3*(Th.^2)/3)+(a4*(Th.^3)/4)+(a5*(Th.^4)/5)+(a6./Th))*R.*Th;
  sh = ((a1.*log(Th))+(a2*Th)+(a3*Th.^2/2)+(a4*Th.^3/3)+(a5*Th.^4/4)+a7)*R;
  %properties for lower temperature range
  cpl = (a8+ (a9*Tl) + (a10*Tl.^2) + (a11*Tl.^3) + (a12*Tl.^4))*R;
  hl = (a8+(a9*Tl/2)+(a10*(Tl.^2)/3)+(a11*(Tl.^3)/4)+(a12*(Tl.^4)/5)+(a13./Tl))*R.*Tl;        
  sl = ((a8.*log(Tl))+(a9*Tl)+(a10*Tl.^2/2)+(a11*Tl.^3/3)+(a12*Tl.^4/4)+a14)*R;

%weight output
w=calculate_weight(mol);
printf('WEIGHT OF %s IS : %d \n',mol,w);

%setting directory
cur_dir=pwd;
mkdir(mol);
cd(mol);

%temperature vs specific heat plot
figure(1);
plot(Tl,cpl,'linewidth',2);
hold on
plot(Th,cph,'linewidth',2);
xlabel('Temperature');
ylabel('Specicic heat');
legend('Low Temperature','High Temperature');
title(sprintf('Temperature vs Specific Heat for %s',mol));

hold off
saveas(figure(1),'Specific Heat.jpg');

%temperature vs entropy plot
figure(2);
plot(Th,sh,'linewidth',2);
hold on
plot(Tl,sl,'linewidth',2);
xlabel('Temperature');
ylabel('Entropy');
legend('Low Temperature','High Temperature');
title(sprintf('Temperature vs Entropy for %s',mol));
hold off
saveas(figure(2),'Entropy.jpg');

%temperature vs enthaply plot
figure(3);
plot(Tl,hl,'linewidth',2);
hold on
plot(Th,hh,'linewidth',2);
xlabel('Temperature');
ylabel('Enthalpy');
legend('Low Temperature','High Temperature');
title(sprintf('Temperature vs Enthaply for %s',mol));
hold off
saveas(figure(3),'Enthalpy.jpg');

cd(cur_dir);
  
end

%function for calculating the molecular weight
function w=calculate_weight(mol)

elements=['A','a','C','c','H','h','N','n','O','o'];
mass=[40 40 12 12 1 1 14 14 16 16];
w=0;

 for i=1:length(mol)  %run for length of the species
   for j=1:length(elements)  %run for length of elements
     if strcmp(mol(i),elements(j))  %comparing with every element
       w=w+mass(j);              
       position=j; 
     end
     end
     n=str2num(mol(i));  %check for multiple number
     if n>1
       w = w+mass(position)*(n-1);
     end  
   end

Command Window Output :

WEIGHT OF CH3 IS : 15
WEIGHT OF CH4 IS : 16
WEIGHT OF CO IS : 28
WEIGHT OF CO2 IS : 44
WEIGHT OF HCO IS : 29
WEIGHT OF CH2O IS : 30
WEIGHT OF CH2OH IS : 31
WEIGHT OF CH3O IS : 31
WEIGHT OF CH3OH IS : 32
WEIGHT OF C2H IS : 25
WEIGHT OF C2H2 IS : 26
WEIGHT OF C2H3 IS : 27
WEIGHT OF C2H4 IS : 28
WEIGHT OF C2H5 IS : 29
WEIGHT OF C2H6 IS : 30
WEIGHT OF CH2CO IS : 42
WEIGHT OF HCCO IS : 41
WEIGHT OF HCCOH IS : 42
WEIGHT OF H2CN IS : 28
WEIGHT OF HCN IS : 27
WEIGHT OF HNO IS : 31
WEIGHT OF N IS : 14
WEIGHT OF NNH IS : 29
WEIGHT OF N2O IS : 44
WEIGHT OF NH IS : 15
WEIGHT OF NH2 IS : 16
WEIGHT OF NH3 IS : 17
WEIGHT OF NO IS : 30
WEIGHT OF NO2 IS : 46
WEIGHT OF HCNO IS : 43
WEIGHT OF HOCN IS : 43
WEIGHT OF HNCO IS : 43
WEIGHT OF NCO IS : 42
WEIGHT OF CN IS : 26
WEIGHT OF HCNN IS : 41
WEIGHT OF N2 IS : 28
WEIGHT OF AR IS : 40
WEIGHT OF C3H8 IS : 44
WEIGHT OF C3H7 IS : 43
WEIGHT OF CH3CHO IS : 44
WEIGHT OF CH2CHO IS : 43

Conclusion :

The co-efficients from the Nasa thermodynamic file was extracted and the plots for enthaply entropy and specific heat capacity vs temperature were generated for the local temperature ranges and saved.

Errors Occured :

There were several syntax errors that occured while developing the code which were cleared with the matlab help library for functions.

A logical error occured with extraction which made the graph inconsistent, this which was resolved by correct assignment of the arguments.