
% Matlab program to illustrate inverse problems with the gravitational 
% potential problem 

% Problem: we discretize the subsurface with blocks and calculate the 
% gravitational potential at the surface for a buried box with gold (the data) 
% In order to find the location of the box, we randomly peturb the model 
% and calculate the corresponding potential we would observe (synthetics).
% Then we comparee observations with our model guess (misfit). We calculate 
% randomly (trial and error method) many models and collect those models with
% a misfit better than "crit" in percent.

% Exercise:
% - Save the original Matlab script in another file 
% - Change the misfit criterion (crit) to smaller or larger values to
% investigate the chance to find "good" solutions
% - Assume that the box consists of two pieces (2x2, 2x2) that are
% distributed in  random way and search for good solutions
% - Further possibilities:
% - change the number of data points (receivers in rec)
% - assume one box with unknown size (1x1 to 4x4) and search for it 
 
rec=[ 1 3 5 7 9 ];  			% Receiver locations
nr=length(rec);				% nr number of receivers
nm=10;							% nm number of acceptable models to collect 
n=10;								% dimension of model space (n x n)
maxiter=1000000;				% max number of iterations
crit=0.1;							% Criterionj for good models (misfit smaller than crit (in %))

model=zeros(n);				% initialization of model matrix

sand=2.2;						% density of sand
gold=19.3;						% density of gold

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE FORWARD PROBLEM

% Calculation of true solution

% initialise true model

model=model+sand;									% bbackground is sand
model(3:5,4:5)=model(3:5,4:5)*0+gold;		%	box with gold

% calculate observed data

rdata=zeros([1 nr]); 		% initalite real data vector
for i=1:nr,
    for j=1:n,
		for k=1:n,
      rdata(i)= rdata(i)+model(j,k)/sqrt( (j-rec(i))^2 + k^2 );   % sum of potential from
      end																			% each block
    end
 end
 
% plot real data
 
subplot(211),plot(1:nr,rdata,'-',1:nr,rdata,'k+'), title(' Observed data ')
subplot(212),pcolor(model'), axis ij, title(' True model')
rmodel=model;

disp(' Press any key to continue with the inversion ... ')
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE INVERSE PROBLEM

% Loop through model space randomly
% collect models if fit is good

ig=0;										% amount of "good, acceptable" models
good_model=zeros([nm,n,n]);		% the actual "good, acceptable" models
good_data=zeros([nm,nr]);			% the data of the "good, acceptable" models

while ig< nm,							% initialize new randomly chosen models
model=model*0+sand;					% while the number of good models ig is smaller  
ix=floor(rand([1 6])*10)+1;		% than the number of desired (nm) good models
iy=floor(rand([1 6])*10)+1;

for j=1:6,
model(ix(j),iy(j))=gold;			% the model location of gold
end

% calculate synthetic data for the random model

data=zeros([1 nr]); 

for i=1:nr,
   for j=1:n,
	  for k=1:n,
	     data(i)= data(i)+model(j,k)/sqrt( (j-rec(i))^2 + k^2 );
     end
   end
end

% calculate the misfit with the real data

%err=sum(abs(rdata-data))/sum(rdata)*100;    % error (L1-norm)
diff=sum((rdata-data).*(rdata-data));			% root mean square error
f=sum(rdata.*rdata) ;
err=1/nr^2*sqrt(diff/f)*100;


% Gather models that are good
if err < crit,
ig=ig+1;
good_model(ig,:,:)=model;
good_data(ig,:,:)=data;
good_err(ig)=err;
disp(sprintf(' Found good model No. %g ',ig));
end

% plot the calculated model whether good or bad

subplot(211),plot(1:nr,rdata,'-',1:nr,rdata,'k+',1:nr,data,'r-',1:nr,data,'r+'), title(sprintf(' Predicted (red) and Observed (blue) data: Diff: %g % ',err)), axis([1 5 40 80])
subplot(212),pcolor(model'), axis ij, title(' Current model')
drawnow

%pause(.05)


end


% at the very end plot the nm "acceptable, good" models collected,
% the data and compare to the real model

ii=0;

for i=1:12,
if i==1, subplot(3,4,i)
 plot(rdata,'b-'),hold on, plot(rdata,'b+')
 for j=1:10, plot(good_data(j,:),'r-'),end
  hold off
elseif i==2, subplot(3,4,i), 
  pcolor(rmodel), axis ij, title(' True model ')
else
  subplot(3,4,i)
  s=zeros(n);for j=1:n, for k=1:n, s(j,k)=good_model(i-2,j,k); end, end
  pcolor(s),axis ij, 
  title(sprintf(' Model %g : Diff: %g %',i-2,good_err(i-2)));
end

end



% End of the story, did you find the gold?