%% Infrasound Array Processing: Infrasound Array Beam Forming
% Script written by Ezer and Brando for NM Tech Course
% June 24, 2009
% Script based on Almendros et al. 2002 "Frequency-slowness analysis"
 
%% Strategy
 
% 1. Identify all pulses of interest and pick arrival times
% 2. Correlate pulses across array
% 3. Stack pulses 
% 4. Highest amplitude in stack = azimuth
% 4. Comparison of relative energies
 
%% Notes on Frequency-slowness analysis
 
% slowness = 1/(apparent velocity of P wave)
% Our frequency band of interest = 0.2 to 10 Hz
% uses overlapping frequency bands
% window = 2.56 seconds (100 samples)
% overlap = 0.2 s
% dominant peak in slowness spectrum for a given freq band rep the most
% coherent component of the wavefield
 
%% READ IN WIGGLES FOR ALL PICKS ALL CHANNELS, WINDOW, FILTER, PLOT
%
clear all;
close all;
cd '/Users/geop572/Desktop/SEGY/'
load Tpicks_Compiled.mat
% load Tpicks_final.mat % This file contains a x vector with start times of picked events
Station_Name = 'TINF';
das = 'A988';

for g = 1:length(X_final)
    

x_convert = X_final(1,g)/60 - 0.02 % converts decimal seconds to minutes
%y = zeros(1,length(X_final));
y=0;
start_time_picks = [Start_time_final(g,1) Start_time_final(g,2) x_convert y];
window = 5; %gets start time n s sooner than picked event. NOTE = picks done on channel 6 which is furthest channel from vent
k=1;
for i = 1 : length(x_convert)
    for j=1:6 % channels 1-6 but channel 4 doesn't exist
        if j==4
            wiggle_final(:,j,k)=zeros(length(wiggle_final),1);
        else
        [wiggle(:,i),hdr] = segycut(start_time_picks(i,:),window,das,j); % Picks where done on Channel 6
        wiggle_final(:,j,k)=wiggle(:,i);
        hdr_final(i)=hdr;
        end
    end
    
    
    sps = hdr.sps;
    ts = (1:length(wiggle))/sps; % create time series vector
    infra_sens = 1*248.843/(0.01*7/12); 
    % Sensitivity of Jeff's sensor = 23 microVolts/ Pa | 7 Voltage control
    % oscillator
    infra_gain = hdr.gain; % programmed gain setting
    atodtoPa= wiggle_final*infra_sens*hdr.atod/32; %atodtoPa analog_to_dig_to_Pascals
    for j=1:6
    wiggle_detrend(:,j,k)=detrend(atodtoPa(:,j,k));
    end
    k=k+1; 
    
% Butter filter data from 0.2 to 10 Hz
npoles = 2;
low = 10; % corner frequency at 10 Hz
high = 0.2; % corner frequency at 0.2 Hz
% signal processing tool box with butter is necessary for following operations
% first do bandpass
 
[B1,A1] = butter(npoles,[high low]/(sps/2)); % The order of high and low MATTERS. lowest number must come first.
% now do highpass
 
bandpass_wiggle = filter(B1,A1,wiggle_detrend);
bandpass_wigglef=bandpass_wiggle(:,1:3,:);
bandpass_wigglef(:,4,:)=bandpass_wiggle(:,5,:);
bandpass_wigglef(:,5,:)=bandpass_wiggle(:,6,:);
 
% To check all data was read in: size(bandpass_wigglef(samples,channels,
% picks))
 
end
% Plot to check to make sure all data read in
% One plot for each pick. Each plot contains all five channels
% This only plots the first 3 picks as a QC
% for i=1:5
%     figure(1), clf
%     subplot(length(start_time_picks(1:3,1)),1,i),plot(bandpass_wigglef(:,:,i))
%     ylabel('Pa')
%     xlabel('time')
%    
% end
%  
figure(1);
    
for i = 1:5
    
    plot(bandpass_wigglef(:,i));
    hold on
    xlabel('time (sec)')
    ylabel('Pa')
    title('Infrasound days 161 to 167')
end
end

%% Crate Slowness and Power Matrices
 
% Create slowness matrix
%vmin=- 360/1000;
%vmax= +360/1000;
%vx=vmin:0.01:vmax;
%vy=vx;
%Sx=1./vx;
%Sy=1./vy;
Smin = -1/0.01;
Smax = 1/0.01;
incr = 10;
Sx = Smin:incr:Smax;
Sy = Sx;
% Where Sx and Sy are slowness (or 1/(apparent velocity))
% Sx and Sy are in km
 
% Calculate power for each value of slowness (Sx,Sy)
power = zeros(length(Sx),length(Sy));
%Creates a matrix of zeroes to be filled later in code
 
%% NETWORK GEOMETRY

infrasound_x = [0, -sind(45)*22.8/1000, -cosd(60)*27.9/1000, 24.9/1000, 0];
infrasound_y = [0, (-cosd(45)*22.8)/1000, (sind(60)*27.9)/1000, 0, -23.3/1000];

