Plot_results_toy_example_ac.m 4.33 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
% This is a matlab script to plot the results of the acoustic toy example of the DENISE code.

clear all
% close all
clc

iteration_no=100;

% general informations about the model
nx=150;		% grid points in x direction
ny=150;		% grid points in y direction
DH=10.0;		% spatial sampling interval
X=(1:1:nx)*DH;	% definition of x axis for plots
Y=(1:1:ny)*DH;	% definition of y axis for plots
FW=30;		% width of PML frame in gridpoints
sources=[350.0, 550.0, 750.0, 950.0, 1150.0];	% x coordinate of source positions

% general settings for plots information
fontsize=12;
x_line=[FW FW nx-FW nx-FW]*DH;
y_line=[1 ny-FW ny-FW 1]*DH;

% reading the true model
true_vp=read_binary_matrix(nx,ny,['../par/model/mod_toy_example_ac_true_vp_it0.bin']);
true_rho=read_binary_matrix(nx,ny,['../par/model/mod_toy_example_ac_true_rho_it0.bin']);

% reading initial model
start_vp=read_binary_matrix(nx,ny,['../par/model/toy_example/mod_toy_example_ac_vp_it0.bin']);
start_rho=read_binary_matrix(nx,ny,['../par/model/toy_example/mod_toy_example_ac_rho_it0.bin']);

% reading inverted vp model
inv_vp=read_binary_matrix(nx,ny,['../par/model/toy_example/mod_toy_example_ac_vp_it' int2str(iteration_no) '.bin']);

% finding minimum and maximum values for colorbar
max_vp=max([max(true_vp(:)) max(start_vp(:)) max(inv_vp(:))]);
min_vp=min([min(true_vp(:)) min(start_vp(:)) min(inv_vp(:))]);

max_rho=max([max(true_rho(:)) max(start_rho(:))]);
min_rho=min([min(true_rho(:)) min(start_rho(:))]);

% Plot true model
figure('units','normalized','outerposition',[0.8 0 0.2 1])
subplot(3,2,1)
imagesc(X,Y,true_vp)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true P-wave velocity model in m/s'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','vp in m/s'); 
caxis([min_vp max_vp])
axis equal tight

subplot(3,2,2)
imagesc(X,Y,true_rho)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['true density model in kg/m^3'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','\rho in kg/m^3'); 
caxis([min_rho max_rho])
axis equal tight

% Plot initial model
subplot(3,2,3)
imagesc(X,Y,start_vp)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['initial P-wave velocity model in m/s'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','vp in m/s'); 
caxis([min_vp max_vp])
axis equal tight

subplot(3,2,4)
imagesc(X,Y,start_rho)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['initial density model in kg/m^3'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','\rho in kg/m^3'); 
caxis([min_rho max_rho])
axis equal tight

% Plot obtained model
subplot(3,2,5)
imagesc(X,Y,inv_vp)
hold on
line(x_line,y_line,'LineStyle','--','Color','k')
plot(sources,10.0*ones(1,5),'*r','MarkerSize',12)
hold off
xlabel('x in m','FontSize',fontsize)
ylabel('y in m','FontSize',fontsize)
title(['inverted P-wave velocity model in m/s',10,'(iteration step ' int2str(iteration_no) ')'],'FontSize',fontsize)
colb=colorbar;
coll=get(colb,'ylabel');
set(coll,'String','vp in m/s'); 
caxis([min_vp max_vp])
axis equal tight

% Plot vertical velocity provfiles for S-wave velocity model
x1=500;
x2=750;

subplot(3,2,6)
plot(true_vp(:,round(x1/DH)),Y,'Color',[0.7 0.7 0.7],'LineWidth',6,'LineStyle','-');
hold on
plot(start_vp(:,round(x1/DH)),Y,'k','LineWidth',3,'LineStyle','--');
plot(inv_vp(:,round(x1/DH)),Y,'b','LineWidth',3);
plot(inv_vp(:,round(x2/DH)),Y,'r','LineWidth',3,'LineStyle','--');
line([min_vp-10 max_vp+10],[(ny-FW)*DH (ny-FW)*DH],'LineStyle','--','Color','k')
hold off
xlim([min_vp-10 max_vp+10])
set(gca,'ydir','reverse')
xlabel('velocity in m/s','fontsize',fontsize)
ylabel('depth in m','fontsize',fontsize)
title('Vertical velocity profiles','fontsize',fontsize)
legend('true','initial',['x=' num2str(x1) 'm'],['x=' num2str(x2) 'm']);
hold off
grid on