function [trace2,fs,sy,jsec]=loadsac(sacfilename)
% read in trace sac file and return vectors and power
%[trace2,fs,sy,jsec]=loadsac(sacfilename)
% inputs 
%    sacfilename='2005.295.00.36.25.0000.YZ.BOOE..BHZ.D.SAC';
%    new ~24 hour data set
% outputs
%    trace2 DC filtered trace
%    fs     sampling rate
%       not currently returned timevec timevector in s
%    sy  power vector
%    jsec seconds timestamp beginning of file
%
[sac_trace_structure] = msac_read(sacfilename);
timevec = sac_trace_structure.time; 
% have determined that is 23.393 hours
trace = sac_trace_structure.trace;  % trace
dt = timevec(2)-timevec(1);
% sampling is 50Hz again


day = sac_trace_structure.nzjday;   % timestamp of beginning
hour = sac_trace_structure.nzhour;
minute = sac_trace_structure.nzmin;
sec = sac_trace_structure.nzsec;
msec = sac_trace_structure.nzmsec;
jsec= dhmstos(day,hour,minute,sec,msec); % 
%jsec is seconds from start of year
fs = 1.0/dt; % sampling rate Hz
sprintf('day%d.wav',day);
wavefilename=ans;
Fs=20000;
clear min;
fac = max(trace)-min(trace)+1;  % 
%wavwrite(trace2/fac,Fs,wavefilename);  % make a sound file

trace = trace-median(trace); % remove DC offset
[sb,sa]=butter(4,1/(200.0*fs),'high'); %high pass at 200s
trace2 = filter(sb,sa,trace); % cut low,  filter out DC
[tb,ta]=butter(4,1.0/(5*fs),'low');  % low pass for power
sy=filter(tb,ta,abs(trace2));  % rectify and smooth
% sy looks real good.  identify all events that go above 0.1e5
% powerdb = 10*log10(sy);  % power as a function of time
%clear trace;
