How to calculate the integral of two areas substracted of each other

Hi,

My question is this. I have a PV-diagram with two areas above each other. How can I calculate the upper area minus the bottom area?

This is the code:

--------------------------------------------------------------------------------------------------- %% % variables b = 68 * 10^-3; % bore [m] s = 54 * 10^-3; % stroke [m] r = s / 2; % crank radius [m] l = 85 * 10^-3; % connecting rod length [m] estimated to be 85 mm CR = 8.5; % compression ratio [-], 8.5:1 V_b = s * (pi * b^2) / 4; % bore volume [m^3] V_c = V_b / (CR - 1); % clearance volume [m^3] TDC_shift = 0.11105; % indices shift to place ref on TDC in fraction of revolution time [C , D] = butter(2, 0.01); % smoothing settings P_ref = 1; % Reference pressure (bar) at 145 CA bTDC for pegging

% select and load one or more text files with format [time, pulse, pressure] files = fileselect;

    h1 = waitbar(0,['Processing 1/',num2str(size(files,2))]);

% convert the data out of the file in to struct per measurement, data(1) is the first measurement, data(2) the second... for i = 1:size(files,2)

    waitbar(i/(size(files,2)),h1,['Processing ',num2str(i),'/',num2str(size(files,2))])

    data(i).file = files{i};

    [data(i).time, data(i).pulse, data(i).pressure] = importfile(files{i});

    % store raw data
    raw(i).time=data(i).time;
    raw(i).pulse=data(i).pulse;
    raw(i).pressure=data(i).pressure;

    % shift data.pulse to a be balanced around 0
    data(i).pulse = data(i).pulse - (max(data(i).pulse)+min(data(i).pulse))/2;
    % determine roughly zero crossing of pulse
    ind1 = sort( find( (data(i).pulse(1:end-1) .* data(i).pulse(2:end)) < 0 ));
    % write zero crossing time, data.crs
    data(i).crs = data(i).time(ind1);
    % smooth pressure data
    data(i).pressure = filtfilt(C,D,data(i).pressure);
    % convert pressure data from V to Bar
    data(i).pressure = ((((data(i).pressure/5)*100)-11.5)*(200/80));
        % Nominaal is de uitgangsspanning 11.5% van de Vs (Voedingsspanning = 5V)
        % De factor (200/80) is de omrekenfactor volgens specificaties van de
        % fabrikant ( per 80% verhoging van de nominaal = 200 bar)

