% sembrance_Height3.m
% 2005.09.01 ryoya
% 2006.05.08 tom
% velocityとDirection Angleを決め撃ちする．
clear all;
availchannel=[1 2 3 4 5 6 7 8 9 10 11 12];
tt=0:0.005:49.995;  %書き換えました
factors=0.0003;
%
%
%
fpath='C:\Itoh\across\Horai\Matlab\Transfer_a_rad_9month\';  %書き換えました
fname=[fpath,'Transferfunction01.mat'];%書き換えました
load(fname);
for ii=1:12;
udcomp(:,ii)=data(:,3*(ii-1)+1);
udnoise(:,ii)=noisedata(3*(ii-1)+1);
end
udcomp(:,2)=data(:,6);
udcomp(:,3)=-udcomp(:,3);
udcomp=udcomp(:,1:12);
udnoise=udnoise(:,1:12);
%
%
% 距離計算
path(path,[pwd,'\geometry\']);
set(0,'defaultuicontrolfontsize',14);
set(0,'defaultaxesfontsize',14);
set(0,'defaulttextfontsize',14);
tonolatlonheight=[35+23/60+3/3600 137+13/60+10/3600 100];
S01_15latlonheight=[
    35.0156111 137.6257500 237.800; % S01 ここから下は正しい??　Handy GPSによる値．
    35.0134722 137.6241944 210.800; % S02
    35.0113333 137.6228056 195.000; % S03
    35.0094444 137.6224722 183.100; % S04
    35.0066944 137.6224722 164.300; % S05
    35.0071389 137.6168611 194.000; % S06 
    35.0043889 137.6223889 162.500; % S07 
    35.0022778 137.6228611 165.700; % S08 
    35.0032222 137.6290556 186.900; % S09 
    35.0011111 137.6283889 168.600; % S10 
    34.9998611 137.6238611 154.600; % S11 
    34.9979167 137.6239444 152.000]; % S12
rad = 6378137;
finv= 298.257223563;
tonoxyz=transWGS842ITRF(tonolatlonheight,rad,finv);
S0115xyz=transWGS842ITRF(S01_15latlonheight(availchannel.',:),rad,finv);
Deltalatlonheight=transDelta_ITRF2WGS84(S0115xyz-ones(size(S0115xyz,1),1)*S0115xyz(1,:),S0115xyz(1,:),rad,finv); % lat,lon,height
%
%
%
tlimit=50; fmax=100; 
udnoise(1,:)=udnoise(1,:)./(max(udcomp(:,:),[],1)-min(udcomp(:,:),[],1));
udcomp(:,:)=(udcomp-ones(size(udcomp,1),1)*mean(udcomp,1))./(ones(size(udcomp,1),1)*(max(udcomp,[],1)-min(udcomp,[],1)));
[B,A]=butter(3,[9/fmax,16/fmax]);
udcomp=filter(B,A,udcomp);
ffud=fft(udcomp);
ffudtilda=ffud(2:tlimit*fmax+1,:);
fmin=1/tlimit; fnum=length(2:tlimit*fmax+1); % 0.02-50Hz
fintvl=1/tlimit;
freq0 = (fmin:fintvl:fmin+fintvl*(fnum-1)).';
freq = (0:fintvl:fmax).';
fn = (length(freq)-1)*2; % 全データ点数
%
%
%
%=====================Sembrance========================
anglefromsource=135;       %決めうち！
degmin=10; degmax=80; degint=2; % 2°毎に，80°まで入射角を取る．35こ．
Pvel=4000;                 %決めうち！
twindowstart=9.7; twindowend=26.5;
twindowwidth=0.2; twindowintval=0.05;
ii=0;
for kkk=degmin:degint:degmax;% 15個
    ii=ii+1;
    deg(ii)=(ii-1)*degint+degmin;
    distancefromS01=(Deltalatlonheight(:,1).*cos(anglefromsource/180*pi)+Deltalatlonheight(:,2).*sin(anglefromsource/180*pi))*sin(deg(ii)/180*pi)+Deltalatlonheight(:,3)*cos(deg(ii)/180*pi); % 北から時計回りに．
    difftime = (distancefromS01).'/Pvel; % 波の到来方向を+に採る．1xMベクトル
    %この分の位相をその波形に足すことで，位相を基準点までずらすことが出来る(遅く到達する波を表現できる)．
    difftimes = ones(fnum,1)*difftime; freqs = freq0*ones(1,length(difftime)); diffphase(:,:) = 2*pi*difftimes.* freqs; % 各周波数成分，各観測点でずらす位相を計算．
    diffdata = ffudtilda.*exp(i*diffphase); % 位相ずらした
    fdatas = zeros(fn,size(diffdata,2));
    fdatas(round(fmin/fintvl)+1:round(fmin/fintvl)+fnum,:) = diffdata(1:fnum,:);
    fdatas(fn+1-round(fmin/fintvl):-1:fn+2-round(fmin/fintvl)-fnum,:) = conj(diffdata(1:fnum,:));
    clear diffdata;
    datas = real(ifft(fdatas)*fn/2);
    for jj=1:(twindowend-twindowstart)/twindowintval+1; % 
        twindow =round((twindowstart-twindowwidth/2+(jj-1)*twindowintval)*fmax*2)+1:round((twindowstart+twindowwidth/2+(jj-1)*twindowintval)*fmax*2)+1;% P波部分 0.2s
        twindow2 =round((twindowstart-twindowwidth/2+(jj-1)*twindowintval)*fmax*2)+1:round((twindowstart+twindowwidth/2+(jj-1)*twindowintval)*fmax*2)+1;% P波部分 0.2s
        SE(ii,jj) = sum(sum(datas(round(twindow),:),2).^2,1)/sum(sum(datas(round(twindow2),:).^2,2),1)/size(datas,2);
    end
    fprintf('Incident angle: %3.3d\n',deg(ii));
end

t=twindowstart:twindowintval:twindowend; % P波の場合
ang=degmin:degint:degmax;
figure;
imagesc(t,ang,SE(:,:));set(gca,'FontSize',20,'CLim',[0 0.5]); 
colorbar;
ylabel('入射角 (deg)');
xlabel('時間 (sec)');
title(sprintf('P波速度 4.0 km/s',Pvel));
save Result3DSemblance_Time.mat SE t ang Pvel;

[Y,I2]=max(SE,[],2); % 30配列のうち最も大きな値を並べた行列の，最大の値は？Mx1のベクトル
[Z,I3]=max(Y,[],1); % その中でも最大の値は？1個のスカラー
maxindicator=[I3,I2(I3)]; % 角度のインデックス，タイムウィンドウのインデックス，スローネスのインデックス，
%
%
%     そのパラメータからずらす位相の値を計算
deg=(maxindicator(1)-1)*degint+degmin; 
twindow=round(((maxindicator(2)-1)*twindowintval+twindowstart-twindowwidth)*fmax*2)+1:round(((maxindicator(2)-1)*twindowintval+twindowstart+twindowwidth)*fmax*2)+1;% P波部分
distancefromS01=(Deltalatlonheight(:,1).*cos(anglefromsource/180*pi)+Deltalatlonheight(:,2).*sin(anglefromsource/180*pi))*sin(deg/180*pi)+Deltalatlonheight(:,3)*cos(deg/180*pi); % 北から時計回りに．
difftime = (distancefromS01).'/Pvel; % 波の到来方向を+に採る．1xMベクトル
%この分の位相をその波形に足すことで，位相を基準点までずらすことが出来る(遅く到達する波を表現できる)．
difftimes = ones(fnum,1)*difftime; freqs = freq0*ones(1,length(difftime)); diffphase(:,:) = 2*pi*difftimes.*freqs; % 各周波数成分，各観測点でずらす位相を計算．
%
%
diffdata = ffudtilda.*exp(i*diffphase); % 位相ずらした
fdatas = zeros(fn,size(diffdata,2));
fdatas(round(fmin/fintvl)+1:round(fmin/fintvl)+fnum,:) = diffdata(1:fnum,:);
fdatas(fn+1-round(fmin/fintvl):-1:fn+2-round(fmin/fintvl)-fnum,:) = conj(diffdata(1:fnum,:));
clear diffdata;

tt=0:1/(fmax*2):50-1/(fmax*2);
datas = real(ifft(fdatas)*fn/2);
figure;
for ii=1:12;
    plot(tt,ii+0.5+datas(:,ii)*factors,'r'); hold on;
end
set(gca,'FontSize',20,'YDir','reverse'); xlim([9.7 16.5]);
fprintf('Angle=%5.1f deg, Time window=%4.2f-%4.2fs, Velocity=%5.2f km/s',deg,(twindow(1)-1)/2/fmax,(twindow(end)-1)/2/fmax,Pvel/1000);
