% GEOSENSOR Project, HI 03-01-2004
%
%
% rot_plot_basic
% Program to visualize seismograms and rotation data
% Calculating crosscorrelation and phase velocity as a function of time

clear
close all

%cd ../NewDataBase_rest               % cd to DataBase
list_of_events                  % Get list of available events

% ievent=[2 4 5 6 16 17 18 19];
% ievent=[5 21 16 18];
ievent=[ 16 ];   % FOR PETER ONLY ONE EVENT WORKS

Tmin=30;        % Low pass filtering cut-off period
np=length(ievent);

ip=0;

for iiii=1:np,                     % Loop over events

    
 %   if iii==1, Tmin=3;end
    ie=ievent(iiii);
    
    
    disp(sprintf(' Processing event %i : %s ',ie,deblank(events(ie,:))));
    eval(['load ' deblank(events(ie,:)) ]);
    
    
% Converting translations to m/s, rename
Transverse_accel = Transverse_accel;
Ttime = translation_time;

% Rename Rotations
RZdata = rotation_Z; 
Rtime = rotation_time;

% Rename Trans
Tdata=Transverse_accel;
Tdata(length(Tdata)+1)=0;

% Get length of tiem window and sampling rate
dts=(Ttime(2)-Ttime(1));
nts=length(Ttime);
dtr=(Rtime(2)-Rtime(1));
ntr=length(Rtime);

% Resample rotation rate to the same as translations
iRZdata=interp1(Rtime,RZdata,Ttime);
% Getting rid of NaNs
iRZdata(find(isnan(iRZdata)==1))=(1:length(find(isnan(iRZdata)==1)))*0;


% Filtering
disp(' Filtering ... ')
Tdata=lowpass(Tdata,dts,1/Tmin,1);
iRZdata=lowpass(iRZdata,dts,1/Tmin,1);


% Quantify semblance between data
d=20;               % length of time window (s)

%if iii==1, d=10; end

id=round(d/dts);
id2=round(id/2);
thres=.8; % CROSS CORRELATION THRESHOLD !!!!
c0=max(Tdata)/max(iRZdata);
dc=5000;

i=id2;
ii=1;
iii=1;

disp(' Calculating semblance ');

cc=[];
tc=[];

while i < nts,
    
    ia=i-id2+1;iz=i+id2; 
    if iz<nts,
    t(ii)=Ttime(i);
    x(ii)=max(xcorr(Tdata(ia:iz),iRZdata(ia:iz),'coeff'));
    
    
    % find appropriate phase velocity 
    % Rationale: loop through  possible phase velocities and find the one
    % that minmizes the difference between the scaled Translations and
    % rotations in the whole time window
    
    if x(ii) > thres,
        
        iii=iii+1;
        
        k=0;
        for c=500:5:20000,     
        k=k+1;
        alpha(k)=1/(2*c*1);
        mis=iRZdata(ia:iz)-alpha(k)*Tdata(ia:iz);
        eps(k)=sum(mis'*mis);
        end
        [c,ki]=min(eps); 
        nalpha=alpha(ki);
        tc(iii)=Ttime(i);
        cc(iii)=abs(1/(2*nalpha));
        
    end
    ii=ii+1;
    
    end
    i=i+id;
end

% Extracting window
ia=round(event_windows3(ie,1)/dts);
iz=round(event_windows3(ie,2)/dts);
ia=1;
iz=nts;
ta=Ttime(ia);
te=Ttime(iz);
tmax=max(abs(Tdata(ia:iz)));
rmax=max(abs(iRZdata(ia:iz)));
fs=6;



figure
fs=9;

subplot(211)
hl = plot(Ttime(ia:iz),Tdata(ia:iz),'k-',Ttime(ia:iz),iRZdata(ia:iz)/rmax*tmax,'r-');
axis([ta te -Inf Inf])
ylabel('Transverse Acc. (m/s^2)','FontSize',fs)
set(gca,'FontSize',fs)

subplot(212)
hl = plot(tc,cc/1000,'+');
axis([ta te -Inf Inf])
set(gca,'FontSize',fs)
xlabel(' Time (s) ','FontSize',fs)
ylabel('Velocity (km/s)','FontSize',fs)

text(ta+.01*(te-ta),3.1,sprintf('%s',deblank(eventsc(ie,:))),'Fontsize',fs)

drawnow

pfile=sprintf('fig_phase_%i.jpg',iiii);
eval([' print -djpeg100 ' pfile ]);

end