% plot(infrasound_x, infrasound_y, 'o')

% Where vectors in order of channel (1,2,3,5,6)
% Where channel 1 is at heart of array
% Why *1000? everything needs to be in km



%% CIRCSHIFT through SLOWNESS MATRIX for each pick, 
% stack channels, calc power
% This takes ~1 minutes for individual pick and mucho minutes for all 30 picks
 
% Goal here is to square waveform and integrate it (or take the sum)
bandpass_wigglef_ch2 = bandpass_wigglef(:,2);
bandpass_wigglef_ch3 = bandpass_wigglef(:,3);
bandpass_wigglef_ch5 = bandpass_wigglef(:,4);
bandpass_wigglef_ch6 = bandpass_wigglef(:,5);

for l = 1:length(Sx) % length Sx = 41
    for m = 1:length(Sy) % length Sy = 41
        shifted_wiggle_channel2 = circshift(bandpass_wigglef_ch2, round(100*(-(Sx(l).*infrasound_x(2)-(Sy(m).*infrasound_y(2))))));
        shifted_wiggle_channel3 = circshift(bandpass_wigglef_ch3, round(100*(-(Sx(l).*infrasound_x(3)-(Sy(m).*infrasound_y(3))))));
        shifted_wiggle_channel5 = circshift(bandpass_wigglef_ch5, round(100*(-(Sx(l).*infrasound_x(4)-(Sy(m).*infrasound_y(4))))));
        shifted_wiggle_channel6 = circshift(bandpass_wigglef_ch6, round(100*(-(Sx(l).*infrasound_x(5)-(Sy(m).*infrasound_y(5))))));
        Tstack= bandpass_wigglef(:,1) + shifted_wiggle_channel2 + shifted_wiggle_channel3 + shifted_wiggle_channel5 + shifted_wiggle_channel6;
        %Tstack = channel 1 + channel 2 + channel 3 + channel 5 + channel 6
        %(NOTE: channel 4 does not exist)
        
        power(m,l) = sum(Tstack.^2);
        % power score = sum of number of samples * [(sum over channels of
        % shifted waveform)^2]
        %end
    end
end



%%
C = max(max(power));
% First max finds max for each row and outputs a column
% Second max finds max for resultant column and outputs a single number
% The max power represents the ideal shift of the waveform across array
[X,Y,C] = find(power==C);
% For C, this function ravages the power matrix and outputs the coordinates
% of C
Sxmax = Sx(X);
Symax = Sy(Y);
 
azimuth = atand(Symax/Sxmax);
% atand = arctangent convert to degrees
% stand(Symac/Sxmax) gives you the highest amplitude in stack which is the azimuth
velocity = 1/(sqrt((Symax^2)+(Sxmax)^2));
 
norm_power = power/max(power(:));
% power(:) converts power matrix to one column

filename = ['model_' num2str(g)]
save(filename)
%save 'model_',num2str(g),'.mat'
%save model_Compiled.mat

end
close all

%% 

cd '/Users/geop572/Desktop/SEGY/'

Sy_final = zeros(length(X_final),length(Sx));
Sx_final = zeros(length(X_final),length(Sx));
norm_power_final = zeros((length(X_final)),length(Sx),length(Sx));
Sxmax_final = zeros(1,length(X_final));
Symax_final = zeros(1,length(X_final));
azimuth_final = zeros(1,length(X_final));
velocity_final = zeros(1,length(X_final));
% These create matrices filled with zeros to be populated with values below

for h = 1:length(X_final)

tt = strcat(['load model_',num2str(h),'.mat']);
eval([tt])

Sy_final(h,:) = Sy;
Sx_final(h,:)  = Sx;
norm_power_final(h,:,:) = norm_power;
Sxmax_final(h,:)  = Sxmax;
Symax_final(h,:)  = Symax;
azimuth_final(h,:)  = azimuth;
velocity_final(h,:)  = velocity;
end

save Model_Compiled.mat Sy_final Sx_final norm_power_final Symax_final Sxmax_final azimuth_final velocity_final

%%
% figure (2), clf
% 
% contour(Sy,Sx,norm_power,[.95:.003:.999]); %Sx is 1x41, Sy 1x41, norm_power 41x41; [] points
% %contour(Sx,Sy,norm_power); %Sx is 1x41, Sy 1x41, norm_power 41x41; [] points
% xlabel('S_x(s/km)');
% ylabel('S_y(s/km)');
% axis equal
% axis square
% z = 30;
% xlim([-z,z])
% ylim([-z,z])
%  
% xtick=-z:5:z;
% ytick=-z:5:z;
% set(gca,'XTick',xtick,'YTick',ytick)
% grid on
% gtext(['Velocity: ' num2str(velocity) 'km/s'])
% gtext(['Azimuth: ' num2str(azimuth) '^o'])