function cs= bsval(iv,k,s0,tau,rfr,coc)
%  cs= bsval(iv,k,s0,tau,rfr,coc)
%
%  Add brief description for function in line above for Matlab help searches
% BSVAL.G  Obtain European option, Black scholes values based on known
%            (possibly smile dependent) implied volatilities
%
%Format: value=bsval(iv,k,s0,tau,rfr,coc)
%
% Where: value is a Nx1 vector of Black Scholes European prices
%        iv  is a Nx1 vector of (implied) volatilities
%        k   is a Nx1 vector of strike prices (negative for puts)
%        tau is a 1x1 scaler or Nx1 vector of time to exp. in yrs
%        rfr is the continuously compounded annual risk free rate
%        coc is the cost to carry
%
%        Note: coc=rfr for a non-dividend paying stock
%              coc=rfr-anndivyield for dividend paying stock
%              coc=rdom-rfor for currency options
%              coc=0 for options on futures
%
%
%




n=max([rows(iv);rows(k);rows(s0);rows(tau);rows(rfr);rows(coc)]);
if rows(s0)==1;
   s0=s0.*ones(n,1);
end;
if rows(iv)==1;
   iv=iv.*ones(n,1);
end;
if rows(tau)==1;
   tau=tau.*ones(n,1);
end;
if rows(rfr)==1;
   rfr=rfr.*ones(n,1);
end;
  if rows(coc)==1;
   coc=coc.*ones(n,1);
 end;
 calldum=k>0;
 ep=10^(-10);
                    d1=(log(abs((s0+ep)./(abs(k)+ep)+ep))+(coc+(iv.^2)...
   ./2).*tau)./(iv.*sqrt(tau)+ep);

 d2=d1-iv.*sqrt(tau);
     cs=calldum.*(s0.*exp((coc-rfr).*tau).*cdfn(real(d1))-(k).*(exp(...
   -rfr.*tau)).*cdfn(real(d2)));

                    cs=cs+(1-calldum).*((abs(k)).*exp(-rfr.*tau).*cdfn...
   (-real(d2))-s0.*exp((coc-rfr).*tau).*cdfn(-real(d1)));


% End of function