classdef GWData < handle
% -- GWData --
% A GWData (Gravitational Wave Data) object stores state information
% about recent data fetching activity.
%
% Methods:
% - Methods for getting data:
% fetch - fetch data from NDS
%
% - Static Methods for GPS time conversion:
% gps_time - get GPS time now, or at some specified time
% gps_to_datenum - convert GPS time now to Matlab date number
%
% - Static Methods for Kerberos Authentication:
% make_kerberos_ready - check Kerberos ticket status (with is_kerberos_ready) and call kinit if needed
%
% How to install:
% - GWData requires the nds2-client library. Downloads and instructions
% for setting up nds2-client are found here:
% https://trac.ligo.caltech.edu/nds2
%
% - After you have nds2-client, add it to Matlab's java path:
% * Run this command in a terminal window to find out where
% nds2-client's java component was installed:
% $ nds-client-config --javaclasspath
%
% * Then create (or edit) the file "javaclasspath.txt" in your Matlab
% startup folder. Paste the nds2-client path as a line in the file.
%
% - Note: if nds2-client is found but not in the path, GWData tries to
% update the path to include it. This lets you use GWData, but it may
% have unintended effects on your Matlab environment, such as clearing
% persistent variables, breakpoints, etc. (See "doc javaaddpath" for
% more details.) Adding nds2-client to the predefined java path avoids
% these side effects.
%
% by Matthew Evans and Chris Wipf, January 2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Properties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
properties (SetAccess = public)
% maximum NDS request size
nds_max_chans = 100;
% configure prefered servers for each site (set in constructor)
site_info = struct('name', [], 'server', [], 'port', 31200);
server_info = struct('server', [], 'port', 31200);
site_list = {};
% Kerberos
kerb_path = ''; % path to Kerberos (kinit and klist)
kerb_srv = 'LIGO.ORG'; % Kerberos service principal
end
properties (SetAccess = protected)
% state information
last_start_time = [];
last_end_time = [];
last_channel_list = {};
% Kerberos
kerb_is_ready = false;
kerb_end_time = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Static Methods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods (Static)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GPS Times
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function t = gps_time(varargin)
% t = gps_time
% get GPS seconds, with sub-second precision
% (see also gps_to_datenum, tzconvert)
%
% if an argument is given, it is used instead of "now"
% (see now, datenum, datevec, etc.)
%
% % Example:
% t_now = GWData.gps_time
%
% % Example:
% t_dec3 = GWData.gps_time('3 Dec 2014 13:00:10')
%
% % Exmaple, using GWData.tzconvert:
% t0 = GWData.gps_time('06 Jan 1980 00:00:00 GMT')
%
% % Example, using full datenum format string:
% tLeapA = GWData.gps_time('Jan 01 00:00:00 GMT 2009', 'mmm dd HH:MM:SS zzz yyyy')
% tLeapB = GWData.gps_time({'Jan 01 00:00:01 GMT 2009', 'mmm dd HH:MM:SS zzz yyyy'})
%
% % Example: consistency check!
% datestr(GWData.gps_to_datenum(GWData.gps_time('3 Dec 2014 13:00:10')))
if nargin == 0
% use current time
dn = now;
else
% convert strings or other things to a date number
if nargin == 1 && isscalar(varargin{1})
dn = varargin{1};
elseif nargin == 1 && ischar(varargin{1})
dn = GWData.tzconvert(varargin{1});
elseif nargin == 1 && iscell(varargin{1})
dn = datenum(varargin{1}{:});
else
dn = datenum(varargin{:});
end
end
% convert date number to GPS time
t = GWData.datenum_to_gps(dn);
end
function t = gps_convert(userGPS, anchorGPS, isAfter, pivotGPS)
% t = gps_convert(userGPS, anchorGPS, isAfter, pivotGPS)
% convert a user specified time into a GPS value
%
% userGPS = time specified by user
% If userGPS is a string, vector, or cell array, it is
% passed to datenum for conversion. Note, datenum uses
% local time, not UTC (see datenum for more information).
%
% anchorGPS = if user time is relative, this is the absolute referece
% (the default is the current GPS)
%
% pivotGPS = numerical GPS values less than this are considered to
% be relative, larger values (corresponding to later times)
% are considered absolute.
% The default is January 1, 2000 (GPS 630730000).
%
% isAfter = relative time is after anchor time? If true
% tGPS = anchorGPS + userGPS. Else relative time
% is before anchor (e.g., tGPS = anchorGPS - userGPS).
% (default is false)
%
% t = userGPS converted to GPS second
%
% Example: GPS second for 3 Dec 2014 13:00:10
% t = gps_convert('3 Dec 2014 13:00:10');
%
% Example: 100 seconds BEFORE the current time.
% t = gps_convert(100);
%
% Example: 100 seconds AFTER 3 Dec 2014 13:00:10.
% t = gps_convert(100, '3 Dec 2014 13:00:10', true);
%
% Example: 100 seconds BEFORE 3 Dec 2014 13:00:10 GMT
% t = gps_convert(100, {'2014 Dec 03 13:00:10 GMT', 'yyyy mmm dd HH:MM:SS zzz'});
% deal with arguments
if nargin < 4
pivotGPS = 630730000;
end
if nargin < 3
isAfter = false;
end
if nargin < 2
anchorGPS = GWData.gps_time();
end
% convert all GPS times to scalars
if ~isscalar(userGPS)
userGPS = GWData.gps_time(userGPS);
end
if ~isscalar(anchorGPS)
anchorGPS = GWData.gps_time(anchorGPS);
end
if ~isscalar(pivotGPS)
pivotGPS = GWData.gps_time(pivotGPS);
end
% check for relativity
if userGPS >= pivotGPS
t = userGPS;
else
if isAfter
t = anchorGPS + userGPS;
else
t = anchorGPS - userGPS;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time Standards
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function t = gps_to_datenum(t_gps)
% t = gps_to_datenum(gpsTime)
% Converts GPS second to a local matlab date number
% (see also gps_time, gps_to_unix, unix_to_datenum)
t = GWData.unix_to_datenum(GWData.gps_to_unix(t_gps));
end
function t = datenum_to_gps(dn)
% t = gps_to_datenum(gpsTime)
% Converts local matlab date number to GPS second
% (see also gps_time, datenum_to_unix, unix_to_gps)
t = GWData.unix_to_gps(GWData.datenum_to_unix(dn));
end
function t = datenum_to_unix(dn)
% t = datenum_to_unix(gpsTime)
% Converts local matlab date number to UNIX time
% (see also gps_time, datenum_to_unix, unix_to_gps)
% use matlab date number (note that milliseconds are lost
% when making javaDate, and added back in later)
[y, m, d, h, mn, s] = datevec(dn);
javaDate = java.util.Date(y - 1900, m - 1, d, h, mn, floor(s));
% UNIX time, noting that java time is in milliseconds
% and that we lost some milliseconds when creating the javaDate
t = javaDate.getTime / 1000 + (s - floor(s));
end
function dn = unix_to_datenum(t_unix)
% t = unix_to_datenum(gpsTime)
% Converts UNIX time to local matlab date number
% (see also gps_time, datenum_to_unix, unix_to_gps)
% This is trickier than it looks!
% Note that the datenum is not continuous across daylight savings
% changes, while UNIX time is, so the following does not work:
% localOffset = datenum('Jan 01 1970 GMT', 'mmm dd yyyy zzz');
% t = (t_unix / 86400) + localOffset;
% go back to datenum through java, reversing datenum_to_unix
javaDate = java.util.Date(t_unix * 1000);
cal = java.util.GregorianCalendar;
cal.setTime(javaDate);
y = cal.get(java.util.Calendar.YEAR);
mo = cal.get(java.util.Calendar.MONTH);
d = cal.get(java.util.Calendar.DAY_OF_MONTH);
h = cal.get(java.util.Calendar.HOUR_OF_DAY);
mi = cal.get(java.util.Calendar.MINUTE);
s = cal.get(java.util.Calendar.SECOND);
ms = cal.get(java.util.Calendar.MILLISECOND);
dn = datenum(y, mo + 1, d, h, mi, s + ms / 1000);
% cross check...
%[y , mo, d, h, mi, s]
%[y , mo, d, h, mi, s] = datevec(dn)
end
function t = gps_to_unix(t_gps)
% t = gps_to_unix(gpsTime)
% Converts GPS second to UNIX time
% (see also gps_time, unix_to_gps, unix_to_datenum)
% this is the offset between Unix time (started Jan 1, 1970)
% and GPS time (started Jan 6, 1980)
gpsOffset = 315964800;
% correct for leap seconds
t = t_gps + gpsOffset - GWData.leap_seconds_in_gps(t_gps);
end
function t = unix_to_gps(t_unix)
% t = unix_to_gps(gpsTime)
% Converts UNIX time to GPS second
% (see also gps_time, gps_to_unix, gps_to_datenum)
% this is the offset between Unix time (started Jan 1, 1970)
% and GPS time (started Jan 6, 1980)
gpsOffset = 315964800;
% convert to GPS time and add leap seconds
t = t_unix - gpsOffset + GWData.leap_seconds_in_unix(t_unix);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Leap Seconds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function leap_num = leap_seconds_in_unix(t_unix)
% leap second correction
% datenum is tied to UNIX time, where leap seconds are removed
% so we have to put them back to convert from GPS to UNIX time
% leap second UNIX times
[~, leap_unix] = GWData.leap_seconds();
% count accumulated leap seconds
leap_num = sum(leap_unix < t_unix);
end
function leap_num = leap_seconds_in_gps(t_gps)
% leap second correction
% datenum is tied to UNIX time, where leap seconds are removed
% so we have to put them back to convert from GPS to Unix time
% leap second GPS times
[leap_gps, ~] = GWData.leap_seconds();
% count accumulated leap seconds
leap_num = sum(leap_gps < t_gps);
end
function [leap_gps, leap_unix, leap_datenum] = leap_seconds()
% get timestamps for when leap seconds happened
% note that leap_gps and leap_unix are timezone independent
% while leap_datenum depends on timezone
% NOTE: execution of this function takes tens of milliseconds!
% tic; [leap_gps, leap_unix, leap_datenum] = GWData.leap_seconds(); toc
% Elapsed time is 0.070834 seconds.
% leap seconds since 1980 (updated January, 2015)
leap_date = [...
'Jul 01 1981 GMT'
'Jul 01 1982 GMT'
'Jul 01 1983 GMT'
'Jul 01 1985 GMT'
'Jan 01 1988 GMT'
'Jan 01 1990 GMT'
'Jan 01 1991 GMT'
'Jul 01 1992 GMT'
'Jul 01 1993 GMT'
'Jul 01 1994 GMT'
'Jan 01 1996 GMT'
'Jul 01 1997 GMT'
'Jan 01 1999 GMT'
'Jan 01 2006 GMT'
'Jan 01 2009 GMT'
'Jul 01 2012 GMT'
'Jul 01 2015 GMT'
'Jul 01 2018 GMT' % estimate
'Jul 01 2021 GMT' % estimate
'Jul 01 2024 GMT' % estimate
'Jul 01 2027 GMT' % estimate
'Jul 01 2030 GMT' % estimate
'Jul 01 2033 GMT' % estimate
'Jul 01 2036 GMT' % estimate
'Jul 01 2039 GMT']; % estimate (~1ms / day => ~0.33s / year)
% this is the offset between Unix time (started Jan 1, 1970)
% and GPS time (started Jan 6, 1980)
gpsOffset = 315964800;
% convert to unix/gps/date numbers
% datenum seems to have DST issues, so use Java functions instead
% leap_datenum = datenum(leap_date, 'mmm dd yyyy zzz');
leap_unix = zeros(size(leap_date, 1), 1);
leap_gps = zeros(size(leap_date, 1), 1);
leap_datenum = zeros(size(leap_date, 1), 1);
dateFormatter = java.text.SimpleDateFormat('MMM dd yyyy z');
for n = 1:size(leap_date, 1)
% get UNIX time for each (comes in milliseconds)
leap_unix(n) = dateFormatter.parse(leap_date(n, :)).getTime/1000;
% remove gps offset, and add accumulated leap seconds
% (can't use unix_to_gps, since that comes here!)
leap_gps(n) = (leap_unix(n) - gpsOffset) + n;
leap_datenum(n) = GWData.unix_to_datenum(leap_unix(n));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time Zone Paring
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function date_num = tzconvert(date_str)
% TZCONVERT Convert date string to serial date number (time zone aware)
% Matlab's builtin DATENUM function ignores time zone specifiers when
% converting date strings (unless an explicit date string format is
% provided). For example:
% >> datenum('01-jan-2015 00:00:00 CST')-datenum('01-jan-2015 00:00:00 PST')
%
% ans =
%
% 0
%
% TZCONVERT works around this shortcoming:
% >> tzconvert('01-jan-2015 00:00:00 CST')-tzconvert('01-jan-2015 00:00:00 PST')
%
% ans =
%
% -0.0833
%
% TZ formats accepted include: time offsets (+0100, -0800, etc),
% US time zone abbreviations (EST, PDT, etc), and UTC or GMT.
%
% See also: DATENUM, DATESTR
%%% Look for time zone specifiers in date_str
tz_re = '(\<(+|-)[0-2][0-9][0-5][0-9]|UTC|GMT|(A|E|C|M|P|AK|HA)(S|D)T\>)';
tz_idx = regexpi(date_str, tz_re, 'tokenExtents');
if numel(tz_idx) > 1
error(['multiple time zone specifiers found in date string: ' date_str]);
end
%%% Convert to datenum
% no time zones => nothing to do
if isempty(tz_idx)
date_num = datenum(date_str);
return;
end
% split the time zone from the rest of the string
tz_idx = tz_idx{1};
tz = date_str(tz_idx(1):tz_idx(2));
date_str = [date_str(1:tz_idx(1)-1) date_str(tz_idx(2)+1:end)];
% rewrite the date string in a known format
date_str = [datestr(date_str, 'dd-mmm-yyyy HH:MM:SS') ' ' tz];
% use java date libraries to make the conversion
% datenum should be able to do this, but it's broken by DST
% for example, the following is not zero (replace PDT with your local
% time zone):
% datenum('24-mar-2015 02:17:00')-datenum('24-mar-2015 02:17:00 PDT', 'dd-mmm-yyyy HH:MM:SS zzz')
dateFormatter = java.text.SimpleDateFormat('dd-MMM-yyyy HH:mm:ss z');
date_num = GWData.unix_to_datenum(dateFormatter.parse(date_str).getTime/1000);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kerberos Authentication
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function env_str = kerberos_env(kerb_path)
% env_str = kerberos_env(kerb_path)
% set LANG and PATH environment variables needed to invoke Kerberos
% commands.
% locale settings influence the printing of dates by klist
env_str = 'LANG=posix LC_TIME=posix';
if ispc
env_str = 'set LANG=posix && set LC_TIME=posix &&';
end
% special path setting for Mac OS X (forces the use of MacPorts)
if ismac && (nargin < 1 || isempty(kerb_path))
kerb_path = '/opt/local/bin';
end
% special path setting for Windows (forces the use of MIT Kerberos)
if ispc && (nargin < 1 || isempty(kerb_path))
kerb_path = 'C:\Program Files\MIT\Kerberos\bin';
end
% apply path setting
path_now = getenv('PATH');
if ~strncmp(path_now, kerb_path, numel(kerb_path))
if ~ispc
env_str = [env_str ' PATH=' kerb_path ':' path_now];
else
env_str = [env_str ' path ' kerb_path ';' path_now ' &&'];
end
end
% check for kerberos (system returns 0 when ok)
if ~ispc
[isNotOk, str] = system([env_str ' which kinit']);
else
[isNotOk, str] = system([env_str ' where kinit']);
end
if isNotOk
error('kinit not found. Unable to setup Kerberos ticket.')
elseif ~strncmp(str, kerb_path, numel(kerb_path))
warning('using kinit at %s, expect it at %s', str, kerb_path)
end
end
function [is_ready, end_time] = make_kerberos_ready(srv, kerb_path)
% [is_ready, end_time] = make_kerberos_ready(srv, kerb_path)
% look for a valid kerberos ticket for a particular service principal
% and use kinit to make a new ticket if none is found.
%
% srv = Kerberos service principal (default = 'LIGO.ORG')
% kerb_path = path for kinit and klist (default = '/opt/local/bin' for Macs)
%
% is_ready = true if a ticket is found or created, false otherwise
% end_time = matlab date number for end of ticket validity (see datenum)
% default service principal
if nargin < 1
srv = 'LIGO.ORG';
end
if nargin < 2
kerb_path = '';
end
env_str = GWData.kerberos_env(kerb_path);
% initialize kerberos
[is_ready, end_time] = GWData.is_kerberos_ready(srv);
if is_ready
% we already have a ticket
return
else
% we need to get a ticket
fprintf('== GW Data: Kerberos authentication ==\n')
if ~ispc
username = input([srv ' user name: '], 's');
if system([env_str ' kinit ' username '@' srv]);
error('kinit failed. Bad password?')
end
else
if system([env_str ' "MIT Kerberos.exe" -kinit']);
error('kinit failed. Bad password?')
end
end
% now get ticket info
[is_ready, end_time] = GWData.is_kerberos_ready(srv);
end
end
function [is_ready, end_time] = is_kerberos_ready(srv, kerb_path)
% [is_ready, end_time] = is_kerberos_ready(srv, kerb_path)
% look for a valid kerberos key for a particular service principal
% this assumes that the path is set properly (see
% make_kerberos_ready)
%
% %%% no cache?
% klist: No credentials cache found (ticket cache FILE:/tmp/krb5cc_501)
%
% %%% found cache...
% Ticket cache: FILE:/tmp/krb5cc_501
% Default principal: matthew.evans@LIGO.ORG
%
% Valid starting Expires Service principal
% 12/09/2014 15:24:34 12/10/2014 15:24:34 krbtgt/LIGO.ORG@LIGO.ORG
if nargin < 2
kerb_path = '';
end
env_str = GWData.kerberos_env(kerb_path);
% first line of ticket info
first_line = 4;
% return values, unless we find otherwise
is_ready = false;
end_time = 0;
% call klist to get list of tickets
[exit_status, result] = system([env_str ' klist']);
% klist failed, so we're definitely not ready
if exit_status ~= 0
return
end
% split result into lines
%result_lines = strsplit(result, {'\n', '\r'});
result_lines = GWData.splitlines(result);
if ~iscell(result_lines) || numel(result_lines) < first_line
% not enough lines
return
end
% look for service principal in tickets
hasSP = strfind(result_lines, srv);
for n = first_line:numel(hasSP)
% if this ticket mentions the SP, look at the date strings
if ~isempty(hasSP{n})
% looking for lines like:
% 12/09/2014 15:24:34 12/10/2014 15:24:34 krbtgt/LIGO.ORG@LIGO.ORG
date_time_strs = GWData.splitwords(result_lines{n});
if numel(date_time_strs) > 4
% get data number for end time, and compare to time now
end_time = datenum([date_time_strs{3} ' ' date_time_strs{4}]);
if end_time > now
% this ticket is ok!
is_ready = true;
return
end
end
end
end
% didn't find a valid ticket
is_ready = false;
end
% same as strsplit(result, {'\n', '\r'});
% for old matlab versions
function strlist = splitlines(str)
nn = strfind(str, char(10)); % \n
nr = strfind(str, char(13)); % \r
nnr = strfind(str, [char(10) char(13)]);
nrn = strfind(str, [char(13) char(10)]);
nAll = unique([nn, nr, nnr, nrn]);
if isempty(nAll)
strlist = str;
else
nLast = 0;
strlist = {};
for n = 1:numel(nAll)
nNext = nAll(n);
if nNext > nLast + 1
strlist(end + 1) = {str((nLast + 1):(nNext - 1))}; %#ok
end
nLast = nNext;
end
end
end
% same as strsplit(result, {' ', '\t'});
% for old matlab versions
function strlist = splitwords(str)
nt = strfind(str, char(9)); % \t
ns = strfind(str, ' ');
nAll = unique([nt, ns]);
if isempty(nAll)
strlist = str;
else
nLast = 0;
strlist = {};
for n = 1:numel(nAll)
nNext = nAll(n);
if nNext > nLast + 1
strlist(end + 1) = {str((nLast + 1):(nNext - 1))}; %#ok
end
nLast = nNext;
end
if nLast < numel(str)
strlist(end + 1) = {str((nLast + 1):end)};
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data Processing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = resample(x, p, q)
% y = my_resample(x, p, q)
% use fitting to resample without end glitches
% if x is too short, just demean and return y
if numel(x) < 10
xm = mean(x);
y = resample(x - xm, p, q) + xm;
return
end
% fit a polynomial to the end points
toff = numel(x) / 2;
tx = (0:(numel(x) - 1))' - toff;
xm = abs(mean(x(1:3) - x(end-2:end))); % shift from start to end
xs = std(x(1:3)) + std(x(end-2:end)); % spread of data
if numel(x) < 30 || xm < 5 * xs
% not may points, or pretty flat data, so use a linear fit
pf = polyfit(tx([1:3, (end-2:end)]), x([1:3, (end-2:end)]), 1);
else
% more points, so use a cubic fit
pf = polyfit(tx([1:10, (end-9:end)]), x([1:10, (end-9:end)]), 3);
end
% make a new x that is zero at both ends
xz = x - polyval(pf, tx);
% resample
yz = resample(xz, p, q);
% add back in polynomial fit
ty = (0:(numel(yz) - 1))' * q / p - toff;
y = yz + polyval(pf, ty);
%plot(tx, [x, xz], ty, [y, yz, polyval(pf, ty)])
%pause
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global Variable Save/Restore
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% These functions need to import global variables into their own
% workspace. So they stash their own local variables in Matlab's
% appdata, to avoid colliding with the global variable names.
function globals = save_globals()
globals = struct();
vars = who('global');
for n = 1:numel(vars)
locals = struct();
locals.n = n;
locals.vars = vars;
locals.globals = globals;
setappdata(0, 'GWData_sg_locals', locals);
setappdata(0, 'GWData_sg_var', vars{n});
eval(['global ' vars{n}]);
setappdata(0, 'GWData_sg_val', eval(getappdata(0, 'GWData_sg_var')));
locals = getappdata(0, 'GWData_sg_locals');
n = locals.n; %#ok
vars = locals.vars;
globals = locals.globals;
globals.(vars{n}) = getappdata(0, 'GWData_sg_val');
end
if isappdata(0, 'GWData_sg_locals')
rmappdata(0, 'GWData_sg_locals');
rmappdata(0, 'GWData_sg_var');
rmappdata(0, 'GWData_sg_val');
end
end
function restore_globals(globals)
vars = fieldnames(globals);
for n = 1:numel(vars)
locals = struct();
locals.n = n;
locals.vars = vars;
locals.globals = globals;
setappdata(0, 'GWData_rg_locals', locals);
setappdata(0, 'GWData_rg_var', vars{n});
setappdata(0, 'GWData_rg_val', globals.(vars{n}));
eval(['global ' vars{n}]);
eval([getappdata(0, 'GWData_rg_var') '= getappdata(0, ''GWData_rg_val'');']);
locals = getappdata(0, 'GWData_rg_locals');
n = locals.n; %#ok
vars = locals.vars;
globals = locals.globals;
end
if isappdata(0, 'GWData_rg_locals')
rmappdata(0, 'GWData_rg_locals');
rmappdata(0, 'GWData_rg_var');
rmappdata(0, 'GWData_rg_val');
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Public Methods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods
function obj = GWData()
% setup kerberos path
if ismac
obj.kerb_path = '/opt/local/bin';
end
% setup site info
obj.site_info(1:3) = obj.site_info;
obj.site_info(1).name = 'L1';
obj.site_info(1).server = 'nds.ligo-la.caltech.edu';
obj.site_info(2).name = 'H1';
obj.site_info(2).server = 'nds2.ligo-wa.caltech.edu';
obj.site_info(3).name = 'C1';
obj.site_info(3).server = 'nds40.ligo.caltech.edu';
% configure backup servers
obj.server_info(1:4) = obj.server_info;
obj.server_info(1).server = 'nds.ligo.caltech.edu';
obj.server_info(2).server = 'ldas-pcdev1.ligo-la.caltech.edu';
obj.server_info(3).server = 'ldas-pcdev1.ligo-wa.caltech.edu';
obj.server_info(4).server = 'nds2.ligo.caltech.edu';
% save list of site names
obj.site_list = {obj.site_info.name};
% check for local environment
ndsserver = getenv('NDSSERVER');
if ~isempty(ndsserver)
ndsserver = regexp(ndsserver, ',', 'split');
ndsserver = regexp(ndsserver{1}, ':', 'split');
server = ndsserver{1};
port = 31200;
if numel(ndsserver) > 1
port = str2double(ndsserver{2});
end
% need to know which site the server is associated with
ifo = getenv('IFO');
if ~isempty(ifo)
site_idx = strcmp(ifo, obj.site_list);
if any(site_idx)
% replace existing site_info entry
obj.site_info(site_idx).server = server;
obj.site_info(site_idx).port = port;
else
% add new site_info entry and update site_list
obj.site_info(end+1).name = ifo;
obj.site_info(end).server = server;
obj.site_info(end).port = port;
obj.site_list = {obj.site_info.name};
end
disp(['Using NDS server ' ndsserver{1} ' for ' ifo]);
else
warning('Ignoring local $NDSSERVER setting for unknown site');
end
end
end
function [data, t, info] = ...
fetch(obj, start_time, end_time, channel_list, data_rate)
% [data, t, info] = gwd.fetch(start_time, end_time, channel_list, data_rate)
%
% Fetch data from NDS2 servers.
%
% == Arguments ==
% start_time = data starts at this GPS second
% If start_time corresponds to a date before January 1, 2000
% (e.g., start_time < 630748814), it is taken as the number
% of seconds BEFORE the current time. To get data from 1 hour
% ago, for instance, use start_time = 3600.
% If start_time is a string, vector, or cell array, it is
% passed to GWData.gps_time for conversion. Note, this
% uses local time, not UTC (see gps_time and tzconvert).
%
% end_time = data ends at this GPS second
% If end_time corresponds to a date before January 1, 2000,
% it is taken as the duration of the data segment.
% If end_time is a string, vector, or cell array, it is
% passed to GWData.gps_time for conversion.
%
% channel_list = channel name, or cell array of channel names to fetch
% In many cases, adding ".mean,s-trend" (or ".mean,m-trend")
% to the end of a channel name is sufficient to get the
% second trend (or minute trend) mean value.
%
% data_rate = rate to resample the data to before returning
% If not specified, or given as [], the rate of the first
% channel in the channel list will be used. If given as
% 'mDV', data will be returned in a structure array with
% one entry for each channel.
%
% == Return Values ==
% data = Nsample x Nchannel array of data
% (i.e., data(:, 2) is the data for the second channel)
% or in 'mDV' format, data = 1xNchannel struct array with fields:
% name = channel name
% data = full data points
% rate = sample rate
% start = actual start GPS
% duration = actual data duration
% t = time vector relative to start_time
% info = structure with the following fields
% start_time = actual start GPS
% end_time = actual end GPS
% duration = actual data duration
% data_buffers = full data returned from NDS fetch function
% conn = NDS connection used to fetch this data
% The connection will be closed, unless fetch is invoked
% with a duration of 0. This may be useful for later calls
% to findChannels, etc.
%
% ====== Examples ======
%
% % start with a GW Data object
% gwd = GWData;
%
% % get the minute trend of IMC input power from a few hours ago
% [data, t] = gwd.fetch(5 * 3600, 1000, 'L1:IMC-PWR_IN_OUT16.mean,m-trend');
% plot(t, data)
%
% % get 16Hz data for Y-arm power and input power
% % start time is 2 hours ago, duration is 1000s
% [data, t] = gwd.fetch(2 * 3600, 1000, ...
% {'L1:ASC-Y_TR_B_SUM_OUT16', 'L1:IMC-PWR_IN_OUT16'});
% plot(t, [data(:, 1) / 1e3, data(:, 2)])
%
% % -- find a channel, then get data --
% % initialize Kerberos (if not already done)
% gwd.make_kerberos_ready;
%
% % connect to server and find some channels
% conn = nds2.connection('nds.ligo-la.caltech.edu', 31200);
% fc = conn.findChannels('L1:LSC-X_TR*OUT*.mean*')
%
% % look through the channels to find the one you want... say #3
% %
% % and then get some data for it
% [data, t, info] = gwd.fetch(5 * 3600, 1000, fc(3).getName);
% t0 = info.start_time; % start time GPS
%
% % change server for H1 (note: H1 is second in GWData.site_list)
% gwd.site_info(2).server = 'nds.ligo.caltech.edu';
%
% ====== Troubleshooting ======
% Start by looking under "Matlab Tools" on the remote data access wiki:
% https://wiki.ligo.org/viewauth/RemoteAccess
%
% --> NDS2 not found <--
% If you installed NDS2, but matlab can't find it try
% >> edit classpath.txt
%
% and add this to the end of the file:
% # NDS2
% /opt/local/lib/java
% (or replace /opt/local/lib/java with your java class path)
%
% --> kinit not found <--
% kinit is used to obtain and cache Kerberos ticket-granting ticket.
% NDS2 needs this for authentication. Try 'which kinit' on the
% command line and Google from there.
%
% --> kinit failed <--
% wrong user name or password? This should be your "albert.einstein"
% username and corresponding ligo.org password.
%
% TODO:
% catch NDS errors and try backup servers
% use CIS to diagnose "No data found" problems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setup dependencies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check for NDS2
if ~exist('nds2.connection', 'class')
% check hardcoded MacPorts path /opt/local/bin, which may not otherwise
% be picked up by the matlab app's PATH environment variable
pathsToCheck = {'', '/opt/local/bin/'};
for n = 1:numel(pathsToCheck)
[status, output] = system([pathsToCheck{n} 'nds-client-config --javaclasspath']);
if ~status
path = deblank(output);
disp(['Found nds2-client library in ' path]);
% javaaddpath is evil: it clears global variables and many
% other aspects of the environment. So try to save and restore
% the globals, and ask the user to install nds2-client in the
% static java path.
warning(['nds2-client library was found in ' path ', which is not in Matlab''s java path.' char(10) ...
'Please see "help GWData" for instructions to complete the installation.']);
globals = GWData.save_globals();
javaaddpath(path);
GWData.restore_globals(globals);
break;
end
end
if ~exist('nds2.connection', 'class')
error('Can''t find nds2-client (use "help GWData.fetch" for more info)');
end
end
% if this was just for initialization...
if nargin == 1
% initialize kerberos
[obj.kerb_is_ready, obj.kerb_end_time] = ...
GWData.make_kerberos_ready(obj.kerb_srv, obj.kerb_path);
% and return empty
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% interpret arguments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% record present time for later use
gps_time_now = floor(GWData.gps_time);
% validate arguments
if ~isa(start_time, 'numeric') || ~isscalar(start_time)
if ~ischar(start_time)
error('start_time must be a GPS time or date string');
elseif any(strncmp(start_time, obj.site_list, 2))
error('start_time looks like a channel name');
end
end
if ~isa(end_time, 'numeric') || ~isscalar(end_time)
if ~ischar(end_time)
error('end_time must be a GPS time or date string');
elseif any(strncmp(end_time, obj.site_list, 2))
error('end_time looks like a channel name');
end
end
if ischar(channel_list)
channel_list = {channel_list};
end
if ~iscellstr(channel_list)
error('channel list must be a cell array of channel names');
end
% start and end times to gps values
start_time = GWData.gps_convert(start_time, gps_time_now, false);
end_time = GWData.gps_convert(end_time, start_time, true);
% make start and end times into integers
% (and ensure that floating point errors don't cause problems)
start_time = floor(start_time + 1e-6);
end_time = ceil(end_time - 1e-6);
% check data rate (default to first channel rate if not given)
if nargin < 5
data_rate = [];
end
% set some background info
num_channel = numel(channel_list);
duration = end_time - start_time;
% check validity of end_time
if duration < 0
error('end_time is before start_time (%d < %d)', end_time, start_time)
end
latency_time = gps_time_now - end_time;
if latency_time < 0
error('end_time is in the future (%d > %d)', end_time, gps_time_now)
end
if latency_time < 30
error('end_time is too close to the present. Try again in 30s.')
end
% if no channels to read, return empty
if num_channel == 0
data = [];
t = [];
info.start_time = start_time;
info.end_time = end_time;
info.duration = duration;
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% look for minute trends, and adjust times
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nMinTrend = strfind(channel_list, 'm-trend');
hasMinTrend = ~isempty([nMinTrend{:}]);
% if there are minute trends in the list, adjust to even minutes
if hasMinTrend
if rem(start_time, 60) ~= 0
start_time = start_time - rem(start_time, 60);
warning('start_time adjusted to even minute for m-trend, now %d', start_time);
end
if rem(end_time, 60) ~= 0
end_time = end_time + 60 - rem(end_time, 60);
warning('end_time adjusted to even minute for m-trend, now %d', end_time);
end
% update the duration
duration = end_time - start_time;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% connect to network data server
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine site
site_name = channel_list{1}(1:2);
site_num = find(strncmp(obj.site_list, site_name, 1), 1);
if isempty(site_num)
error('Unknown site name %s in channel %s.', site_name, channel_list{1});
end
server = obj.site_info(site_num).server;
port = obj.site_info(site_num).port;
% check others
for n = 2:num_channel
if ~strcmp(site_name, channel_list{2}(1:2))
error('All channels must be from the same site.');
end
end
% initialize kerberos, unless connecting to the NDS1 port 8088
if port ~= 8088
[obj.kerb_is_ready, obj.kerb_end_time] = ...
GWData.make_kerberos_ready(obj.kerb_srv, obj.kerb_path);
end
% connect to server
disp(['Connecting to NDS server ' server]);
conn = nds2.connection(server, port);
% if no data to read, return empty
if duration == 0
data = [];
t = [];
info.start_time = start_time;
info.end_time = end_time;
info.duration = duration;
info.conn = conn;
return
end
s = '';
if numel(channel_list) > 1
s = 's';
end
disp(['Fetching ' num2str(numel(channel_list)) ' channel' s ...
', start GPS ' num2str(start_time), ', duration ' num2str(duration) ' sec']);
cleanup_conn = onCleanup(@() conn.close());
% Provide a status bar when many channels are requested
h_win = 0;
if numel(channel_list) > obj.nds_max_chans
h_win = waitbar(0, ['Fetching ' num2str(numel(channel_list)) ' channels...'], ...
'CreateCancelBtn', 'setappdata(gcbf, ''canceling'', 1)', 'Name', 'GWData');
setappdata(h_win, 'myX', 0);
setappdata(h_win, 'canceling', 0);
cleanup_window = onCleanup(@() delete(h_win));
end
% Define a callback to update the status bar
function cb(inc)
if h_win ~= 0
if getappdata(h_win, 'canceling')
error('GWData:userCancelled', 'NDS data request cancelled')
end
x = getappdata(h_win, 'myX') + inc/numel(channel_list);
setappdata(h_win, 'myX', x);
waitbar(x, h_win);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get data for each channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data_buffers = obj.rfetch(conn, channel_list, start_time, end_time, @cb);
clear cleanup_conn;
if h_win ~= 0
clear cleanup_window;
drawnow;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% resample and repackage data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% validate
if isempty(data_buffers)
error('Unable to fetch data from %s.', server);
end
% determine rate, if not already known
if isempty(data_rate)
data_rate = data_buffers(1).getLength / duration;
end
% output format selection
if strcmp(data_rate, 'mDV')
t = [];
data = struct('name', {}, 'data', {}, 'rate', {}, 'start', {}, 'duration', {});
for n = 1:numel(data_buffers)
data(n) = struct('name', data_buffers(n).getChannel.getName, ...
'data', data_buffers(n).getData, ...
'rate', data_buffers(n).getChannel.getSampleRate, ...
'start', data_buffers(n).getGpsSeconds + data_buffers(n).getGpsNanoseconds/10^9, ...
'duration', data_buffers(n).getLength/data_buffers(n).getChannel.getSampleRate);
end
else
% initialize time vector and data space
num_sample = duration * data_rate;
t = (0:(num_sample - 1))' / data_rate;
data = zeros(num_sample, num_channel);
% loop over channels to resample, if needed
for n = 1:num_channel
this_data = double(data_buffers(n).getData);
if data_buffers(n).getLength == num_sample
% no need to resample
data(:, n) = this_data;
else
% resample the data
data(:, n) = GWData.resample(this_data, ...
num_sample, data_buffers(n).getLength);
end
end
end
% for the record
info.start_time = start_time;
info.end_time = end_time;
info.duration = duration;
info.conn = conn;
info.data_buffers = data_buffers;
obj.last_start_time = start_time;
obj.last_end_time = end_time;
obj.last_channel_list = channel_list;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Protected Methods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods (Access = protected)
function bufs = rfetch(obj, conn, chan_list, start_time, end_time, cb)
% recursive fetch helper function
bufs = {};
if numel(chan_list) < 1
return;
end
if numel(chan_list) > obj.nds_max_chans
for n = 1:obj.nds_max_chans:numel(chan_list)
endN = min(numel(chan_list), n+obj.nds_max_chans-1);
bufsN = obj.rfetch(conn, chan_list(n:endN), start_time, end_time, cb);
bufs = [bufs bufsN]; %#ok
end
return;
end
try
bufs = conn.fetch(start_time, end_time, chan_list);
cb(numel(chan_list));
return;
catch exc
if strcmp(exc.identifier, 'GWData:userCancelled')
rethrow(exc);
end
if numel(chan_list) == 1
disp(['Failed to fetch channel ' chan_list{1}]);
cb(numel(chan_list));
rethrow(exc);
end
chan_list1 = chan_list(1:floor(end/2));
chan_list2 = chan_list(floor(end/2)+1:end);
bufs = [obj.rfetch(conn, chan_list1, start_time, end_time, cb)...
obj.rfetch(conn, chan_list2, start_time, end_time, cb)];
end
end
end
end