function [delta, gamma, theta, rho, vega] = eurcallgreeks(S0, K, r, T, t, sigma, D0)
% usage:  [delta, gamma, theta, rho, vega] = eurcallgreeks(S0, K, r, T, t, sigma, D0)
% description: Returns the greek sensitivities for a (vector of) European
%     call options with continuous dividends
% parameters:
%     S0 - vector of initial prices
%     K - strike price
%     r - risk-free return rate r (in decimal format)
%     T - time of the call
%     t - current time
%     sigma - standard deviation
%     D0 - continuously compounded dividend rate
% returns: sensitivities to changes in the portfolio
%     delta - changes in stock price
%     gamma - higher order changes in stock price
%     theta - changes in time
%     rho - changes in interest rate
%     vega - changes in volatility

% Contact: Brian Holley, bholley@unlnotes.unl.edu

% local variables:
%     d1, d2 - helper variables for Black-Scholes
%     cdfd1, cdfd2 - normcdf values at d1 and d2

% verify parameters
if (nargin ~= 7)
  error('eurcallgreeks: need 7 arguments: S0, K, r, T, t, sigma, D0');
end
if (length(S0) < 1)
    error('eurcallgreeks: S0 must have at least 1 entry');
end
if (K < 0)
    error('eurcallgreeks: K must be positive');
end
if (r < 0)
    error('eurcallgreeks: r must be positive');
end
if (r > 1)
    warning('eurcallgreeks: r may not be in decimal form');
end
if (T < 0)
    error('eurcallgreeks: T must be positive');
end
if (t < 0)
    error('eurcallgreeks: t must be positive');
end
if (D0 < 0)
    error('eurcallgreeks: D0 must be positive');
end

% start with some helper variables
d1 = (log(S0./K) + (r - D0 + (sigma ^ 2) / 2) * (T-t)) / (sigma * sqrt(T-t));
d2 = (log(S0./K) + (r - D0 - (sigma ^ 2) / 2) * (T-t)) / (sigma * sqrt(T-t));
cdfd1 = 0.5 * erfc(-d1 / sqrt(2));
cdfd2 = 0.5 * erfc(-d2 / sqrt(2));

% begin computing the greeks
% formulas used from http://www.riskglossary.com/articles/merton_1973.htm
% later found to be incorrect(?), adapting from Prof. Steve Dunbar's notes
% available at http://www.math.unl.edu/~sdunbar/Teaching/MathematicalFinance/Lessons/BlackScholes/Greeks/greeks.shtml
delta = (exp(-D0 * (T-t)) * cdfd1)';
%gamma = ((cdfd1 .* exp(-D0 * (T-t))) ./ (S0 .* sigma .* sqrt(T-t)))';
gamma = (1 ./ (S0 .* sqrt(2 * pi) .* sigma .* sqrt(T-t)) .* exp(-d1 .^ 2 / 2))';
%theta = (-(S0 .* exp(-D0 * (T-t)) .* cdfd1 .* sigma) ./ (2 * sqrt(T-t)) + (D0 .* S0 .* exp(-D0 * (T-t)) .* cdfd1) - (r .* K .* exp(-r * (T-t)) .* cdfd2))';
theta = (-(S0 .* exp(-d1 .^ 2 / 2) .* sigma) / (2 .* sqrt(2 * pi) .* sqrt(T-t)) - (r-D0) .* K .* exp(-(r-D0) * (T-t)) .* cdfd2)';
rho = (K .* (T-t) .* exp(-r * (T-t)) .* cdfd2)';
%vega = (S0 .* exp(-D0 * (T-t)) .* cdfd1 .* sqrt(T-t))';
vega = ((S0 .* sqrt(T-t)) ./ sqrt(2 * pi) .* exp(-d1 .^ 2 / 2))';