%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 
%%% Implied volatility calculator
%%% 
%%% 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 [sigma_call, sigma_put] = ImpliedVolatility(S, K, T, r, q, P_call, P_put)

if ((length(S) ~= size(P_call,1)) || (length(K) ~= size(P_call,2)))
    error('size of K and P_call do not match');
end
if ((length(S) ~= size(P_put,1)) || (length(K) ~= size(P_put,2)))
    error('size of K and P_put do not match');
end

sigma_call = zeros(length(S),length(K));
sigma_put = zeros(length(S),length(K));

for iter_s = 1:length(S)
    for iter_k = 1:length(K)
        sigma_call(iter_s, iter_k) = blsimpv(S(iter_s),K(iter_k),r,T,P_call(iter_s, iter_k),[], q, [], {'Call'});
        sigma_put(iter_s, iter_k) = blsimpv(S(iter_s),K(iter_k),r,T,P_put(iter_s, iter_k),[], q, [], {'Put'});
    end
end