%%%% Input arguments are date/time after 1st of month.  Date/time can be listed
%%%% as decimal Julian day (i.e. 2.00 is January 2nd 12:00 AM) or as [Jday hour
%%%% min sec msec].  Following input arguments are: window, das, ch.
%%%%
%%%% To make this script work you must be located in Segy data directory where
%%%% the dayfiles are located (e.g., R001.01 must be in the directory)

%%%% function [seis,hdr]=segycut(jtime,window,das,ch,data_stream)

function [seis,hdr]=segycut(jtime,window,das,ch,data_stream)

cur_dir = pwd;
%if str2num(cellstr(cur_dir(length(cur_dir)))) == 1; cd ..; end
%if str2num(cellstr(cur_dir(length(cur_dir)))) == 2; cd ..; end
%if str2num(cellstr(cur_dir(length(cur_dir)))) == 3; cd ..; end

if length(jtime)==2
    jtime = jtime(1) + jtime(2)/24;
elseif length(jtime)==3
    jtime = jtime(1) + jtime(2)/24 + jtime(3)/24/60;
elseif length(jtime)==4
    jtime = jtime(1) + jtime(2)/24 + jtime(3)/24/60 + jtime(4)/24/60/60;
end
if length(window) == 2; % this means time before and time after the jtime
    jtime = jtime - window(1)/24/60/60;
    window = sum(window);
end
if nargin<4; ch = 1; end
if nargin<5; data_stream = 1; end

seis = []; hdr=[];

%%%% automatically find directory where data is located
cd_directories = [
    {'G:\DATA\REVENTADOR\SEGY'} 
    {'F:\REVENTADOR\SEGY'}
    {'G:\DATA\REVENTADOR\SEGY'}
    {'F:\SG_DATA\SEGY'}
    {'\\earth\Public\lyons\Segy_data'}] % Fuego data
    %{'E:'}]
    %keyboard
for i=1:size(cd_directories,1)
    cd_directory = char(cd_directories(i));
    exist_directory = dir(cd_directory);
    if size(exist_directory,1)>0 % directory exists
        %cd(cd_directory); i=0; break;
    end
end
%if i==size(cd_directories,1);
    %error('directory not found'); return; 
%end
cur_dir = pwd;

%file_open=jhms2secs(jday,hour,min,sec,msec);
file_open = ((jtime-1)*24*60*60);
file_close = file_open + window + 1;

%%% determine direct, which is the directory where the segy files are located
directory = floor(jtime);
%data_stream = str2num(cur_dir(length(cur_dir))); 
stream = ['.0' num2str(data_stream) '/'];
direct = ['R' num2str(directory) stream];
if directory<100; direct = ['R0' num2str(directory) stream]; end
if directory<10; direct = ['R00' num2str(directory) stream]; end
if size(dir(direct),1) == 0
    %keyboard
    disp('directory not found; data not acquired on this day'); cd ..; return; 
end
cd(direct)
filenames = dir(['*.' num2str(das) '.*' num2str(ch)]);
%keyboard
%filenames = ls(['*.' num2str(das) '.*' num2str(ch)]);

%%% determine e_direct (or extra dirctory), which is the directory where the previous day's
%%% segy files are located
e_directory = floor(jtime)-1;
e_direct = ['../R' num2str(e_directory) stream];
if e_directory<100; e_direct = ['../R0' num2str(e_directory) stream]; end
if e_directory<10; e_direct = ['../R00' num2str(e_directory) stream]; end
%extras = ls([e_direct '*.' das '.*' num2str(ch)]);
extras = dir([e_direct '*.' das '.*' num2str(ch)]);
%keyboard
if size(extras,1)==0 % no files in previous directory
    first_file=1;
else
    extra_filename = [e_direct extras(size(extras,1)).name];
    first_file=0;
end
%
numfiles=size(filenames,1);
for i=first_file:numfiles
    if i==0; filename = extra_filename; % descend into the previous day's directory
    else filename = filenames(i).name; end
    [seis,hdr] = segy2mat(filename,'hdr'); % return only the header value and not data
    if isempty(hdr)==0
        fileopen = jhms2secs(hdr.day,hdr.hour,hdr.min,hdr.sec,hdr.msec);
        fileclose = fileopen + hdr.nsamp/hdr.sps;
        if file_open > fileopen & file_open < fileclose
            disp(['reading segy data from: ' filename])
            delay = file_open - fileopen;
            [seis,hdr] = segy2mat(filename,'l',delay,window);
            fileclose_tmp = fileclose;
            while file_close>fileclose_tmp & file_open<fileclose_tmp
                i=i+1;
                if i>size(filenames,1) %%% too many filenames 
                    disp(['file length exceeds files in directory: ' num2str(directory)])
                    cd(cur_dir)
                    hdr.nsamp = length(seis); 
                    hdr.channel = ch;
                    hdr.das = das;
                    hdr.stream = data_stream;
                    return
                end
                filename = filenames(i).name;
                disp(['reading segy data from: ' filename])
                [seis_tmp,hdr_tmp] = segy2mat(filename,'hdr');
                
                fileopen_tmp = jhms2secs(hdr_tmp.day,hdr_tmp.hour,hdr_tmp.min,hdr_tmp.sec,hdr_tmp.msec);
                nans = round((fileopen_tmp - fileclose_tmp)*hdr.sps); % create nans to span potential data gap
                fileclose_tmp = fileopen_tmp + hdr_tmp.nsamp/hdr.sps;
                seis_tmp(1:nans,1) = nan; 

                window_tmp = file_close - fileopen_tmp;
                delay_tmp = 0;
                
                [seis_tmp,hdr_tmp] = segy2mat(filename,'l',delay_tmp,window_tmp);
                seis = [seis; seis_tmp];
            end
            hdr.nsamp = length(seis);
            hdr.nsamp = hdr.sps*window;
            hdr.channel = ch;
            hdr.das = das;
            hdr.stream = data_stream;
            seis = seis(1:hdr.nsamp);
            cd(cur_dir)
            return
        end
    end
end