

%
% this program reads in the real data (rotation rate around z axis and transverse acceleration 
% for the  M8.3 25 Sep 2003 event in Hokkaido (it s actually called Tokachi-Oki)
% Heiner Igel, Dec. 09

clear 
close all


% Tasks for ROSE scholars:

% 1. Uncomment lines 71,72 to see the effects of using the zero-lag
% coefficient what is the difference?

% 2. Change the length of the sliding time window d to larger and smaller
% values and look at the results (line 56)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REAL DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load hokkaido_rotations_observed
rot_obs_z=rotation_Z;              % Rotation rate in  rad/s
rot_obs_time=rotation_time*3600;   % Time in seconds

load hokkaido_translations_observed
tra_obs_time=nTime;                 % time for transverse acceleration
tra_obs_t=nTdata;                   % tra for transverse acceleration

% plot real data
c=2*4400;      % scaled with phase velocity
figure
plot(rot_obs_time,rot_obs_z*c,'r',tra_obs_time,tra_obs_t,'k-'), axis([3500 7000 -Inf Inf])
title(' REAL DATA: Rotation - red, Acceleration - black ')
ylabel(' Acceleration (m/s^2) ')
xlabel(' Time (s) ')
pause 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Resample to sampling rate of accelerationr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rrot_obs_z=interp1(rot_obs_time,rot_obs_z,tra_obs_time);
rrot_obs_time=tra_obs_time;
figure
plot(rrot_obs_time,rrot_obs_z*c,'r',tra_obs_time,tra_obs_t,'k-'), axis([3500 7000 -Inf Inf])
title(' RESAMPLED: Rotation - red, Acceleration - black ')
ylabel(' Acceleration (m/s^2) ')
xlabel(' Time (s) ')
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cross correlation along time axis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d=10; % length of time window (s)
dt=tra_obs_time(2)-tra_obs_time(1);
id=round(d/dt);
id2=round(id/2);
nts=length(tra_obs_time);

i=id2;
ii=1;
disp(' Calculating semblance ');
while i < nts,
    ia=i-id2+1;ie=i+id2;
    if ie<nts,
    t(ii)=tra_obs_time(i);
    x(ii)=max(xcorr(tra_obs_t(ia:ie),rrot_obs_z(ia:ie),'coeff'));
    % take zero lag value
    %tmp=xcorr(tra_obs_t(ia:ie),rrot_obs_z(ia:ie),'coeff');
    %x(ii)=tmp(length(tra_obs_t(ia:ie)));
    ii=ii+1;
    end
    i=i+id;
end


figure
% plot results
subplot(211),plot(rrot_obs_time,rrot_obs_z*c,'r',tra_obs_time,tra_obs_t,'k-'), axis([3500 7000 -Inf Inf])
xlabel(' Time (s)')
ylabel(' Acceleration (m/s^2) ')

subplot(212),plot(t,x)
axis([3500 7000 0 1])
xlabel(' Time (s)')
ylabel(' Max Cross-corr. norm. ')


