%   NymFys - Example 3.1 - Martian atmosphere
%   18 December 2012

clear all;

R=8.314;    % universal gas constant
g=3.71;     % gravitational constant on Mars (at ground)
M=0.04334;  % mean molecular weight near ground (on Mars) in kg/mole

z0=0.0;     % defining the interval (ground level)
z1=15.0;    % defining the interval (height; in km)
N=100;      % defining the amount of points for Simpson's rule
h=(z1-z0)/(N-1);    % defining the length of each interval

% computing the integral
F=0; 

for i=1:N
     z=z0+(i-1)*h;
     f=1/(234.0-2.25*z+14.0*exp(-2.0*z));  
     if (i==1) | (i==N)     % if we are at the left or right end of the main
                            % interval
        F=F+f;  
     else
        if mod(i,2)==0     % if we are at an "even" point
           F=F+4*f;     
        else               % if we are at an "odd" point
           F=F+2*f; 
        end;
     end;
 end;
 
integral=(1000.0*h/3)*(-g*M/R)*F % evaluating the integral 
                                 % (factor 1000 due to km to m conversion) 

% computing the change in pressure 
percentage=100*exp(integral) % pressure at z1 compared to z0 in per cent