
clear 
close all

% Hypo2002_exacttheory
% a programme to use probabilistic inverse theory 
% to determine hypocenter locations (highly simplified)
% here simplified to a 2D problem, but extendable to 3D

% model and data
ns=4;								% number of stations (and observations)
sx=[0 10 15 20]; 	    			% station coordinates x (km)
sy=[0 0 0 0];  	                    % station coordinates y (km)
sz=[0 0 0 0];	                    % station coordinates z (km)
v=5.;								% homogeneous velocity (km/s)
xs=20;                              % source location x
ys=0;                               % source location y
zs=8;                               % source location z
t0=20;                              % origin time
ns=length(sx);						% number of stations (and observations)


% Calculate exact data
for i=1:ns,
    tobs(i)=t0+1/v*sqrt( (sx(i)-xs)^2 + (sy(i)-ys)^2 + (sz(i)-zs)^2);
end
t=tobs*0;         % vector used later

% Define uncertainties in data using variances (length of vector is number
% of data points
%sobs=[ 1 2 1 1];
%sobs=[ .5 1 .5 .5];
sobs=[ .5 .3 .2 .4];
if length(sobs) ~= ns, disp(' Please adjust vector length (sobs) '), pause, end
%sobs=[ 0.5 1 0.5 0.5]; 

% Define model space and increments for model exploration
xa=0; xe=35; xinc=.5;			% upper and lower bounds in x0
ya=0; ye=0; yinc=.5;            % y0
za=0; ze=25; zinc=1.;           % z0
Ta=10; Te=30; Tinc=.25;         % t0
nx=(xe-xa)/xinc+1;              % number of sample in dimension x0
ny=(ye-ya)/yinc+1;              % y0
nz=(ze-za)/zinc+1;              % z0
nT=(Te-Ta)/Tinc+1;              % t0

% Now let s solve the inverse problem
% Some parameters
ntot=nx*ny*nz*nT;               % total number of forward modellings to do
nm=round(ntot/100);                    % info after ntot/100 samples 
nn=0;

disp(sprintf(' Proceed to calculate a posteriori pdf for %i x %i x %i x %i = %i samples ? ',nx,ny,nz,nT,ntot))
disp(' Press any key to continue ')
pause
disp(' Starting exploration of model space ... ')



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

sigD=zeros([nx ny nz nT]);      % initilization of pdf to 0

% Loop through model space

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
            
            for is=1:ns,    % loop over number of observables
            t(is)=T+1/v*sqrt((x-sx(is))^2+(y-sy(is))^2+(z-sz(is))^2);
            end
        
            % calculate the associated probability
            
                term=1;
            for is=1:ns,    % loop over number of observables
                term=term*exp(-1/(2*sobs(is)^2)*(t(is)-tobs(is))^2);
            end
            sigM(ix,iy,iz,iT)=term;         % 
          
                                      nn=nn+1;
                                      if rem(nn,nm)==0, 
                                         disp(sprintf(' Calculated %g %%',round(nn/ntot*100)))
                                      end
                                      
                          end
    %                 end
                  end
            end
               
               
               
% Calculate the marginal a posteriori probability

% plane x-z

sigMxz=zeros([ix iz]);

for i=1:ix,
   for j=1:iz,
      sigMxz(i,j)=sum(sum(sigM(i,:,j,:)));
      end
end

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

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

%  marginal a posteriori for T

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

figure
imagesc(interp(x,4),interp(z,4),interp2(sigMxz',4)),axis equal, axis([xa xe za ze])
xlabel(' x(km) ')
ylabel(' z(km) ')
title(' A posteriori probability: hypocenter ')
figure
plot(Ta:Tinc:Te,sigMT), 
xlabel(' Origin time (s) ')
ylabel(' Probability')
title(' A posteriori probability: origin time ')


