Blog: Stable Atmosphere: process_lineout.m

clear

syms r

rCrit = 3.587091;
rZero = 1.386534;
rOb = 1.279937;
gamma = 1.66666667;
dx = 0.03125;

%Read in pressure and density information
%%%%%%
density_file = fopen('density_equilibrium.vtk');
text = textscan(density_file,'%s');
text = text{1};

index = find(strcmp(text,'POINTS')) + 1;
npoints = str2double(text{index});
index = find(strcmp(text,'double')) + 1;

for point = 1:npoints
   y(point) = str2double(text{index});
   density(point) = str2double(text{index+1});
   index = index + 3;
end

pressure_file = fopen('pressure_equilibrium.vtk');
text = textscan(pressure_file,'%s');
text = text{1};

index = find(strcmp(text,'POINTS')) + 1;
npoints = str2double(text{index});
index = find(strcmp(text,'double')) + 1;

for point = 1:npoints
   pressure(point) = str2double(text{index+1});
   index = index + 3;
end
%%%%%%%

for point = 1:npoints
    r = y(point);
    analytic_density(point) = real((rCrit*(1./r - 1./rZero))^(1./(gamma-1.)));
    analytic_pressure(point) = real((rCrit*(1./r - 1./rZero))^(gamma/(gamma-1.)));
end

% Plot density and pressure on separate graphs (can insert more plots here, change to log, etc.)
plot(y,density,'k','LineWidth',2)
hold on;
plot(y,analytic_density,'--')
xlabel('y (CU)')
ylabel('Density (CU)')

figure

semilogy(y,pressure,'k','LineWidth',2)
hold on;
semilogy(y,analytic_pressure,'--')
xlabel('y (CU)')
ylabel('Pressure (CU)')