% % linear interpolation of crossing to get the exact time of the zero crossing, data.crse % data(i).crse = data(i).crs; % for j=1:length(data(i).crse) % if abs(S(ind(j))) > eps(S(ind(j))) % % interpolate only when data point is not already zero % NUM = (data(i).time(ind(j)+1) - data(i).time(ind(j))); % DEN = (data(i).pulse(ind(j)+1) - data(i).pulse(ind(j))); % data(i).crse(j) = data(i).crse(j) - data(i).pulse(ind(j)) * NUM / DEN; % end % end

    % determin the ref point
    ind2 = sort(find ( (data(i).crs(2:end-1) - data(i).crs(1:end-2)) ...
                        > (mean(data(i).crs(2:end-1) - data(i).crs(1:end-2))*1.5) ... 
                     & (data(i).crs(3:end)   - data(i).crs(2:end-1)) ...
                        < (mean(data(i).crs(3:end)   - data(i).crs(2:end-1))*1.5) ) );

    % determine the ref time, data.ref
    data(i).ref = data(i).crs(ind2);

    % determine indices in time vector of the ref points
    [~, data(i).indtime] = ismember(data(i).ref, data(i).time);

    % determine engine revolutions per min        
    for k = 1:size(data(i).ref,1) - 1
         data(i).rpm(k,:) = 60 / (data(i).ref(k+1) - data(i).ref(k));
    end

    % determine mean engine revolutions per min
    data(i).meanrpm = mean(data(i).rpm);

    % determine volume as function of crank angle per full revolution, data.rev.vol
    for m = 1:size(data(i).indtime,1) - 1
        for n = 1:data(i).indtime(m+1) - data(i).indtime(m)
            % determine crank angel for every data point between two ref points
            theta(m,n) = (n-1) * (2 * pi) / (data(i).indtime(m+1) - data(i).indtime(m)-1);
            x(n)     = r .* cos(theta(m,n)) + sqrt(l^2 - r^2 .* sin(theta(m,n)).^2);
            data(i).rev(m).vol(n,:) = ((r + l) - x(n)) * (pi * b^2 / 4) +V_c;
        end

        % determine shift of pressure 
        ind_shift = round( TDC_shift * (data(i).indtime(m+1)-data(i).indtime(m)));
        data(i).ind_shift_plot(m,:) = data(i).indtime(m) + ind_shift;

        % shift pressure to match volume, data.rev.pressure
        data(i).rev(m).pressure = data(i).pressure(data(i).indtime(m)+ind_shift:data(i).indtime(m+1)-1+ind_shift);
    end

    % Pegging: P_ave is used to determine the average of measured pressure at the reference point (here 145 bTDC)
    P_ave = 0;
    N = 0;
    for m = 1:size(data(i).indtime,1) - 1
        for n = 1:data(i).indtime(m+1) - data(i).indtime(m)
            if  theta(m,n) > 3.67 && theta(m,n) < 3.84      %  3.75 rad corresponds to 145 bTDC (3.67-> 150, 3.84->140 bTDC)
                P_ave = data(i).rev(m).pressure(n,1)+P_ave;
                N = N + 1;
            end
        end
    end   
    P_ave = P_ave/N;  % define the avarage of all the cycles
    data(i).pressure = data(i).pressure- P_ave + P_ref;   % correct the original values with the pegging values. 

        % plot P-V diagram of all revolutions on top of each other
    for m = 1:size(data(i).indtime,1) - 1
        data(i).rev(m).pressure = data(i).rev(m).pressure - P_ave + P_ref;
        figure(i*10)
        set(i*10,'units','pixel','Position',[scrsz]);
        hold on
        plot(data(i).rev(m).vol,data(i).rev(m).pressure,'Color','blue','LineWidth',1)
        title('P-V Diagram')
        xlabel('volume [m^{3}]')
        ylabel('pressure [Bar]')
    end

    % combine an even and an odd full revolution to get a complete cycle 
    % for the pressure and the volume, data.cyl.pressure and data.cyl.vol
    if mod(size(data(i).rev,2),2) == 0
        p = size(data(i).rev,2)/2;
    else
        p = (size(data(i).rev,2)-1)/2;
    end 

    for q = 1:p
        data(i).cyc(q).vol       = [data(i).rev(q*2-1).vol;      data(i).rev(q*2).vol];
        data(i).cyc(q).pressure  = [data(i).rev(q*2-1).pressure; data(i).rev(q*2).pressure];
        cycle_length(q)          = length(data(i).cyc(q).vol);
        data(i).cyc(q).cycle_ind = (1:length(data(i).cyc(q).vol))';
    end

    % average of volume and pressure over the cylces
    % find mean, round lenght of cycles
    mean_cycle_length = round(mean(cycle_length));

    % fit smoothingspline through volume and pressure and normalize over 
    % cycle length, data.cyc.fit_vol and data.cyc.fit_pressure
    for w = 1:p
        data(i).cyc(w).fit_vol = fit(data(i).cyc(w).cycle_ind/cycle_length(w), data(i).cyc(w).vol, 'smoothingspline');
        data(i).cyc(w).fit_pressure = fit(data(i).cyc(w).cycle_ind/cycle_length(w), data(i).cyc(w).pressure, 'smoothingspline');
    end

    % calculate volume and pressure at the same normalized point for all
    % cycles and calculate the mean, data.mean_vol and data.mean_pressure
    vol = []; pressure = [];
    for g = 1 : mean_cycle_length
        for h = 1:p
            vol(h) = feval(data(i).cyc(h).fit_vol,g/mean_cycle_length);
            pressure(h) = feval(data(i).cyc(h).fit_pressure,g/mean_cycle_length);
        end
        data(i).mean_vol(g)         = mean(vol);
        data(i).mean_pressure(g)    = mean(pressure);
    end

% save data, data.mean_vol and data.mean_pressure mean_vol = data(i).mean_vol; mean_pressure = data(i).mean_pressure; % save(['results_',data(i).file],'mean_vol','mean_pressure','-ascii','-tabs')

-------------------------------------------------------------------------------------

My guess is it has something to do with "strapz", but this only gives me the two areas summed together.