MainInput.m

% This version intentionally follows the original processing logic:
%   1. Use EW1 files as station references.
%   2. Read station metadata from EW1 headers.
%   3. Use EW1 sampling rate and scale factor for all three components.
%   4. Apply linear detrending to EW and UD only.
%   5. Align records to the earthquake origin time.
%   6. Allow matrices to expand dynamically to the longest trace.
%   7. Remove zero and abnormal stations using the original criteria.
% =========================================================================
% Stroing-Motion Data could be downloaded from NIED:
%                 https://www.kyoshin.bosai.go.jp/en/
% You need to register an account and set up an email address:
% e.g.,
% UserName: Yuxiao
% Password: ***
% =========================================================================
% For KiK-net, each station has two sets of three-component instruments:
% e.g.,
% EW1: East-west acceleration record from borehole / underground
% EW2: East-west acceleration record from surface / terrestrial
% Use underground station waveforms to minimize the impact of near-surface
% site effects and urban noise on inversion
% =========================================================================
% The earthquake event time and the recording times at each station should 
% be in JST format (not UTC)
clc;
close all;
clear;
%% =========================================================================
% Input parameters
% File: Input
%         |-- Datafunc: processing functions
%         |-- DataWave: initial waves (e.g., EW1,NS1)
%         |-- OutSta: output station info and location
%         |-- OutWave: output wave.dat (each has 3 components)
%         |-- Info.txt: EQ Information
% =========================================================================
PATH = 'D:\Matlab\MyMechanismFunc\FaultInverseAuto\Input\'; % Root Input
DirWaveIn = [PATH,'DataWave'];
DirWaveOut= [PATH,'OutWave'];
DirStaOut = [PATH,'OutSta'];
DirFunc   = [PATH,'Datafunc'];
DirInfo   = [PATH,'Info.txt'];
addpath(DirFunc);
srate = 4;                      % Target sampling rate in Hz
UnitGal = 1e-2;    % 1 gal = 0.01 m/s^2
fprintf('[Now] %s\n',datetime('now'));
fprintf('[1] Folders are ready \n');

%% =========================================================================
% Read earthquake origin time and source parameters
% =========================================================================
if ~exist(DirInfo, 'file')
    error('Info.txt was not found: %s', DirInfo);
end
fid = fopen(DirInfo, 'r');
if fid < 0
    error('Unable to open Info.txt.');
end
% Read the first description line.
fgetl(fid);
% Read the earthquake origin time.
dateori = strtrim(fgetl(fid));
% Read and discard the parameter-header line.
fgetl(fid);
% Read numerical source parameters.
Strings = fscanf(fid, '%f');
fclose(fid);
if mod(numel(Strings), 11) ~= 0
    error('The number of source parameters is not a multiple of 11.');
end
Strings = reshape(Strings, 11, [])';
% =========================================================================
% IMPORTANT COORDINATE CONVENTION
% The original code uses:
%   epi = [ear_all(:,2), ear_all(:,1)]
%
% The da_zh function expects:
%   loca = [latitude, longitude]
%   epi  = [latitude, longitude]
% =========================================================================
epicenter = [Strings(:, 2), Strings(:, 1)];     % [Latitude, Longitude]
depth = Strings(:, 3);                    % Earthquake depth in km
macha1 = Strings(:, 4:6);                 % First focal mechanism
macha2 = Strings(:, 7:9);                 % Second focal mechanism
fprintf('[2] Earthquake-Info are ready \n');

%% =========================================================================
% Find the downloaded strong-motion event DirWave or just given simply
% =========================================================================
% file_SM = dir(DirInput);
% isSubFolder = [file_SM.isdir];
% folderNames = {file_SM(isSubFolder).name}';
% 
% folderNames(ismember(folderNames, {'.', '..'})) = [];
% 
% if isempty(folderNames)
%     error('No event DirWave was found under: %s', DirInput);
% end
% 
% % Keep the original logic: use the first event DirWave.
% Dname = folderNames{1};
% DirWave = fullfile(DirInput, Dname);
% if ~isfolder(DirWave)
%     error('Event DirWave does not exist: %s', DirWave);
% end
% fprintf('# Successfully find the Info.txt');

