 
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=[ .2 .4 .2 .2];
%sobs=[ 0.5 1 0.5 0.5]; % standard deviations observations
xa=0; xe=35; xinc=0.5;			% upper and lower bounds
ya=0; ye=0; yinc=1;
za=0; ze=25; zinc=0.5;
Ta=10; Te=30; Tinc=0.2;
nx=(xe-xa)/xinc+1;
ny=(ye-ya)/yinc+1;
nz=(ze-za)/zinc+1;
nT=(Te-Ta)/Tinc+1;


% calculation of a priori proabailtiy in the data space
ta=28; te=31; tinc=0.5; 
i1=0; i2=0; i3=0; i4=0;
nt=(te-ta)/tinc+1;
disp(sprintf(' Proceed to calculate a priori pdf for %i samples ? ',nt))
pause

rhoD=zeros([nt nt nt nt]);
for t1=ta:tinc:te;
   i1=i1+1;
   i2=0;
   for t2=ta:tinc:te;
      i2=i2+1;
      i3=0;
      for t3=ta:tinc:te;
         i3=i3+1;
         i4=0;
         for t4=ta:tinc:te;
            i4=i4+1;
            rhoD(i1,i2,i3,i4)=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);
         end
      end
   end
end

rhoD=rhoD/max(max(max(max(rhoD))));

disp(' Continue to calculate marginals ? ')
pause

% Calculate the marginal probabilities and display in 3D space

% plane t1-t2

rhoD12=zeros([nt nt]);
rhoD13=zeros([nt nt]);
rhoD23=zeros([nt nt]);

for i=1:nt,
   for j=1:nt,
      rhoD12(i,j)=sum(sum(rhoD(i,j,:,:)));
      rhoD13(i,j)=sum(sum(rhoD(i,:,j,:)));
      rhoD23(i,j)=sum(sum(rhoD(:,i,j,:)));
   end
end

% normalize marginals for plotting
rhoD12=rhoD12/max(max(rhoD12));
rhoD13=rhoD13/max(max(rhoD13));
rhoD23=rhoD23/max(max(rhoD23));

% show information on data
t=ta:tinc:te;
t2=t'*ones([1 nt]);

surf(t2,t2',t2*0+ta,rhoD12),hold on
surf(t2,t2*0+te,t2',rhoD13),
surf(t2*0+te,t2,t2',rhoD23),
xlabel(' t_1 (s) ')
ylabel(' t_2 (s) ')
zlabel(' t_3 (s) ')
shading faceted
hold off



% Now let s solve the inverse problem

disp(sprintf(' Proceed to calculate a priori 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 psoteriori probability

% plane t1-t2

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;
imagesc(x,z,sigDxz');

% plot 
               

