%
% euler.m - MATLAB code
%
% Solve the "Why does not the atmosphere fall down"-problem
% from Lecture 1 in the Computational Physics course at
% Lund University, FYTN03. The Euler algorithm is used.
%
% To run the program simply type 'euler' at the MATLAB-prompt.
%
% by Tobias Ambjornsson, 2010
%

clear all;      
format long;    

% --- Input parameters ---
y0=18000;             % characteristic length, units of m
h=0.1*y0;             % step length, units of m 
rho0=1.0;             % density at ground level, units of atmospheres 
N=80;                 % number of euler steps

% --- Initialize vectors ---  
rho=zeros(N,1);       % initialize vector containing rho
yvec=zeros(N,1);      % initialize vector containing heights
%

% --- The Euler algorithm ---
rho(1)=rho0;          % density at the ground
yvec(1)=0;            % height at the ground = 0
for n=1:1:N-1         % loop over different heights
   yvec(n+1)=yvec(n)+h;
   rho(n+1)=rho(n)-rho(n)*h/y0;    
end;

% --- Plot ---
set(gca,'Fontsize',15); 
plot(yvec,rho,'bo','Linewidth',2); hold on;  % plot with marks
plot(yvec,rho,'b-','Linewidth',2); hold on;  % plot a continuous curve
xlabel('y (m)','Fontsize',15);               % set labels on x-axis
ylabel('\rho (atmospheres)','Fontsize',15);  % set label on y-axis
hold off;