%% =========================================================================
% Find KiK-net borehole waveform files (index=1) or surface (index=2)
% =========================================================================
fileEW = dir(fullfile(DirWaveIn, '*.EW1'));
fileNS = dir(fullfile(DirWaveIn, '*.NS1'));
fileUD = dir(fullfile(DirWaveIn, '*.UD1'));
if isempty(fileEW)
    error('No EW1 borehole waveform files were found in: %s', DirWaveIn);
end
if numel(fileEW) ~= numel(fileNS) || numel(fileEW) ~= numel(fileUD)
    warning(['The numbers of EW1, NS1, and UD1 files are inconsistent. ', ...
             'The code assumes matching file order and station count.']);
end
Nsta = numel(fileEW);
fprintf('[3] Stations are ready \n');
fprintf('     Counts: %d \n',Nsta);

%% =========================================================================
% Initialize variables
% MATLAB dynamically expands them when longer records are assigned.
% Therefore, the final sample count equals the longest processed trace.
% =========================================================================
loca = zeros(Nsta, 2);      % [Latitude, Longitude]
mm = repmat(' ', Nsta, 1);  % Station names
ts = zeros(1, Nsta);        % Record-start time minus origin time
EW = zeros(1, Nsta);
NS = zeros(1, Nsta);
UD = zeros(1, Nsta);
ob = zeros(1, Nsta * 3);
% =========================================================================
% Read three-component strong-motion records
% =========================================================================
for i = 1:Nsta
    ewFile = fullfile(DirWaveIn, fileEW(i).name);
    % ---------------------------------------------------------------------
    % Read station metadata from the EW file header
    % ---------------------------------------------------------------------
    fid = fopen(ewFile, 'r');
    if fid < 0
        warning('Unable to open EW file: %s', ewFile);
        continue;
    end
    % Read the station-name line.
    for k = 1:6
        ch = fgets(fid);
    end
    dex = find(ch == 'e', 1, 'first');
    if isempty(dex)
        stationName = strtrim(ch);
    else
        stationName = ch(dex + 1:end);
        stationName(isspace(stationName)) = [];
    end
    mm(i, 1:numel(stationName)) = stationName;

    % Read station latitude.
    ch = fgets(fid);
    dex = find(ch == 'L', 1, 'first');
    if isempty(dex)
        error('Unable to parse latitude from: %s', ewFile);
    end
    locaLat = str2double(ch(dex + 4:end));
    % Read station longitude.
    ch = fgets(fid);
    dex = find(ch == 'L', 1, 'first');
    if isempty(dex)
        error('Unable to parse longitude from: %s', ewFile);
    end
    locaLon = str2double(ch(dex + 5:end));
    loca(i, :) = [locaLat, locaLon];
    % Read the sampling frequency
    for k = 1:3
        ch = fgets(fid);
    end
    dex = strfind(ch, ')');
    if isempty(dex)
        error('Unable to parse sampling frequency from: %s', ewFile);
    end
    srate0 = str2double(ch(dex + 2:end - 3));
    % Read the scale factor
    % Notice: Here, it is simplified to assume that the amplitude scaling
    % factor for each individual station should be the same; however, 
    % a more rigorous approach would involve reading the scaling factor 
    % for each station individually
    % Format:   Scale Factor      2940(gal)/6170270
    for k = 1:3
        ch = fgets(fid);
    end
    dex1 = strfind(ch, 'r');
    dex2 = strfind(ch, '(');
    if isempty(dex1) || isempty(dex2)
        error('Unable to parse scale factor from: %s', ewFile);
    end
    faca = str2double(ch(dex1 + 7:dex2 - 1));
    facb = str2double(ch(dex2 + 6:end));
    % Read the record-start time.
    for k = 1:2
        ch = fgets(fid);
    end
    dex1 = strfind(ch, 'n');
    if isempty(dex1)
        error('Unable to parse record time from: %s', ewFile);
    end
    timestr = ch(dex1 + 4:end);
    daterec = datenum(timestr) * 24 * 3600;
    ts(i) = daterec - datenum(dateori) * 24 * 3600;
    % Skip the final header line before reading waveform samples.
    fgets(fid);
    % ---------------------------------------------------------------------
    % Read EW component
    % ---------------------------------------------------------------------
    ob0 = fscanf(fid, '%f');
    % Convert digital counts using the EW scale factor.
    ob0 = ob0 * faca / facb;
    % Remove the mean value.
    ob0 = ob0 - mean(ob0);
    % Remove linear drift from EW
    obpoly = legendre_fit(ob0, 1:numel(ob0), 1, 0, 2);
    ob0 = ob0 - obpoly;
    % Resample to the target sampling rate.
    ob0 = resample(ob0, srate, srate0);
    % Align the waveform beginning to the earthquake origin time.
    if ts(i) > 0
        ob0 = [zeros(round(ts(i) * srate), 1); ob0];
    else
        ob0 = ob0(1 - round(ts(i) * srate):end);
    end
    obEW = ob0;
    % Convert gal to m/s^2.
    EW(1:numel(obEW), i) = obEW * UnitGal;
    fclose(fid);

    % ---------------------------------------------------------------------
    % Read NS component
    %
    % The original code uses the EW sampling rate and EW scale factor.
    % That behavior is retained here for compatibility.
    % ---------------------------------------------------------------------
    nsFile = fullfile(DirWaveIn, fileNS(i).name);
    fid1 = fopen(nsFile, 'r');
    if fid1 < 0
        warning('Unable to open NS file: %s', nsFile);
        continue;
    end
    % Skip the 17-line NIED ASCII header.
    for k = 1:17
        fgets(fid1);
    end
    ob0 = fscanf(fid1, '%f');
    % Keep the original logic: use EW scale factor.
    ob0 = ob0 * faca / facb;
    % Remove the mean value.
    ob0 = ob0 - mean(ob0);
    % Keep the original logic: no linear detrending for NS.
    ob0 = resample(ob0, srate, srate0);
    % Align the waveform beginning to the earthquake origin time.
    if ts(i) > 0
        ob0 = [zeros(round(ts(i) * srate), 1); ob0];
    else
        ob0 = ob0(1 - round(ts(i) * srate):end);
    end
    obNS = ob0;
    % Convert gal to m/s^2.
    NS(1:numel(obNS), i) = obNS * UnitGal;
    fclose(fid1);

    % ---------------------------------------------------------------------
    % Read UD component
    %
    % The original code uses the EW sampling rate and EW scale factor.
    % That behavior is retained here for compatibility.
    % ---------------------------------------------------------------------
    udFile = fullfile(DirWaveIn, fileUD(i).name);
    fid2 = fopen(udFile, 'r');
    if fid2 < 0
        warning('Unable to open UD file: %s', udFile);
        continue;
    end
    % Skip the 17-line NIED ASCII header.
    for k = 1:17
        fgets(fid2);
    end
    ob0 = fscanf(fid2, '%f');
    % Keep the original logic: use EW scale factor.
    ob0 = ob0 * faca / facb;
    % Remove the mean value.
    ob0 = ob0 - mean(ob0);
    % Remove linear drift from UD, following the original code.
    obpoly = legendre_fit(ob0, 1:numel(ob0), 1, 0, 2);
    ob0 = ob0 - obpoly;
    % Resample to the target sampling rate.
    ob0 = resample(ob0, srate, srate0);
    % Align the waveform beginning to the earthquake origin time.
    if ts(i) > 0
        ob0 = [zeros(round(ts(i) * srate), 1); ob0];
    else
        ob0 = ob0(1 - round(ts(i) * srate):end);
    end
    obUD = ob0;
    % Convert gal to m/s^2.
    UD(1:numel(obUD), i) = obUD * UnitGal;
    fclose(fid2);
    
    % ---------------------------------------------------------------------
    % Construct the observed data matrix
    %
    % Column order:
    % EW1 NS1 UD1 EW2 NS2 UD2 ...
    % ---------------------------------------------------------------------
    if numel(obEW) ~= numel(obNS) || numel(obEW) ~= numel(obUD)
        error(['The three components at station %s do not have identical ', ...
               'lengths after processing.'], strtrim(mm(i, :)));
    end
    ob0 = [obEW, obNS, obUD] * UnitGal;
    ob(1:numel(obEW), 3 * (i - 1) + 1:3 * i) = ob0;
    fprintf('   Processed %d/%d: %s;   Samples=%d\n', ...
        i, Nsta, strtrim(mm(i, :)), numel(obEW));
