head	3.0;access;symbols;locks; strict;comment	@// @;3.0date	2000.06.13.19.17.59;	author gilles;	state Exp;branches;next	;desc@Release 3@3.0log@*** empty log message ***@text@function [datatot]=qc_dsi(datain,headw1,tstart,tend,comp1,comp2)%% Please note: Software does not check for reasonable parameters or dead traces%%%by G. Perron (Nov 15th, 1996)%%based on rot3c_dirp from S. Guest and D. Eaton+%$Id:$%$Log:$  %%Copyright (C) 1998 Seismology and Electromagnet+ic Section/%Continental Geosciences Division/Geological Survey of Canada%%This library is free software; you can redistribute it and/or%modify it under the terms of the GNU Library General Public%License as published by the Free Software Foundation; either%version 2 of the License, or (at your option) any later version.%%This library is distributed in the hope that it will be useful,%but WITHOUT ANY WARRANTY; without even the implied warranty of%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU%Library General Public License for more details.%%You should have received a copy of the GNU Library General Public%License along with this library; if not, write to the%Free Software Foundation, Inc., 59 Temple Place - Suite 330,%Boston, MA  02111-1307, USA.%%DSI Consortium%Continental Geosciences Division%Geological Survey of Canada%615 Booth St.%Ottawa, Ontario%K1A 0E9%%email: dsi@@cg.nrcan.gc.cadatarot=rot3c(datain,headw1,tstart,tend,comp1,comp2);disp('angle rot:');datarot.th{1}(5,:)disp('S2R azimuth:');datarot.th{1}(9,:)for i=1:datarot.fh{13}     if (datarot.th{1}(9,i) <= 180)        datarot.th{1}(9,i)=360-(180-datarot.th{1}(9,i));        datarot.th{2}(9,i)=datarot.th{1}(9,i);        datarot.th{3}(9,i)=datarot.th{1}(9,i);     else        datarot.th{1}(9,i)=(datarot.th{1}(9,i)-180);        datarot.th{2}(9,i)=datarot.th{1}(9,i);        datarot.th{3}(9,i)=datarot.th{1}(9,i);      end	        if datarot.th{1}(9,i)-datarot.th{1}(5,i)<=0        h1pos(i)=(360-datarot.th{1}(5,i))+datarot.th{1}(9,i);      else        h1pos(i)=datarot.th{1}(9,i)-datarot.th{1}(5,i);   endendsroff=datarot.th{1}(53,:);%sroff=sqrt(((datarot.th{1}(29,:)-datarot.th{1}(35,:)).^2) + ((datarot.th{1}(31,:)-datarot.th{1}(37,:)).^2));sraz=datain.th{1}(9,:);relev=datarot.th{1}(39,:);seisdata=datain;disp('[dataout]=qc_dsi(datain,headw1,tstart,tend,comp1,comp2)');w=pi/180;trclength=datain.fh{7}; %number of points per tracesmpint=datain.fh{8}; %smpint is the sampling interval%dataout=datain;%check to make sure data is separated into components%ntr is the number of traces in each recordfor COUNT=3:-1:1 %get number of traces in each component	ntr(COUNT)=datain.th{COUNT}(12,1);end %forif (ntr(1)~=ntr(2)) | (ntr(1)~=ntr(3))	error('check data format - different number of traces in components');end%ifif length(datain.dat)~=3	error('data must have only 3 records - one for each of x, y and z');end %if%*************************************************************************%create a look-up table of sin and cos valuesfor COUNT=0:360   cosang(COUNT+1)=cos(COUNT*w);   sinang(COUNT+1)=sin(COUNT*w);end %for COUNTntr=ntr(1);angmx(1:ntr)=0; %initialize angmx vector for storing rotation anglessheet=qcsheet;%hold onhr=findobj(sheet,'Tag','rotations');ha=findobj(sheet,'Tag','amplitudes');hs{1}=findobj(sheet,'Tag','seis1');hs{2}=findobj(sheet,'Tag','seis2');hs{3}=findobj(sheet,'Tag','seis3');hs{4}=findobj(sheet,'Tag','seis4');hs{5}=findobj(sheet,'Tag','seis5');hpol=findobj(sheet,'Tag','Polar1');a={'-r','-g','-b','-m','-k','-c'};b={'*r','*g','*b','*m','*k','*c'};c={'-*r','-*g','-*b','-*m','-*k','-*c'};ledg=num2str(datain.th{1}(26,:)');seisdata=sortrec_new(seisdata,13);datarot=sortrec_new(datarot,13); datatot=seisdata;for i=1:datatot.fh{12}   datatot.dat{i}(:,4:6)=datarot.dat{i}(:,:);   datatot.th{i}(:,4:6)=datarot.th{i}(:,:);   datatot.th{i}(12,:)=datatot.th{i}(12,:)+3;enddatatot.fh{1}=datatot.fh{1}*2;datatot.fh{13}=datatot.fh{13}+3;for COUNT1=1:ntr samp1=round((datain.th{1}(headw1,COUNT1)-datain.fh{9})/smpint)-round(tstart/smpint)+1; %calulates start of interval to be analyzed samp2=round((datain.th{1}(headw1,COUNT1)-datain.fh{9})/smpint)+round(tend/smpint)+1; %calulates end of interval to be analyzed %this following loops over specified angles and returns the angle within %those specified that maximizes the radial component cmax=0; %initializing the max. energy of a component datawin=samp1:samp2; %the window over which the data is to be analyzed  lendatawin=length(datawin); %length of the data window xy=[datain.dat{comp1}(datawin,COUNT1), datain.dat{comp2}(datawin,COUNT1)]; %the window of data to be analyzed z=datain.dat{3}(datawin,COUNT1);  tr1=1; tr2=6; t1=(samp1*seisdata.fh{8}); t2=(samp2*seisdata.fh{8}); smp=seisdata.fh{8}; datarotin=subset(datatot,tr1,tr2,t1,t2,COUNT1,COUNT1);   vert=sum(z.*z); horiz1=sum(xy(:,1).*xy(:,1)); horiz2=sum(xy(:,2).*xy(:,2)); horiztot=horiz1+horiz2; ratio(COUNT1)=sroff(COUNT1)/relev(COUNT1); ratio_obs(COUNT1)=horiztot/vert; allamps=[horiz1 horiz2 vert horiztot]; rt=zeros(lendatawin,2);  axes(hs{COUNT1}); hold on plot(datarotin.th{1}(15,:),'.r'); seisplot2(datarotin.dat{1},t1,t2,1,6,smp,1,0,1,a{COUNT1},t1); for COUNT2=0:359 %checks from 0 to 90 degrees  %the following applies the rotation matrix and sums over each component  rt(:,1)=xy(:,1)*cosang(COUNT2+1)+xy(:,2)*sinang(COUNT2+1);  rt(:,2)=-xy(:,1)*sinang(COUNT2+1)+xy(:,2)*cosang(COUNT2+1);  c1rms=sum(rt(:,1).*rt(:,1));  c1rms2(COUNT2+1)=c1rms;  c2rms=sum(rt(:,2).*rt(:,2));  c2rms2(COUNT2+1)=c2rms;  ang(COUNT2+1)=COUNT2;end %for COUNT2ang=datain.th{1}(9,COUNT1)-ang;for i=1:length(ang)if ang(i) <= 0   ang(i)=360+ang(i);endendl=find(c1rms2==max(c1rms2));w=[ang' c1rms2'];w=sortrows(w,1);axes(hr)hold onsemilogy(w(:,1),w(:,2),a{COUNT1});xlabel('Azimuth of H2')ylabel('RMS Amplitude')box ongrid on axes(ha)hold onsemilogy(allamps,c{COUNT1})hold onxlabel('H1 H2 Z Htot')ylabel('RMS Amplitude')box ongrid on end %for COUNT1axes(hr)legend(ledg(1),ledg(2),ledg(3),ledg(4),ledg(5));%legend(ledg)%hold offh1poscompinfo=[h1pos' round(sroff)' [1 2 3 4 5]' round(sraz)'];cinfosort=sortrows(compinfo,2);fcisort=flipud(cinfosort);for i=1:datarot.fh{12}   axes(hpol);   newplot;   fcisort(i,1)=fcisort(i,1).*pi/180;   [x,y]=pol2cart(fcisort(i,1),fcisort(i,2));   theta=fcisort(i,4)*(pi/180)+pi;   rho=fcisort(i,2);   compass(x,y,a{fcisort(i,3)});   hold on   polar(theta,rho,b{fcisort(i,3)});   hold onendset(gca,'view',[90 -90]);h=num2str(datain.th{1}(56,1));suptitle(['Level ',h])%hratio=figure(2);%plot(ratio,'-*r')%hold on%plot(ratio_obs,'-*b')%xlabel('Trace Number')%ylabel('HSUM/VERT')%grid on%legend(gca,'Theorical Amplitude Ratio','Observed Amplitude Ratio');@