Commit 36b61444 by Simone Butzer

### add Matlab plot programs to mfiles

parent 29d95e46
 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]);
mfiles/qapprox.m 0 → 100644
 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]);
mfiles/qflt.m 0 → 100644
 % 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);
mfiles/qgsls.m 0 → 100644
 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;
mfiles/qplot.m 0 → 100644
 %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]);
mfiles/resamp.m 0 → 100644
 % 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
mfiles/seismo.m 0 → 100644