end
fprintf('[4] Successfully loaded observations (Three-components) \n');
clear obEW obNS obUD ob0 obpoly;

%% =========================================================================
% Remove stations whose EW component sum is zero
% =========================================================================
ba = sum(EW, 1);
ca = find(ba(:) == 0);
a = ca';
for i = 1:numel(a)
    ob(:, (a(i) * 3 - 2):a(i) * 3) = 0;
    EW(:, a(i)) = 0;
    NS(:, a(i)) = 0;
    UD(:, a(i)) = 0;
    loca(a(i), :) = 0;
    mm(a(i), :) = 0;
end
ob(:, sum(abs(ob), 1) == 0) = [];
EW(:, sum(abs(EW), 1) == 0) = [];
NS(:, sum(abs(NS), 1) == 0) = [];
UD(:, sum(abs(UD), 1) == 0) = [];
loca(sum(abs(loca), 2) == 0, :) = [];
mm(sum(abs(mm), 2) == 0, :) = [];

% =========================================================================
% Remove stations with abnormal differences among components
%
% This section intentionally follows the original calculation:
% max(EW), max(NS), and max(UD) are used without abs().
% =========================================================================
mEW = max(EW);
mNS = max(NS);
mUD = max(UD);
s1 = mEW ./ mNS;
s2 = mNS ./ mEW;
[a1, b1] = max(max(s1, s2));
c1 = mEW ./ mUD;
c2 = mUD ./ mEW;
[a2, b2] = max(max(c1, c2));

