function eqlocation1(N,plots)
% EQLOCATION1(N,plots)
%
% First find N random station locations
% Pick random points on a sphere
% Not just phi uniform on 0->2pi and theta uniform on 0->pi
% If that was the case, the probability of falling in patch
% tended by dtheta by dphi is sin(theta)dtheta*dphi
% (infinitesimal solid angle, area element...).
% So the smaller theta, the larger the probability, which leads
% to a clustering at the poles.
% But if pick theta=acos(1-2u) where u is uniform on 0->1
% then dtheta=-1/(sqrt(1-u^2))du, 
% sin(theta)=sin(acos(1-2u))=sqrt(1-u^2) so
% sin(theta)*dthat=-du
% hence probability of falling in this patch is -du*dphi/4pi
% which is clearly uniform.

defval('N',30);
defval('plots',0);
defval('Pspeed',8);

clf
% Start from slightly more favorable conditions
Slon=rand(N,1)*100;
Slat=rand(N,1)*45;
Tloc=randsph(1);
Gloc=randsph(1);
% Random origin time shift
Ttime=rand(1)*3000;


subplot(211)
plot(Slon,Slat,'+')
Glon=mean(Slon);
Glat=mean(Slat);
GTime=0;
hold on
Tlon=Glon-5;
Tlat=Glat-5;
b=plot(Tlon,Tlat,'o','MarkerS',6,'MarkerF','r');
axis equal

% Calculate true travel time from event to stations
[d,azt,tobs]=distat([Tlon Tlat],[Slon Slat],Pspeed,Ttime);
tobs=tobs;

% Build in convergence criterion
index=0;
tpred=1e99;
while mean(abs(tpred-tobs))>1
  index=index+1;
  disp(sprintf('Iteration %i',index))
  % Calculate distance and azimuths from the guessed location to the stations
  [delta,az,tpred]=distat([Glon Glat],[Slon Slat],Pspeed,GTime);
  % In first iteration, fix origin time
  GTime=GTime+mean(tobs-tpred);
  tpred=tpred+GTime*(index==1);
  delres=(tpred-tobs)*Pspeed;

  subplot(211)
  plot(Glon,Glat,'s','MarkerS',6,'MarkerF','g')
  subplot(212)
  sortplot(az,delres)
  [md,mdi]=max(delres);
  % Shift Glon and Glat by md in the direction az(mdi)  
  [Glon,Glat]=grcazim([Glon Glat],md,360-az(mdi));
  pause
end

% Plot residuals
subplot(211)
f=plot(Glon,Glat,'s','MarkerS',10,'MarkerF','g');
b=plot(Tlon,Tlat,'o','MarkerS',8,'MarkerF','r');
subplot(212)
[delta,az,tpred]=distat([Glon Glat],[Slon Slat],Pspeed,GTime);

if plots==1
  subplot(211)  
  plotonearth(1,1)
  hold on
  xS=sin(theta).*cos(phi);
  yS=sin(theta).*sin(phi);
  zS=cos(theta);
  a=plot3(xS,yS,zS,'v');
  [xT,yT,zT]=sph2cart(Tlon*pi/180,Tlat*pi/180,1);
  b=plot3(xT,yT,zT,'o');
  view(-Tlon,Tlat)
  set(b,'MarkerS',10,'MarkerF','k')
  for index=1:length(Slon)
    gces=grcircle([Slon(index) Slat(index)]/180*pi,[Tlon,Tlat]/180*pi,20);
    [x,y,z]=sph2cart(gces(:,1),gces(:,2),1);
    c(index)=plot3(x,y,z,'-');
  end  
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [del,azm,tim]=distat(epicenter,stations,speed,To)
del=grcdist(epicenter,stations);
azm=sazimuth(epicenter,stations);
tim=del/speed+To;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sortplot(az,delres)
[azs,azi]=sort(az);
plot(az(azi),delres(azi),'+-')


