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.

clf

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

% Random earthquake location
Tloc=randsph(1);
Ttime=rand(1)*3000;

% Random initial guess
Gloc=randsph(1);
Gtime=0;

% In longitude and latitude
[Tlon,Tlat]=cart2sph(Tloc(1),Tloc(2),Tloc(3));
[Glon,Glat]=cart2sph(Gloc(1),Gloc(2),Gloc(3));

% In degrees
Tlon=Tlon*180/pi;
Tlat=Tlat*180/pi;
Glon=Glon*180/pi;
Glat=Glat*180/pi;

% No negative longitudes
Glon=360*(Glon<0)+Glon;
Tlon=360*(Tlon<0)+Tlon;

% Make random distribution of stations
u=rand(N,1);
v=rand(N,1);
phi=2*pi*u;
theta=acos(2*v-1);
Slon=phi*180/pi;
Slat=90-theta*180/pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if plots==1
  ah(1)=subplot(221);
  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');
  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
  view(Tlon,Tlat)
  axis off
  ah(2)=subplot(222);
  plot(Slon,Slat,'+')
  hold on
  b=plot(Tlon,Tlat,'o','MarkerS',10,'MarkerF','r');
  axis equal
  axis([0 360 -90 90])
  box on
end

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

% Iterate until error is acceptable 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;
  if plots==1
    axes(ah(2))
    plot(Glon,Glat,'s','MarkerS',10,'MarkerF','g')
    xl(1)=xlabel('Longitude');
    yl(1)=ylabel('Latitude');
    ah(3)=subplot(212);
    sortplot(az,delres)
    tit=title(sprintf(...
	'%5.3fE, %5.3fN, %5.1fs\n%5.3fE, %5.3fN, %5.1fs, RMSE= %5.3f',...
	Tlon,Tlat,Ttime,...
	Glon,Glat,Gtime,rms(tpred-tobs)));
    set(tit,'HorizontalA','Center')
    xl(2)=xlabel('Azimuth');
    yl(2)=ylabel('Distance Difference');
    box on
    set([ah(2:end) xl yl tit],'FontS',15)
    if index==1; delete(ah(1)); end
    set(ah(2),'Position',[0.1300    0.5811    0.7750    0.3439])
    shrink(ah(3),1,1.2)
    shrink(ah(2),1,1.2)
    orient portrait
    figdisp
    keyboard
  end
    [md,mdi]=max(delres);
    % Shift Glon and Glat by md in the direction az(mdi)  
    [Glon,Glat]=grcazim([Glon Glat],md,360-az(mdi));
end


% Plot residuals
if plots==1
  axes(ah(2))
  f=plot(Glon,Glat,'s','MarkerS',10,'MarkerF','g');
  b=plot(Tlon,Tlat,'o','MarkerS',8,'MarkerF','r');
  axes(ah(3))
  [delta,az,tpred]=distat([Glon Glat],[Slon Slat],Pspeed,Gtime);
  tit=title(sprintf(...
      '%5.3fE, %5.3fN, %5.1fs\n%5.3fE, %5.3fN, %5.1fs, RMSE= %5.3f',...
      Tlon,Tlat,Ttime,...
      Glon,Glat,Gtime,rms(tpred-tobs)));
  set(tit,'HorizontalA','Center')
  box on
  set([ah(2:end) xl yl tit],'FontS',15)
end

disp(sprintf(...
    '%5.3fE, %5.3fN, %5.1fs\n%5.3fE, %5.3fN, %5.1fs, RMSE= %5.3f',...
    Tlon,Tlat,Ttime,...
    Glon,Glat,Gtime,rms(tpred-tobs)))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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),'+-','LineW',2)
grid on
xlim([0 360])

