%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Option pricers driver for European option under Black-Scholes-Merton
%%% model
%%% 
%%% Copyright (C) 2008  Vladimir Surkov
%%%
%%%
%%% This program is free software: you can redistribute it and/or modify
%%% it under the terms of the GNU General Public License as published by
%%% the Free Software Foundation, either version 3 of the License, or
%%% (at your option) any later version.
%%%
%%% This program is distributed in the hope that it will be useful,
%%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%%% GNU General Public License for more details.
%%% 
%%% You should have received a copy of the GNU General Public License
%%% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Driver_European_BlackScholesMerton(S, K, T, sigma, r, q)

P = 4;
times = zeros(1,P+1);

[c1, p1, times(P+1)] = ClosedForm_European_BlackScholesMerton(S, K, T, sigma, r, q);

for pricer = 1:P
    if pricer == 1
        [c2, p2, times(pricer)] = ConvolutionQuadrature_European_BlackScholesMerton(S, K, T, sigma, r, q);
    elseif pricer == 2
        [c2, p2, times(pricer)] = FourierProductQuadrature_European_BlackScholesMerton(S, K, T, sigma, r, q);
    elseif pricer == 3
        [c2, p2, times(pricer)] = CarrMadanFFT_European_BlackScholesMerton(S, K, T, sigma, r, q);
    else
        [c2, p2, times(pricer)] = FST_European_BlackScholesMerton(S, K, T, sigma, r, q);
    end
        
    if ((length(S) > 1) && (length(K) > 1))
        subplot(P,2,pricer*2-1);
        surf(S, K, log(abs(c1-c2))');
        subplot(P,2,pricer*2);
        surf(S, K, log(abs(p1-p2))');
    elseif (length(S) > 1)
        subplot(P,2,pricer*2-1);
        plot(S, log(abs(c1-c2)));
        subplot(P,2,pricer*2);
        plot(S, log(abs(p1-p2)));
    elseif (length(K) > 1)
        subplot(P,2,pricer*2-1);
        plot(K, log(abs(c1-c2)));
        subplot(P,2,pricer*2);
        plot(K, log(abs(p1-p2)));
    else
        log(abs(c1-c2))
        log(abs(p1-p2))
    end
end

figure;
bar(times);