Commit 576d1dc8 authored by tilman.metz's avatar tilman.metz

Merge branch 'develop'

Conflicts:
	src/ifos3d.c
	src/steplength.c
parents 01c47bcf a65059e3
clear all
%close all
% Definition von Parametern
% Homogener Vollraum
npts=10000;
delta=0.0001;
t=(0:delta:(npts-1)*delta);
% Parameter: Homogener Vollraum
%T_d=0.004;
%T_d=0.0033;
T_d=1/300.0;
%T_d=0.002;
F0=1;
npts_start=1;
npts_end=T_d/delta;
wavelet=zeros(1,npts);
% Definition des Quellwavelets
wavelet(1:npts_end)=F0*(sin((pi*t(1:npts_end))/T_d)).^3;
%figure
%plot(t,wavelet)
%xlabel('Zeit in s')
WAVELET=delta*fft(wavelet);
maximum=max(max(abs(WAVELET)))
WAVELET=WAVELET./maximum;
frequenz=[0:fix(npts/2),-ceil(npts/2)+1:-1]'/(npts*delta); %Definition des kompletten Frequenzvektors
f=[0:fix(npts/2)]'/(npts*delta);
%{
figure
subplot(211)
plot(f(1:fix(npts/2)+1),abs(WAVELET(1:fix(npts/2)+1)))
xlab=xlabel('Frequenz in Hz')
set(gca,'FontSize',14)
tit=title(['Laenge des Quellsignals: ' num2str(T_d) ' s']);
set(xlab,'FontSize',14)
subplot(212)
plot(f(1:fix(npts/2)+1),abs(WAVELET(1:fix(npts/2)+1)))
xlab=xlabel('Frequenz in Hz')
set(gca,'FontSize',14)
ylim([0 0.01])
%set(gca,'XTick',[0 30 50 70 90 110 150 200 250])
%set(gca,'XTickLabel',{'0';'30';'50';'70';'90';'110';'150';'200';'250'})
set(xlab,'FontSize',14)
%}
figure
plot(t,wavelet,'-k')
xlab=xlabel('Time in s');
ylab=ylabel('s(t) (dimensionless)');
set(gca,'FontSize',18)
set(xlab,'FontSize',20)
set(ylab,'FontSize',20)
xlim([0 1.5*T_d])
figure
plot(f(1:fix(npts/2)+1),abs(WAVELET(1:fix(npts/2)+1)),'-k')
xlab=xlabel('Frequency in Hz');
ylab=ylabel('|S(\omega)| in 1/Hz');
set(gca,'FontSize',18)
set(xlab,'FontSize',20)
set(ylab,'FontSize',20)
%tit=title(['Laenge des Quellsignals: ' num2str(T_d) ' s']);
%set(xlab,'FontSize',14)
xlim([0 600])
%==========================================================================
% BINREAD by André Kurzmann
% read 2D arrays from binary files
%==========================================================================
function A = binread(fname,ny,nx,fmt)
if nargin < 4
fmt = 'float32';
end
f = fopen(fname,'r');
A = fread(f,inf, fmt);
fclose(f);
if nargin == 2
nx = length(A)/ny;
A = reshape(A,ny,nx);
end
if nargin == 3
A = reshape(A,ny,nx);
end
function structname=create_su_struct();
% This function creates a matlab struct containing the
% su header words as fields and an additional field 'trace'
% for the actual seismogram trace.
%
% All fields are set to zero apart from
% .trid=1 (seismic trace)
% .scalel=1
% .scalco=1
% .counit=1
% .gain=1
% .igc=1
% .igi=1
% .corr=1 (uncorrelated traces)
structname = struct('tracl',0,...
'tracr',0,...
'fldr',0,...
'tracf',0,...
'ep',0,...
'cdp',0,...
'cdpt',0,...
'trid',1,...
'nvs',0,...
'nhs',0,...
'duse',0,...
'offset',0,...
'gelev',0,...
'selev',0,...
'sdepth',0,...
'gdel',0,...
'sdel',0,...
'swdep',0,...
'gwdep',0,...
'scalel',1,...
'scalco',1,...
'sx',0,...
'sy',0,...
'gx',0,...
'gy',0,...
'counit',1,...
'wevel',0,...
'swevel',0,...
'sut',0,...
'gut',0,...
'sstat',0,...
'gstat',0,...
'tstat',0,...
'laga',0,...
'lagb',0,...
'delrt',0,...
'muts',0,...
'mute',0,...
'ns',0,...
'dt',0,...
'gain',1,...
'igc',1,...
'igi',1,...
'corr',1,...
'sfs',0,...
'sfe',0,...
'slen',0,...
'styp',0,...
'stas',0,...
'stae',0,...
'tatyp',0,...
'afilf',0,...
'afils',0,...
'nofilf',0,...
'nofils',0,...
'lcf',0,...
'hcf',0,...
'lcs',0,...
'hcs',0,...
'year',0,...
'day',0,...
'hour',0,...
'minute',0,...
'sec',0,...
'timbas',0,...
'trwf',0,...
'grnors',0,...
'grnofr',0,...
'grnlof',0,...
'gaps',0,...
'otrav',0,...
'd1',0,...
'f1',0,...
'd2',0,...
'f2',0,...
'ungpow',0,...
'unscale',0,...
'ntr',0,...
'mark',0,...
'shortpad',0,...
'unass',0,...
'trace',0);
%plots misfit normalised to initial misfit value
clear all;
load /data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/misfit_toy.txt ;
misfit1=misfit_toy./misfit_toy(1)';
x=1:60;
figure(44)
plot(x,misfit1,'b')
xlabel('iteration');
ylabel('normalised misfit');
set(get(gca,'Ylabel'),'FontSize',14);
set(get(gca,'Ylabel'),'FontWeight','normal');
set(get(gca,'Xlabel'),'FontSize',14);
set(get(gca,'Xlabel'),'FontWeight','normal');
%calculates the error between the true and inverted model
clear all;
nx=300; ny=350; nz=200;
cuty1=40;
cuty2=290;
outx=1; outy=1; outz=1;
dh=0.8;
nx=nx/outx;ny=ny/outy;nz=nz/outz;
file_inp1='/workag13/FWI/random_KTB/Model/random7.vs';
file_inp2='/workag13/FWI/random_KTB/2D_10sources_vert/model/random72D_vert.vs_it73';
fid_rot=fopen(file_inp1,'r','ieee-le');
rot1=zeros(1,(nx*ny*nz));
rot1=fread(fid_rot,(nx*ny*nz),'float');
%rot=reshape(rot1,(nx*ny*nz));
size(rot1)
fid_div=fopen(file_inp2,'r','ieee-le');
div=zeros(nz/outz,nx/outx,ny/outy);
div1=fread(fid_div,(nx*ny*nz),'float');
div=reshape(div1,nz/outz,nx/outx,ny/outy);
error=0.0;
for i=1:(nx*ny*nz)
error=error+abs((rot1(i)-div1(i))/rot1(i));
end
error=(error/(nx*ny*nz))
% plots vertical velocity profile of real, inverted and starting model
% S. Butzer 2013
clear all;
%close all;
nx=320; ny=160; nz=320;
cuty1=1;
cuty2=320;
outx=1; outy=1; outz=1;
dh=0.8;
nx=nx/outx;ny=ny/outy;nz=nz/outz;
fignum=44;
file_inp1='/data14/sdunkl/surface/results12/model/surface11_real.vs';
file_inp2='/data14/sdunkl/surface/results12/model/surface11_start.vs';
file_inp3='/data14/sdunkl/surface/results12/model/surface12.vs_it170';
px=160; %location of profile
pz=200;
%--------------------------------------------------------------------------
nxm=nx*dh+dh; nym=ny*dh+dh; nzm=nz*dh+dh;
xp1=dh; xp2=nx*dh; yp1=dh; yp2=ny*dh; zp1=dh; zp2=nz*dh;
x=xp1:dh*outx:xp2*outx;
y=yp1:dh*outy:yp2*outy;
z=zp1:dh*outz:zp2*outz;
fid_rot=fopen(file_inp1,'r','ieee-le');
rot1=zeros(nz/outz,nx/outx,ny/outy);
rot1=fread(fid_rot,(nx*ny*nz),'float');
MEAN=mean(rot1)
%rot1=rot1-div1;
%rot1=1./((rot1+5e-7));
%min(min(rot))
%max(max(rot1))
%rot1=rot1./max(max(rot1));
%
%norm222=norm(rot1)
%rot1=rot1.*div1;
%rot1=rot1./max(max(abs(rot1)));
%rot1=log10(rot1);
MAX=max(max(rot1))
MIN=min(min(rot1))
rot=reshape(rot1,nz/outz,nx/outx,ny/outy);
depth=rot(pz,px,:);
depth1=depth(:);
A=size(depth1)
%--------------------------------------------------------------------------
fid_rot=fopen(file_inp2,'r','ieee-le');
rot1=zeros(nz/outz,nx/outx,ny/outy);
rot1=fread(fid_rot,(nx*ny*nz),'float');
MEAN=mean(rot1)
%rot1=rot1-div1;
%rot1=1./((rot1+5e-7));
%min(min(rot))
%max(max(rot1))
%rot1=rot1./max(max(rot1));
%
%norm222=norm(rot1)
%rot1=rot1.*div1;
%rot1=rot1./max(max(abs(rot1)));
%rot1=log10(rot1);
MAX=max(max(rot1))
MIN=min(min(rot1))
rot=reshape(rot1,nz/outz,nx/outx,ny/outy);
depth=rot(pz,px,:);
depth2=depth(:);
A=size(depth2)
%--------------------------------------------------------------------------
fid_rot=fopen(file_inp3,'r','ieee-le');
rot1=zeros(nz/outz,nx/outx,ny/outy);
rot1=fread(fid_rot,(nx*ny*nz),'float');
MEAN=mean(rot1)
%rot1=rot1-div1;
%rot1=1./((rot1+5e-7));
%min(min(rot))
%max(max(rot1))
%rot1=rot1./max(max(rot1));
%
%norm222=norm(rot1)
%rot1=rot1.*div1;
%rot1=rot1./max(max(abs(rot1)));
%rot1=log10(rot1);
MAX=max(max(rot1))
MIN=min(min(rot1))
rot=reshape(rot1,nz/outz,nx/outx,ny/outy);
depth=rot(pz,px,:);
depth3=depth(:);
A=size(depth3)
figure(fignum)
plot(depth1,y,'k-');hold on;
plot(depth2,y,'b-');hold on;
plot(depth3,y,'r-');
xlabel('velocity in m/s');
ylabel('depth in m');
set(gca, 'xDir','normal')
set(gca, 'yDir','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);
xlim([350 1000]);
ylim([1*dh 140*dh]);
close all
clear all
% Q-approximation using an improved tau-method (Blanch et al., 1995)
% otimization routine: Marquardt-Levenberg
global L w Qf1 Qf2
%-------------------------INPUT-PARAMETERS----------------------------
Q0=50.0; % constant Q to be approximated
fp1=5.0;fp2=200;df=0.1; % within frequency range fp1,..., fp2
L=10; % number of relaxation mechanisms
fl_st=[0.01 0.1 1 5 10 70 200 500 1000 10000]; % L starting values for the relaxation frequencies
t=2/Q0; % starting value for tau
%-------------------------END: INPUT-PARAMETER-----------------------
f=fp1:df:fp2;
w=2*pi*f;
Qf1=Q0+f*0; % here constant Q-approximation, arbitrary
% frequency dependency of Q possible, define
% Q as function of frequency here
Qf1=Q0+sin(2*pi*f/fp2)*Q0;
% defining option for otimization (see 'help foptions')
options(1)=1;
options(2)=0.1;
options(3)=0.1;
options(5)=0; % =0 : Levenberg-Marquardt Method, =1: Gauss-Newton Method
options(14)=500;
x=[fl_st t];
% if optimization toolbox is installed use 'leastq'
[x,options]=leastsq('qflt',x,options);
% output of results:
RELAXATIONSFEQUENZ=x(1:L),
TAU=x(L+1),
sigma=sqrt(sum((Qf2-Qf1).*(Qf2-Qf1))/size(Qf2,2))*100/Q0;
GEW_REL_STANDARDABW=sigma,
% plot Q as function of frequency:
plot(f,Qf1,f,Qf2,'r--');
xlabel('frequency [hz]');
ylabel('quality factor (Q)');
%axis([fp1 fp2 10 100]);
% This function computes the difference between a
% frequency indepent quality factor and Q as function of
% relaxation frequencies and tau (written for optimization with leastsq).
function delta=qstd(x)
global L w Qf1 Qf2
fl=x(1:L);
t=x(L+1);
% computing telaxation times and relaxation frequencies
ts=1./(2*pi*fl);
% Q for a generalized standard linear solid:
sumzQ=0;sumnQ=0;
for l=1:L,
d=1+w.*w*ts(l)*ts(l);
sumnQ=(w*ts(l)*t./d)+sumnQ;
sumzQ=((w.*w*ts(l)*ts(l)*t)./d)+sumzQ;
end
Qf2=(1+sumzQ)./sumnQ;
delta=(Qf2-Qf1);
function q=qgsls(te,ts,L,w)
% Q fuer den L-fachen standard linear solid:
sumzQ=0;sumnQ=0;
for l=1:L,
d=1.0+w.*w*ts(l)*ts(l);
sumzQ=((1.0+w.*w*te(l)*ts(l))./d)+sumzQ;
sumnQ=(w*(te(l)-ts(l))./d)+sumnQ;
end
% Qualitaetsfaktor als Funktion der Frequenz fuer dem L-fachen SLS
q=(1.0-L+sumzQ)./sumnQ;
%Plotting of Q(f,L,f_l,tau) for a generalized SLS
% for plotting
Q0=50.0;
fp1=1.0;fp2=200;df=0.1;
%------------------------------------------------------------------
% Relaxation frequencies and tau od
L=1;
tau1=2.0/Q0;
fl1=[70.0];
%----------------------------------------------------------
f=fp1:df:fp2;
w=2*pi*f;
ts=1.0./(2.0*pi*fl1);
te=ts*(tau1+1.0);
Qft1=qgsls(te,ts,1,w);
;
p_h1=plot(f,Qft1,'k-');
xlabel('frequency [hz]');
ylabel('quality factor (Q)');
axis([fp1 fp2 10 200]);
% Load ASCII-data
load /seismics/data1/SUPGE053/MILET/C/tr94_24c.asc
trace=tr94_24c;
clear tr94_24c;
% Old sample intervall
dt=0.001;
n=size(trace,1);
t=dt:dt:n*dt;
trace=trace/max(trace);
plot(t,trace,'+')
%pause
% new sampling rate (corresponds to FD time step):
dtnew=0.00025;
tnew=dtnew:dtnew:n*dt;
% spline interpolation
trace_new=spline(t,trace,tnew);
hold on
plot(tnew,trace_new,'b')
% write interpolated data into new file:
fid=fopen('/seismics/data1/SUPGE053/FDVISKO/MILET2D/signals/tr94_24c_25.asc','w');
n_new=size(trace_new,2);
for i=1:n_new, fprintf(fid,'%10.5e\n',trace_new(i)); end
fclose(fid);
0.0000 0.0000 1.0000
0.0310 0.0310 0.9966
0.0621 0.0621 0.9931
0.0931 0.0931 0.9897
0.1241 0.1241 0.9862
0.1552 0.1552 0.9828
0.1862 0.1862 0.9793
0.2172 0.2172 0.9759
0.2483 0.2483 0.9724
0.2793 0.2793 0.9690
0.3103 0.3103 0.9655
0.3414 0.3414 0.9621
0.3724 0.3724 0.9586
0.4034 0.4034 0.9552
0.4345 0.4345 0.9517
0.4655 0.4655 0.9483
0.4966 0.4966 0.9448
0.5276 0.5276 0.9414
0.5586 0.5586 0.9379
0.5897 0.5897 0.9345
0.6207 0.6207 0.9310
0.6517 0.6517 0.9276
0.6828 0.6828 0.9241
0.7138 0.7138 0.9207
0.7448 0.7448 0.9172
0.7759 0.7759 0.9138
0.8069 0.8069 0.9103
0.8379 0.8379 0.9069
0.8690 0.8690 0.9034
0.9000 0.9000 0.9000
0.9500 0.9500 0.9500
1.0000 1.0000 1.0000
1.0000 1.0000 1.0000
0.9500 0.9500 0.9500
0.9000 0.9000 0.9000
0.9034 0.8690 0.8690
0.9069 0.8379 0.8379
0.9103 0.8069 0.8069
0.9138 0.7759 0.7759
0.9172 0.7448 0.7448
0.9207 0.7138 0.7138
0.9241 0.6828 0.6828
0.9276 0.6517 0.6517
0.9310 0.6207 0.6207
0.9345 0.5897 0.5897
0.9379 0.5586 0.5586
0.9414 0.5276 0.5276
0.9448 0.4966 0.4966
0.9483 0.4655 0.4655
0.9517 0.4345 0.4345
0.9552 0.4034 0.4034
0.9586 0.3724 0.3724
0.9621 0.3414 0.3414
0.9655 0.3103 0.3103
0.9690 0.2793 0.2793
0.9724 0.2483 0.2483
0.9759 0.2172 0.2172
0.9793 0.1862 0.1862
0.9828 0.1552 0.1552
0.9862 0.1241 0.1241
0.9897 0.0931 0.0931
0.9931 0.0621 0.0621
0.9966 0.0310 0.0310
1.0000 0.0000 0.0000
nt=3334; dt=0.000045;
nrec=459;
%nrec=130;
file_inp1='/workag13/FWI/random_KTB/3D/su_obs/obs_random7_vz_it1.bin.shot2filtLP320';
file_inp3='/workag13/FWI/random_KTB/3D/su/cal_random71_vz_it73.bin.shot2filtLP320';
file_inp2='/workag13/FWI/random_KTB/3D/su/cal_random7_vz_it1.bin.shot2filtLP320';
file_inp39='/workag13/FWI/random_KTB/2D_10sources/su/cal_random723_vz_it73.bin.shot2filtLP320';
file_inp29='/workag13/FWI/random_KTB/2D_10sources/su/cal_random723_vz_it1.bin.shot2filtLP320';
file_inp19='/workag13/FWI/random_KTB/2D_10sources/su_obs/obs_random72a_vz_it1.bin.shot2filtLP320';
fig=55;
%--------------------------------------------------------------------------
SEIS1=binread(file_inp1,nt,nrec);
SEIS2=binread(file_inp2,nt,nrec);
t=dt:dt:nt*dt;
trace=SEIS1(:,3);
size(trace)
figure(1)
plot(t,trace);
figure(fig);
imagesc(1:nrec, t, (SEIS1-SEIS2)/max(max(SEIS1)));
size(SEIS1)
fpal = f_colpal(77, 'linear', 4, 256);
colormap(fpal);
xlabel('receiver number');
ylabel('time in s');
set(get(gca,'Ylabel'),'FontSize',12);
set(get(gca,'Ylabel'),'FontWeight','normal');
set(get(gca,'Xlabel'),'FontSize',12);
set(get(gca,'Xlabel'),'FontWeight','normal');
set(gca,'FontSize',12);
set(gca,'FontWeight','normal');
set(gca,'Linewidth',1.0);
cb = colorbar('vertical');
zlab = get(cb,'ylabel');
set(zlab,'String','residuals normalized to observed data','fontsize',16)
set(gca,'ydir','reverse')
ylim([0.03 0.14]);
caxis([-0.1 0.1])
\ No newline at end of file
%plot of seismograms along two perpendicular receiver lines through shot
%position
clear all;
close all;
file_inp1='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/su/obs_toy_vy_it1.su.shot4';
file_inp2='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/su/obs_toy_vz_it1.su.shot4';
file_inp3='/data14/sdunkl/3DAWAIT/trunk_JURECA/results_toy/su/obs_toy_vy_it1.su.shot4';
%file_inp3='/data14/sdunkl/surface/results6/su/cal_surface5_vz_it70.su.shot10_int';
fignum=17;
shot1=su2matlab(file_inp1);
shot2=su2matlab(file_inp2);
shot3=su2matlab(file_inp3);
dt=shot1(1).dt*1e-6
nt=shot1(1).ns
t=dt:dt:nt*dt;
[shot1(1).sx]
[shot1(1).sy]
numgx=find([shot1.gx] == [shot1(1).sx])
for i=1:13
shot1cutx(i)=shot1(numgx(i));
shot2cutx(i)=shot2(numgx(i));
shot3cutx(i)=shot3(numgx(i));
offset(i)=([shot1cutx(i).gy]-[shot1cutx(i).sy])/100;
end
%max(max(shot1cutx(1).trace))
%max((shot2cutx(1).trace))
figure(fignum)
for i=1:13
plot(t,(shot1cutx(i).trace)/abs(max(shot2cutx(i).trace))+0.01*offset(i),'k-','LineWidth',1);hold on
%plot(t,(shot2cutx(i).trace)/abs(max(shot1cutx(i).trace))+0.05*offset(i),'b-','LineWidth',1);hold on
%plot(t,(shot3cutx(i).trace)/abs(max(shot1cutx(i).trace))+0.03*offset(i),'r-','LineWidth',1);
hold on
end
xlabel('time in s');
ylabel('offset in m');
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(gca,'YTick',0:5:5)