%Black-Scholes formula for the standard European call and put
%   *****************************************************
function price = bsoption(S,K,r,del,sigma,texp,c_p)
%   *****************************************************
%S: current stock price, 
%K: strike, 
%r: risk-free interrest rate,
%del: devidend rate,
%sigma: stock volatility
%texp: time to expiration (in years)
%   *****************************************************
eps = 1e-10;
if S/K < eps || sigma*sqrt(texp) < eps
switch c_p
    case 'call'
        price = 0.5*(S-K + abs(S-K));
    case 'put'
        price = 0.5*(K-S + abs(K-S));
    case 'dcall'
        price = heaviside(S-K);
    case 'dput'
        price = heaviside(K-S);
    otherwise
        error('c_p must be either call or put');
end
return
end

d1 = (log(S/K) + (r-del+sigma^2/2)*texp)/(sigma*sqrt(texp));
d2 = d1 - sigma*sqrt(texp);

switch c_p
    case 'call'
        price = 0.5*S*exp(-del*texp)*erfc(-d1/sqrt(2)) - ...
                0.5*K*exp(-r*texp)*erfc(-d2/sqrt(2));
    case 'put'
        price = 0.5*K*exp(-r*texp)*erfc(d2/sqrt(2)) - ...
                0.5*S*exp(-del*texp)*erfc(d1/sqrt(2)) ;
    case 'dcall'
        price = 0.5*exp(-r*texp)*erfc(-d2/sqrt(2));
    case 'dput'
        price = 0.5*exp(-r*texp)*erfc(d2/sqrt(2));
    otherwise
        error('c_p must be either call or put');
end
