function [sigma] = BSMiVolCall(S, E, r, tau, Cgiven)

% [sigma] = BSMCallIVol(S,E,r,tau,Cgiven)
% finds implied volatility of a call option by applying Newton's method
%

assert(tau > 0, 'Time to expiration must be positive');

%starting value
sigma = sqrt(2*abs( (log(S/E) + r*tau)/tau ) );

% desired precision
tol = 1e-8;

% The function N
N = @(x) (1+erf(x/sqrt(2))) / 2;
dN = @(x) 1/sqrt(2*pi) * exp(-x.^2/2);

% Solving equation C(..., sigma) - Cgiven = 0
for k = 1:100,
    % Calculate approx Call value and vega
    d1 = (log(S/E) + (r + 0.5*sigma^2)*(tau))/(sigma*sqrt(tau));
    d2 = d1 - sigma*sqrt(tau);
    C = S.*N(d1) - E*exp(-r*(tau))*N(d2);
    Cvega = S*sqrt(tau)*dN(d1);
    
    increment = (C-Cgiven)/Cvega;
    sigma = sigma - increment;
    if abs(increment) < tol,
        % Method converged
        return; 
    end;
end;

% Did not converge
sigma = -1;

end