function [sigma_min_dirac,sigma_max_dirac,epsilon] = Dirac_delta_bound(S0,K,r,q,tau,C_real,N,N2,sigma_down,sigma_up,epsilon0,k1)
% N: the number of abscissas for calculating epsilon
% N2: the number of abscissas for calculating boundary
% sigma_down,sigma_up: naive boundary
% epsilon0: naive limiting parameter
% k1 = 1e-3: adjusted coefficient related to boundary

% generate discrete points
sigma_all = (sigma_down:(sigma_up - sigma_down)/N2:sigma_up)';
C_all = European_call(S0,K,r,tau,sigma_all,q);

% locate the positions to compute epsilon according to formula 11
temp00 = 1:round(N2/N):N2;
temp0 = sum(C_all(temp00,:)<C_real);
temp1 = [temp0.*4;temp0.*4+1];

% compute the vega and its derivative
vega1 = (exp((-((tau.^2.*(4.*(q-r).^2+4.*(q+r).*(sigma_all(temp1)).^2+(sigma_all(temp1)).^4)+...
    4.*log(S0./K).^2)./(8.*tau.*(sigma_all(temp1)).^2)))).*K.*(S0./K).^(1./2+(q-r)./(sigma_all(temp1)).^2).*...
    sqrt(tau))./sqrt(2.*(pi));
vvega1 = -(1./(4.*sqrt(2.*(pi)).*sqrt(tau).*(sigma_all(temp1)).^3)).*exp((-((tau.^2.*...
    (4.*(q-r).^2+4.*(q+r).*(sigma_all(temp1)).^2+(sigma_all(temp1)).^4)+4.*log(S0./K).^2)./(8.*tau.*...
    (sigma_all(temp1)).^2)))).*K.*(S0./K).^(1./2+(q-r)./(sigma_all(temp1)).^2).*(tau.*(2.*q-2.*r+(sigma_all(temp1)).^2)-...
    2.*log(S0./K)).*(tau.*(-2.*q+2.*r+(sigma_all(temp1)).^2)+2.*log(S0./K));
vvvega1 = (1./(16.*sqrt(2.*(pi)).*tau.^(3./2).*(sigma_all(temp1)).^6)).*exp(-((tau.^2.*(4.*(q-r).^2+4.*...
    (q+r).*(sigma_all(temp1)).^2+(sigma_all(temp1)).^4)-4.*log(S0./K).^2)./(8.*tau.*(sigma_all(temp1)).^2))).*K.*(S0./K).^(1./2+...
    q./(sigma_all(temp1)).^2-(r.*tau+log(S0./K))./(tau.*(sigma_all(temp1)).^2)).*(tau.^3.*(16.*(q-r).^4.*tau-48.*(q-...
    r).^2.*(sigma_all(temp1)).^2-8.*(q-r).^2.*tau.*(sigma_all(temp1)).^4-4.*(sigma_all(temp1)).^6+tau.*(sigma_all(temp1)).^8)-8.*(2.*(q-r).*...
    tau-log(S0./K)).*log(S0./K).*(tau.*(4.*(q-r).^2.*tau-6.*(sigma_all(temp1)).^2-tau.*(sigma_all(temp1)).^4)+2.*...
    log(S0./K).*(2.*(-q+r).*tau+log(S0./K))));

% compute the reasonable epsilon according to formula 13
epsilon = min(min(abs(sigma_all(temp1).*vega1.^3./(4.*vvega1 + ...
    2.*sigma_all(temp1).*vvvega1)).*k1),epsilon0);

% recompute the vega and its derivative
vega_all = (exp((-((tau.^2.*(4.*(q-r).^2+4.*(q+r).*(sigma_all).^2+(sigma_all).^4)+...
    4.*log(S0./K).^2)./(8.*tau.*(sigma_all).^2)))).*K.*(S0./K).^(1./2+(q-r)./(sigma_all).^2).*...
    sqrt(tau))./sqrt(2.*(pi));
vvega_all = -(1./(4.*sqrt(2.*(pi)).*sqrt(tau).*(sigma_all).^3)).*exp((-((tau.^2.*...
    (4.*(q-r).^2+4.*(q+r).*(sigma_all).^2+(sigma_all).^4)+4.*log(S0./K).^2)./(8.*tau.*...
    (sigma_all).^2)))).*K.*(S0./K).^(1./2+(q-r)./(sigma_all).^2).*(tau.*(2.*q-2.*r+(sigma_all).^2)-...
    2.*log(S0./K)).*(tau.*(-2.*q+2.*r+(sigma_all).^2)+2.*log(S0./K));

% compute the intergrand
f1 = 1./sqrt(2.*pi.*epsilon).*exp(-(C_all - C_real).^2./(4.*epsilon)).*vega_all + ...
    sigma_all./sqrt(2.*pi.*epsilon).*exp(-(C_all - C_real).^2./(4.*epsilon)).*vega_all.*...
    (-2.*(C_all - C_real)./(4.*epsilon).*vega_all) + ...
    sigma_all./sqrt(2.*pi.*epsilon).*exp(-(C_all - C_real).^2./(4.*epsilon)).*vvega_all;

% compute Dirac delta bound according to formula 14
% the abscissas corresponding to the maximum and minimum values
[temp4,~] = find(f1==max(f1));
sigma_min_dirac  = sigma_all(temp4)';
[temp4,~] = find(f1==min(f1));
sigma_max_dirac  = sigma_all(temp4)';



end