function curiv= euriv(cp,k,s0,tau,rfr,coc,niters)
%  curiv= euriv(cp,k,s0,tau,rfr,coc,niters)
%
%
%   procedure to calculate implied volatilities of a set of European option 
%   quotes using a Newton Raphson type Method.
%
%   See AMERIV for American option equivalent
%
%
%Format: iv=EURIV(cp,k,s0,tau,rfr,coc,niters)
%
% Where: iv  is a Nx1 vector or scaler of implied volatilities
%
%        cp  is a Nx1 vector or scaler of call prices (or put prices)
%
%        k   is a Nx1 vector or scaler of strike prices (negative of strike
%            for puts)
%
%        s0  is a Nx1 vector or scaler for current level of the underlying
%            asset
%
%        tau is a Nx1 vector or sclaer of time to exp. in yrs
%
%        rfr is a Nx1 vector or scaler for the continuously compounded
%            annualized risk free rate
%            Note: if you have an effective annual yield (annyield) simply
%            pass ln(1+annyield) to the procedure.
%
%        coc is the annualized rate for cost-to-carry
%            for non-div. paying stock this is simply the rfr
%            for div paying stock this is rfr-divrate
%            for currency options this is rdom-rfor
%
%     niters is a scaler indicating the number Newton-Raphson iterations
%            to conduct (3 or 4 yields reasonable results)
%
%      note: You may want to order the command disable before calling this
%            function, for particularly wild data. Be sure to issue the
%            command enable afterwards, to reactivate error checking.
%
%
rcp=rows(cp); rs0=rows(s0);
rs=rows(k); rt=rows(tau); rr=rows(rfr); rc=rows(coc);
rowvec=([rcp;rs0;rs;rt;rr;rc]);
n=max(rowvec);
if n>1 & min(rowvec)~= n;
    if rs0==1;
      s0=s0.*ones(n,1);
   end;
    if rs==1;
      k=k.*ones(n,1);
   end;
    if rt==1;
      tau=tau.*ones(n,1);
   end;
    if rr==1;
      rfr=rfr.*ones(n,1);
   end;
    if rc==1;
      coc=coc.*ones(n,1);
   end;
    if rcp==1;
      cp=cp.*ones(n,1);
   end;
end;  curiv=0.2.*ones(rows(cp),1);
% initial guess 
 ep=10^(-10);
 ep2=10^(-200);
iter=1;     while iter <= niters;

 [cfit]=bsval(curiv,k,s0,tau,rfr,coc);
 [cfit2]=bsval(curiv+ep,k,s0,tau,rfr,coc);
 vega=(cfit2-cfit)./ep;
 curiv=curiv+(cp-cfit)./(vega+ep2);
 iter=iter+1;
end;

% End of function