v1 = mNS ./ mUD;
v2 = mUD ./ mNS;
[a3, b3] = max(max(v1, v2));
s1 = [a1, a2, a3];
s2 = [b1, b2, b3];

dex = find(s1 > 10);
ndex = s2(dex);
ndex = unique(ndex);

a = ndex;
for i = 1:numel(a)
    ob(:, (a(i) * 3 - 2):a(i) * 3) = 0;
    EW(:, a(i)) = 0;
    NS(:, a(i)) = 0;
    UD(:, a(i)) = 0;
    loca(a(i), :) = 0;
    mm(a(i), :) = 0;
end
ob(:, sum(abs(ob), 1) == 0) = [];
EW(:, sum(abs(EW), 1) == 0) = [];
NS(:, sum(abs(NS), 1) == 0) = [];
UD(:, sum(abs(UD), 1) == 0) = [];

loca(sum(abs(loca), 2) == 0, :) = [];
mm(sum(abs(mm), 2) == 0, :) = [];
clear a a1 a2 a3 b1 b2 b3 c1 c2 v1 v2 s1 s2;
fprintf('[5] Successfully clean observations (Three-components) \n');

%% =========================================================================
% Plot waveforms sorted by epicentral distance
% =========================================================================
% da_zh requires:
%   loca = [Latitude, Longitude]
%   epi  = [Latitude, Longitude]
da = da_zh2(loca, epicenter(1, :), 1);
% Keep only epicentral distance in km.
da = da(:, 1);
% Repeat each station distance for EW, NS, and UD.
da = repmat(da', 3, 1);
da = da(:);
% Sort all waveform component columns by epicentral distance.
[da1, ndex] = sortrows(da);
% sortedObservedData corresponds to the original variable ob1.
Sortob = ob(:, ndex);


% =================================
%              Plot
% figure('Color', 'w');
% mplot(Sortob, da1/10);
% ylabel('Epicentral distance /10');
% xlabel('Time Sample');
Fun_SaveVars('ob','Sortob','da','da1','srate','DataObs');
% Fun_PlotWave(Sortob, da1, srate, 60);

% =================================





%% =========================================================================
%                     S1. Save station coordinates
% [Longitude Latitude]
% =========================================================================
staloc = [loca(:, 2), loca(:, 1)];
stafile = fullfile(DirStaOut, '\Sta.txt');
fid = fopen(stafile, 'wt');
if fid < 0
    error('Unable to write Sta.txt.');
end
fprintf(fid, '%.4f %.4f\n', staloc');
fclose(fid);
fprintf('[S1] Successfully save stations location \n');

% =========================================================================
%                 S2. Save three-component waveform files
% Each station file has three columns:
%       [East North Up]
% Unit:
% m/s^2
% =========================================================================
numberOfValidStations = size(loca, 1);
for i = 1:numberOfValidStations
    stationName = strtrim(mm(i, :));
    outputFile = fullfile(DirWaveOut, [stationName, '.dat']);
    data = [EW(:, i), NS(:, i), UD(:, i)];
    fid = fopen(outputFile, 'wt');
    if fid < 0
        warning('Unable to write waveform file: %s', outputFile);
        continue;
    end
    % Write one EW-NS-UD triplet per line.
    fprintf(fid, '%.7f %.7f %.7f\n', data');
    fclose(fid);
end
fprintf('[S2] Successfully save wave data \n');

% =========================================================================
%                   S3. Save station information file
% =========================================================================
stationInfoFile=fullfile(DirStaOut, 'StaInfo.txt'); % 'txt' is ok or 'dat'
fid = fopen(stationInfoFile, 'wt');
if fid < 0
    error('Unable to create StaInfo');
end
Dname='DataWave';
fprintf(fid, 'Date:%s\n', Dname);
fprintf(fid, ...
    '# earthquake hypocenter: latitude[deg] longitude[deg] depth[km]\n');
fprintf(fid, '%.4f %.4f %d\n', epicenter(1, 1), epicenter(1, 2), depth(1));
fprintf(fid, '# Reference time: JST-format\n');
fprintf(fid, '%s\n', dateori);
fprintf(fid, ...
    '# Nsta Wave-Type(0/1/2=dis/vel/acc) Time-interval [s]\n');

fprintf(fid, '%d %d %.4f\n', ...
    size(loca, 1), 2, 1 / srate);

fprintf(fid, '# column of component\n');
fprintf(fid, '# east north up\n');
fprintf(fid, '1 2 3\n');

fprintf(fid, '# Strike1 Dip1 Rake1 Strike2 Dip2 Rake2\n');
fprintf(fid, '%d %d %d %d %d %d\n', ...
    macha1(1, 1), macha1(1, 2), macha1(1, 3), ...
    macha2(1, 1), macha2(1, 2), macha2(1, 3));

fprintf(fid, ...
    '# station latitude longitude Tstart Nsample\n');
numberOfSamples = size(ob, 1);
for i = 1:size(loca, 1)
    stationName = strtrim(mm(i, :));
    stationLatitude = loca(i, 1);
    stationLongitude = loca(i, 2);
    startTime = 0.0;
    fprintf(fid, '%s   %.4f   %.4f   %.0f   %d\n', ...
        stationName, ...
        stationLatitude, ...
        stationLongitude, ...
        startTime, ...
        numberOfSamples);
end
fclose(fid);
fclose('all');
fprintf('[S3] Successfully save stations Info \n');

%% =========================================================================
% Display summary
% =========================================================================
fprintf('\nProcessing completed.\n');
fprintf('Number of valid stations: %d\n', size(loca, 1));
fprintf('Number of samples: %d\n', size(ob, 1));
fprintf('Observed data matrix size: %d x %d\n', size(ob, 1), size(ob, 2));
fprintf('Sorted data matrix size: %d x %d\n', ...
    size(Sortob, 1), size(Sortob, 2));
fprintf('Sampling rate: %.2f Hz\n', srate);
fprintf('[Final] %s \n',datetime('now'));