function [error,torigin] = ...
         cal_errorinlo(stla,stlo,time,baz,evla,evlo,trefgvel)

% function designed to calculate error by comparing difference in time and 
% back azimuth.
% Preferred norm is selected by referring to Kennett(1996)
% Update; Normaization for epicentral distance. (08/27/2003)

Erad = 6371;

% Calculate origintime by average origin time
% (time is matched with julian days)

tmpepi0(1) = cal_epi1(stla(1),stlo(1),evla,evlo);
tmpepi0(2) = cal_epi1(stla(2),stlo(2),evla,evlo);
tmpepi0(3) = cal_epi1(stla(3),stlo(3),evla,evlo);

dtorigin = tmpepi0'*Erad*pi/180/trefgvel./(24*3600);
origin = time'-dtorigin;
torigin = mean(origin);

% Epicentral distance from observed time in km.
obsepi(1)=(time(1)-torigin)*24*3600*trefgvel;
obsepi(2)=(time(2)-torigin)*24*3600*trefgvel;
obsepi(3)=(time(3)-torigin)*24*3600*trefgvel;

% Epicentral distances from trial epicenter in km.
tmpepi=tmpepi0.*Erad*pi/180;

% Back azimuths from trial epicenter.
tmpbaz(1)=cal_baz(stla(1),stlo(1),evla,evlo);
tmpbaz(2)=cal_baz(stla(2),stlo(2),evla,evlo);
tmpbaz(3)=cal_baz(stla(3),stlo(3),evla,evlo);

% Calculate Error by combining difference in epicenter and backazimuth.

alpha = 0.8;
n1=3;
errepi=(1-alpha).*((obsepi'-tmpepi').*tmpepi0'./90).^2./n1;

tempin=zeros(size(baz,2),3);
out = zeros(size(baz,2),1);
tempin(:,1) = tmpbaz'-baz';
tempin(:,2) = -(tmpbaz'+360-baz');
tempin(:,3) = baz'+360-tmpbaz';

[dumm,imin]=min(abs(tempin'));

for k=1:size(baz,2)
    out(k,1)=tempin(k,imin(k));
end

n2=3;

ttang=tmpepi0'*pi/180;
errbaz=alpha.*(out.*pi/180.*1000.*sin(ttang)).^2./n2;
error=real(sum(sqrt((errepi+errbaz)).^(1.25)));
