Commit 69e447f9 authored by tilman.metz's avatar tilman.metz

some changes to rheinstetten sytnthetic test

parent 4afe7b37
......@@ -16,6 +16,8 @@ nx=220; ny=100; nz=220; %ny:vertical
outx=1; outy=1; outz=1;
dh=0.15;
FW=10;
Free_Surf=1;
nx=nx/outx;ny=ny/outy;nz=nz/outz;
X=dh:dh*outx:nx*outx*dh;
......@@ -53,11 +55,18 @@ REC=textscan(fp, '%f %f %f ', 'delimiter','\t', 'commentstyle','#','MultipleDeli
fclose(fp);
end
% Lines to mark PML position
fontsize=14;
x_line=[FW FW nx-FW nx-FW FW]*dh;
y_line=[FW ny-FW ny-FW FW FW]*dh;
z_line=[FW nz-FW nz-FW FW FW]*dh;
if(Free_Surf==1)
y_line=[0 ny-FW ny-FW 0 0]*dh;
else
y_line=[FW ny-FW ny-FW FW FW]*dh;
end
%m --> model indices
......@@ -187,9 +196,9 @@ axis tight
hold on
%Plot sources if in model slice
for ii=1:length(SRC{1})
if SRC{2}(ii)==sly;
% if SRC{2}(ii)==sly;
plot(SRC{1}(ii),SRC{3}(ii),'*blue','MarkerSize',15,'LineWidth',2.0);
end
% end
end
%Plot receiver if in model slice
......
%3D plot of model
clear all;
nx=220; ny=100; nz=220; %ny:vertical
outx=1; outy=1; outz=1;
dh=0.15;
nx=nx/outx;ny=ny/outy;nz=nz/outz;
fignum=12;
%file_inp1='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy_real.vp';
%file_inp2='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/model/toy_real.vp';
file_inp1='../../par/model/rheinstetten.vs_it0';
file_inp2='../../par/model/rheinstetten.vs_it60';
phi1=0; %rotation angles to x-z plane of first and second plane
phi2=90;
rotaxes=[1,0,0]; %direction rotation axes [0,1,0] rotation of plane around vertical axis
% [1,0,0] rotation of plane around x-axis
%rotpoint=[120.8 140.8 80.8] ; %turning point [x y z] in meter
rotpoint=[(nx/4)*dh (15)*dh (3*nz/4)*dh] ;
%rotpoint=[56 60 60] ;
rotpoint2=[110 130 70] ;
%viewpoint=[122,26];
viewpoint=[-20,12];
caxis_value_1=0;%vs
caxis_value_2=300;%vs
%caxis_value_1=0;%vp
%caxis_value_2=2000;%vp
xcut1=1; xcut2=nx;
ycut1=1; ycut2=ny;
zcut1=1; zcut2=nz;
fid_rot=fopen(file_inp1,'r','ieee-le');
rot1=zeros(ny/outy,nx/outx,nz/outz);
rot1=fread(fid_rot,(nx*ny*nz),'float');
% fid_div=fopen(file_inp2,'r','ieee-le');
% div1=zeros(nz/outz,nx/outx,ny/outy);
% div1=fread(fid_div,(nx*ny*nz),'float');
%rot1=-rot1./max(max(rot1));
%rot1=rot1-div1;
%max(max(rot1))
%rot1=1./((rot1+3e-10));
%rot1=rot1.*div1*50000;
%rot1=-rot1./max(max(abs(rot1)));
%rot1=log10(rot1);
rot=reshape(rot1,ny/outy,nx/outx,nz/outz);
rot=rot(ycut1:ycut2,xcut1:xcut2,zcut1:zcut2);
nx=xcut2-xcut1+1;
ny=ycut2-ycut1+1;
nz=zcut2-zcut1+1;
%size(rot)
xp1=xcut1*dh; xp2=xcut2*dh; yp1=ycut1*dh; yp2=ycut2*dh; zp1=zcut1*dh; zp2=zcut2*dh;
x=xp1:dh*outx:xp2*outx;
y=yp1:dh*outy:yp2*outy;
z=zp1:dh*outz:zp2*outz;
%x=xcut1:dh*outx:xcut2*outx;
%y=ycut1:dh*outy:ycut2*outy;
%z=zcut1:dh*outz:zcut2*outz;
figure(12)
[Z,X,Y]=meshgrid(z,x,y);
xmin = min(X(:)); ymin = min(Y(:)); zmin = min(Z(:));
xmax = max(X(:)); ymax = max(Y(:)); zmax = max(Z(:));
%zplane=zeros(nx,nz);
%zplane(:,:)=rotpoint2(2);
%hslicez = surf(z,x,zplane);
%rotate(hslicez,rotaxes,phi1,rotpoint);
%xd1 = get(hslicez,'XData');
%yd1 = get(hslicez,'YData');
%zd1 = get(hslicez,'ZData');
%zplane=zeros(ny,nx);
%zplane(:,:)=rotpoint(3);
%hslicez = surf(x,y,zplane);
%rotate(hslicez,rotaxes,phi1,rotpoint);
%xd = get(hslicez,'XData');
%yd = get(hslicez,'YData');
%zd = get(hslicez,'ZData');
xd2 = rotpoint(1);
yd2 = rotpoint(2);
zd2 = rotpoint(3);
%zplane=zeros(nz,ny);
%zplane(:,:)=rotpoint(3);
%hslicez = surf(y,z,zplane);
%rotate(hslicez,rotaxes,phi2,rotpoint);
%xd2 = get(hslicez,'XData');
%yd2 = get(hslicez,'YData');
%zd2 = get(hslicez,'ZData');
%zplane=zeros(ny,nx);
%zplane(:,:)=rotpoint(3);
%hslicez = surf(x,y,zplane);
%rotate(hslicez,rotaxes,phi2,rotpoint);
%xd2 = get(hslicez,'XData');
%yd2 = get(hslicez,'YData');
%zd2 = get(hslicez,'ZData');
zplane=zeros(ny,nx);
zplane(:,:)=rotpoint(3);
hslicez = surf(x,y,zplane);
rotate(hslicez,[0 1 0],phi2,rotpoint);
xd3 = get(hslicez,'XData');
yd3 = get(hslicez,'YData');
zd3 = get(hslicez,'ZData');
%zplane=zeros(ny,nx);
%zplane(:,:)=rotpoint2(3);
%hslicez = surf(x,y,zplane);
%rotate(hslicez,[1 0 0],90,rotpoint2);
%xd4 = get(hslicez,'XData');
%yd4 = get(hslicez,'YData');
%zd4 = get(hslicez,'ZData');
rot=permute(rot,[2,3,1]);
%rott=size(rot)
%size(rot)
%size(Z)
figure(fignum)
%h3 = slice(Z,X,Y,rot,zd1,xd1,yd1);set(h3,'FaceColor','interp',...
% 'EdgeColor','none',...
% 'DiffuseStrength',.8)
% alpha(.5)
%hold on
%h = slice(Z,X,Y,rot,zd,xd,yd);set(h,'FaceColor','interp',...
% 'EdgeColor','none',...
% 'DiffuseStrength',.8)
%hold on
h2 = slice(Z,X,Y,rot,zd2,xd2,yd2);set(h2,'FaceColor','interp',...
'EdgeColor','none',...
'DiffuseStrength',.8)
hold on
% h2 = slice(Z,X,Y,rot,zd3,xd3,yd3);set(h2,'FaceColor','interp',...
% 'EdgeColor','none',...
% 'DiffuseStrength',.8)
hold on
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut1*dh],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut2*dh],[ycut2*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut2*dh xcut2*dh],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut2*dh],[ycut1*dh ycut1*dh],'-black','LineWidth',0.5);
hold on
plot3([zcut1*dh zcut2*dh],[xcut1*dh xcut1*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
plot3([zcut2*dh zcut2*dh],[xcut1*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut2*dh],[xcut2*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut1*dh],[xcut1*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'-black','LineWidth',0.5);
hold on
plot3([zcut1*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut1*dh],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[ycut2*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([zcut1*dh zcut1*dh],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
plot3([zcut2*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut2*dh],'-black','LineWidth',0.5);
hold on
plot3([zcut1*dh zcut2*dh],[rotpoint(1) rotpoint(1)],[rotpoint(2) rotpoint(2)],'k--','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[xcut1*dh xcut2*dh],[rotpoint(2) rotpoint(2)],'k--','LineWidth',0.5);
plot3([rotpoint(3) rotpoint(3)],[rotpoint(1) rotpoint(1)],[ycut1*dh ycut2*dh],'k--','LineWidth',0.5);
% hold on
% hold on
% for(i=1:13)
% for(kk=1:13)
% plot3([20*dh+(kk-1)*10*dh 20*dh+(kk-1)*10*dh],[20*dh+(i-1)*10*dh 20*dh+(i-1)*10*dh],[24 24],'+blue','MarkerSize',3,'LineWidth',2.0);
% end
% end
% hold on
hold off
xlabel('z in m');
ylabel('x in m');
zlabel('y in m');
set(gca, 'xDir','normal')
set(gca, 'yDir','normal')
set(gca, 'zDir','reverse')
set(get(gca,'Ylabel'),'FontSize',14);
set(get(gca,'Ylabel'),'FontWeight','normal');
set(get(gca,'Xlabel'),'FontSize',14);
set(get(gca,'Xlabel'),'FontWeight','normal');
set(get(gca,'Zlabel'),'FontSize',14);
set(get(gca,'Zlabel'),'FontWeight','normal');
set(gca,'FontSize',14);
set(gca,'FontWeight','normal');
set(gca,'Linewidth',1.0);
colormap('jet')
cb = colorbar('vert');
xlabel(cb, 'v_p in m/s');
%set(cb,'YTick',[5800:200:6600]); % v_p
set(gca,'FontSize',14);
%load('MyColormapsgrad','mycmap');
% set(figure(fignum),'Colormap',mycmap)
view(viewpoint);
xlim([xcut1*dh xcut2*dh]);
ylim([zcut1*dh zcut2*dh]);
zlim([ycut1*dh ycut2*dh]);
caxis([caxis_value_1 caxis_value_2])
daspect([1,1,1]);
% axis tight
box on
% exportfig(fignum, 'toy_model_real_vp.eps','bounds','tight', 'color','rgb', ...
% 'preview','none', 'resolution',200, 'lockaxes',1);
......@@ -10,7 +10,7 @@
8.25 0.15 17.25 0.000 31.5 1.0
8.25 0.15 21.75 0.000 31.5 1.0
8.25 0.15 26.25 0.000 31.5 1.0
3.75 0.15 3.75 0.000 31.5 1.0
12.75 0.15 3.75 0.000 31.5 1.0
12.75 0.15 8.25 0.000 31.5 1.0
12.75 0.15 12.75 0.000 31.5 1.0
12.75 0.15 17.25 0.000 31.5 1.0
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment