
clear 
clf

% Hypo2002
% a programme to use probabilistic inverse theory 
% to determine hypocenter locations (highly simplified)

% model and data
ns=4;								% number of stations (and observations)
sx=[5 10 15 20]; 				% station coordinates (km)
sy=[0 0 0 0];  
sz=[0 0 0 0];
v=5.;								% homogeneous velocity (km/s)
tobs=[ 30.3 29.4 ...			% observations (travel times in s)
      28.6 28.3];
%sobs=[ 1 2 1 1];
%sobs=[ .5 1 .5 .5];
sobs=[ .2 .4 .2 .2];
%sobs=[ 0.5 1 0.5 0.5]; % standard deviations observations
xa=0; xe=35; xinc=.5;			% upper and lower bounds
ya=0; ye=0; yinc=.5;
za=0; ze=25; zinc=1.;
Ta=10; Te=30; Tinc=.25;
nx=(xe-xa)/xinc+1;
ny=(ye-ya)/yinc+1;
nz=(ze-za)/zinc+1;
nT=(Te-Ta)/Tinc+1;


% Now let s solve the inverse problem

disp(sprintf(' Proceed to calculate a posteriori pdf for %i x %i x %i x %i samples ? ',nx,ny,nz,nT))
pause

ntot=nx*ny*nz*nT;
nm=ntot/1000;
nn=0;

% We now loop through the model space (x,y,z, and T)

sigD=zeros([nx ny nz nT]);

ix=0;
for x=xa:xinc:xe,
   ix=ix+1;
   iy=0;
 %  for y=ya:yinc:ye,
 %     iy=iy+1;
 			iz=0;
 			iy=1;
			 y=0;
      for z=za:zinc:ze,
         iz=iz+1;
         iT=0;
        for T=Ta:Tinc:Te,
            iT=iT+1;
            % calculate forward problem (z=0)
            t1=T+1/v*sqrt((x-sx(1))^2+(y-sy(1))^2+z^2);
            t2=T+1/v*sqrt((x-sx(2))^2+(y-sy(2))^2+z^2);
            t3=T+1/v*sqrt((x-sx(3))^2+(y-sy(3))^2+z^2);
            t4=T+1/v*sqrt((x-sx(4))^2+(y-sy(4))^2+z^2);
            
            %disp(sprintf(' T, x,z, t: %g %g %g %g %g %g %g ',T,x,z,t1,t2,t3,t4))
            
            % calculate the associated probability
            sigD(ix,iy,iz,iT)=exp(-1/(2*sobs(1)^2)*(t1-tobs(1))^2)*...
               					exp(-1/(2*sobs(2)^2)*(t2-tobs(2))^2)*...
                              exp(-1/(2*sobs(3)^2)*(t3-tobs(3))^2)*...
                              exp(-1/(2*sobs(4)^2)*(t4-tobs(4))^2);
          
                                      nn=nn+1;
                                      if rem(nn,nm)==0, 
                                         disp(sprintf(' Calculated %g \%%',nn/ntot*100))
                                      end
                                      
                                         
                           
                          end
    %                 end
                  end
            end
               
               
               
% Calculate the marginal a posteriori probability

% plane x-z

sigDxz=zeros([ix iz]);

for i=1:ix,
   for j=1:iz,
      sigDxz(i,j)=sum(sum(sigD(i,:,j,:)));
      end
end

% normalize marginals for plotting
sigDxz=sigDxz/max(max(sigDxz));

x=xa:xinc:xe;
z=za:zinc:ze;

% interpolate and plot

imagesc(interp(x,4),interp(z,4),interp2(sigDxz',4))

%  marginal a posteriori for T

for i=1:nT,sigDT(i)=sum(sum(sum(sigD(:,:,:,i)))); end
s=sum(sigDT);
sigDT=sigDT/s;

subplot(121)
imagesc(interp(x,4),interp(z,4),interp2(sigDxz',4))
xlabel(' x(km) ')
ylabel(' z(km) ')
title(' A posteriori probability: hypocenter ')
subplot(122)
plot(Ta:Tinc:Te,sigDT)
xlabel(' Origin time (s) ')
ylabel(' Probability')
title(' A posteriori probability: origin time ')


