maci/ 0000755 0000765 0000765 00000000000 11165024663 010117 5 ustar tah tah maci/tsview/ 0000755 0000765 0000765 00000000000 11165417772 011447 5 ustar tah tah maci/tsview/AvgRes.m 0000644 0000765 0000765 00000004415 11156217426 013012 0 ustar tah tah function AvgRes(ha, AvgLength,AvgTime) % Function to plot averaged Residuals. Plotting is done by replacing the X % and Y data on the plots themselves. % MOD TAH 090312: Added boundary time AvgTime (zero if not set) % % Find the residuals associated with this set of axes axes(ha); %he = findobj(ha,'Tag','PlotErrb'); hr = findobj(ha,'Tag','PlotRes'); % % Remove error bars from plot % delete(he); % % Now get the data rall = get(hr(1),'UserData'); times = rall{1}; res = rall{2}; err = rall{3}; delete(hr); % % Now average the data start = min(times); stop = max(times); % dyr = AvgLength/365.25; % MOD TAH 090312: See if need to adjust for AvgTime (boundary should fall % on requested day). if AvgTime ~= 0 starto = start; start = floor((start-AvgTime)/dyr)*dyr + AvgTime; fprintf('Average start: Data %f , boundary %f Step %f yr\n',starto,start,dyr); end % avT = []; avR = []; avS = []; for t = start:dyr:stop sel = find(times >= t & times < t+dyr); ats = times(sel); ars = res(sel); ass = err(sel); % Now find the weighted means num = length(ats); if( num > 1 ) w = 1./ass.^2; Wmean = w*ars'/sum(w); r = ars - Wmean; nrms = sqrt(sum(r.^2.*w)/(num-1)); werr = sqrt(1/sum(w)); Tmean = w*ats'/sum(w); avT = [avT Tmean]; avR = [avR Wmean]; avS = [avS werr]; end end % % OK plot the averaged values he = errorbar(avT,avR,avS,'ks-') ; set(he,'MarkerSize',9,'MarkerFaceColor',[1 0 0],'LineWidth',2); % Save the data so that we can average later if requested set(he,'Tag','PlotRes'); % set(he(1),'Tag','PlotErrb'); set(he,'UserData',{times,res,err}); ymn = min(avR); ymx = max(avR); dy = ymx-ymn; set(ha,'Ylim',[ymn-0.2*dy ymx+0.2*dy]); % % Now compute the WRMS and NRMS num = length(avT); if( num > 0 ) w = 1./avS.^2; nrms = sqrt(sum(avR.^2.*w)/(num-1)); wrms = nrms*sqrt(num/sum(w)); txt = sprintf('Av %6.1f days, # %d, WRMS %6.2f mm NRMS %6.1f', ... AvgLength, num, wrms, nrms); % See if Average text already present. If it is delete it ht = findobj(ha,'Tag','AvgText'); if ~isempty(ht) delete(ht); end ht = text(0.01,1.050,txt, ... 'VerticalAlignment','bottom','Units','normalized','FontWeight','bold', ... 'Tag','AvgText'); end maci/tsview/CalcRealSigma.m 0000644 0000765 0000765 00000004733 11155777336 014267 0 ustar tah tah function nrmsR = CalcRealSigma(Data,res,PlotTitle) % Function to compute realistic sigmas % Algorithm here uses the changes in the nrms with increasing averaging % interval to predict infinite averaging NRMS % % First get some bounds on how much averaging we can do. times = Data(1,:); err = Data(3,:); start = min(times); stop = max(times); minav = 7; maxav = (stop-start)*365.25/10; if maxav < minav; minav = maxav/4; end numav = fix(maxav/minav); % % zero the arrays we need savsta = []; timsta = []; % Get the "white-noise" part of the error budget by differincing adjucent % data points. nd = length(res); dres = res(1:nd-1)-res(2:nd); dvar = err(1:nd-1).^2+err(2:nd).^2; dchi = sqrt(sum(dres.^2./dvar)/nd); fprintf('RealSigma white dchi %10.6f\n',dchi); err = err*dchi; % Loop over the averaging intervals for i = 1:numav dyr = i*minav/365.25; % Compute the averaged residals for this interval chi2 = 0; num = 0; avR = []; avS = []; for t = start:dyr:stop sel = find(times >= t & times < t+dyr); ars = res(sel); ass = err(sel); if length(ars) > 2 w = 1./ass.^2; Wmean = w*ars'/sum(w); Werr = 1/sum(w); avR = [avR Wmean]; avS = [avS Werr]; end end % Fit we have more than 2 values get the chi^2 num = length(avR); if num > 2 savsta = [savsta sum(avR.^2./avS)/num]; timsta = [timsta i*minav]; end end % For the moment write the values %for i = 1:numav % fprintf(1,'Stats: %6.1f day, Chi^2 %10.2f\n',minav*i,savsta(i)); %end % % Now do the fitting taus = [1 2 4 8 16 32 64 128 256 512 1024 2048 5096 10192]; ntau = length(taus); nc = length(savsta); for i = 1:ntau ef = (1-exp(-timsta/taus(i))); alpha(i) = mean(savsta./ef); res = savsta-(alpha(i)*ef); rmsa(i) = std(res); end [mnr, ir] = min(rmsa); nrmsR = sqrt(alpha(ir))*dchi; fprintf(1,'NRMS Realistic %7.2f; Correlation time %8.2f days\n',nrmsR,taus(ir)); % ef = (1-exp(-timsta/taus(ir))); fit = ef*alpha(ir); % See if we should display the fit he = findobj(gcf,'Tag','DisplayFit'); DisplayFit = get(he,'Value'); if DisplayFit figure; ht = plot(timsta,savsta,'rs',timsta,fit,'g^-'); %set(gca,'Title','Chi^2 Fit as function of averaging time', ... % 'XLabel','Averging time (days)', ... % 'YLabel','Chi^2 of Residuals') title(strcat(PlotTitle,' Chi^2 Fit as function of averaging time')); xlabel('Averging time (days)'); ylabel('Chi^2 of Residuals'); end maci/tsview/checkbr.m 0000644 0000765 0000765 00000002527 11155777336 013240 0 ustar tah tah function unbrk = checkbr(breaks) % CHECKBR Check for duplicate times with breaks % unbrk = []; if ~isempty(breaks) nb = length(breaks(:,1)); ed = ones(1,nb); for i = 1:nb-1 for j = i+1:nb if breaks(j,1) == breaks(i,1) && breaks(j,2) == breaks(i,2) % If one of these is stanard break remove that one if ed(j) == 1 ed(j) = 0; if breaks(j,2) == 0 && breaks(j,6) == 0 fprintf(1,' Removing duplicate regular break %d Date %4d %2d %2d\n',j,YrToDate(breaks(j,1),3)) elseif breaks(j,2) == 0 && breaks(j,6) == 3 fprintf(1,' Removing duplicate earthquake break %d Date %4d %2d %2d\n',j,YrToDate(breaks(j,1),3)) elseif breaks(j,2) > breaks(j,1) fprintf(1,' Removing duplicate exponent break %d Date %4d %2d %2d tau %f6.2 days \n', ... j,YrToDate(breaks(j,1),3),(breaks(j,2)-breaks(j,2))*365 ) else fprintf(1,' Removing duplicate log break %d Date %4d %2d %2d tau %6.2f days \n', ... j,YrToDate(breaks(j,1),3),breaks(j,2)*365 ) end end end end end ed = logical(ed); unbrk = breaks(ed,:); end maci/tsview/datenuml.m 0000644 0000765 0000765 00000007356 11155777336 013455 0 ustar tah tah function n = datenuml(y,mo,d,h,mi,s) %DATENUM Serial date number. % N = DATENUM(S) converts the string S into a serial date number. % Date numbers are serial days where 1 corresponds to 1-Jan-0000. % The string S must be in one of the date formats 0,1,2,6,13,14, % 15,16 (as defined by DATESTR). Date strings with 2 character years % are interpreted to be within the 100 years centered around the % current year. % % N = DATENUM(S,PIVOTYEAR) uses the specified pivot year instead % of the current year minus 50 years. % % N = DATENUM(Y,M,D) returns the serial date number for % corresponding elements of the Y,M,D (year,month,day) arrays. % Y,M, and D must be arrays of the same size (or any can be a scalar). % % N = DATENUM(Y,M,D,H,MI,S) returns the serial date number for % corresponding elements of the Y,M,D,H,MI,S (year,month,hour, % minute,second) arrays values. Y,M,D,H,MI,and S must be arrays of % the same size (or any can be a scalar). Values outside the normal % range of each array are automatically carried to the next unit (for % example month values greater than 12 are carried to years). % % Examples: % n = datenum('19-May-1995') returns n = 728798. % n = datenum(1994,12,19) returns n = 728647. % n = datenum(1994,12,19,18,0,0) returns n = 728647.75. % % See also NOW, DATESTR, DATEVEC. % Author(s): C.F. Garvin, 3-20-95 % Copyright (c) 1984-98 by The MathWorks, Inc. % $Revision: 1.14 $ $Date: 1998/06/09 16:33:31 $ if nargin == 0 error('Not enough input arguments.'); elseif nargin == 1 if isstr(y) c = datevec(y); y=c(:,1);mo=c(:,2);d=c(:,3);h=c(:,4);mi=c(:,5);s=c(:,6); else n = y; return end elseif nargin == 2 if isstr(y) c = datevec(y,mo); % datenum(y,pivotyear) y=c(:,1);mo=c(:,2);d=c(:,3);h=c(:,4);mi=c(:,5);s=c(:,6); else error('Wrong number of input arguments.'); end elseif nargin ==4 | nargin == 5 error('Wrong number of input arguments.'); elseif nargin == 3 % Do three-way scalar expansion. if length(y)==1, y = y(ones(size(d))); end if length(mo)==1, mo = mo(ones(size(y))); end if length(d)==1, d = d(ones(size(mo))); end if length(y)==1, y = y(ones(size(d))); end % Required to be sure. h = zeros(size(y));mi = h;s = h; end if nargin == 6 % Do six-way scalar expansion. if length(y)==1, y = y(ones(size(s))); end if length(mo)==1, mo = mo(ones(size(y))); end if length(d)==1, d = d(ones(size(mo))); end if length(h)==1, h = h(ones(size(d))); end if length(mi)==1, mi = mi(ones(size(h))); end if length(s)==1, s = s(ones(size(mi))); end % Required to be sure. if length(y)==1, y = y(ones(size(s))); end if length(mo)==1, mo = mo(ones(size(y))); end if length(d)==1, d = d(ones(size(mo))); end if length(h)==1, h = h(ones(size(d))); end if length(mi)==1, mi = mi(ones(size(h))); end end sizes = [size(y);size(mo);size(d);size(h);size(mi);size(s)]; row_index = find(sizes(:,1) ~= sizes(1,1) & sizes(:,1) ~= 1); col_index = find(sizes(:,2) ~= sizes(2,2) & sizes(:,2) ~= 1); if ~isempty(row_index) | ~isempty(col_index) error('Y,M,D,H,MI,and S must all be the same size.') end % Make sure mo is in the range 1 to 12. mo(mo==0) = 1; y = y + floor((mo-1)/12); mo = rem(mo-1,12)+1; % Running total of days per month cumdpm = cumsum([0;31;28;31;30;31;30;31;31;30;31;30;31]); % result = (365 days/year)*(number of years) + number of leap years ... % + days in previous months + days in this month + fraction of a day. n = 365.*y + ... % Convert year, month, day to date number ceil(y/4)-ceil(y/100)+ceil(y/400) + reshape(cumdpm(mo),size(mo)) + ... ((mo > 2) & ((rem(y,4) == 0 & rem(y,100) ~= 0) | rem(y,400) == 0)) + d; if any(h~=0) | any(mi~=0) | any(s~=0) n = n + (h.*3600+mi.*60+s)./86400; end maci/tsview/DateToYr.m 0000644 0000765 0000765 00000000770 11155777336 013330 0 ustar tah tah function Yr = DateToYr(date,n) % Yr = DateToYr(date,n): Converts calender Date to fractional year % ser = datenuml(date(1),date(2),date(3)); serbeg = datenuml(date(1),1,1); % See if year and therefore how many days in year if mod(date(1),4) == 0 & mod(date(1),100) ~= 0 | mod(date(1),400 )== 0 loy = 366; else loy = 365; end Yr = date(1)+(ser-serbeg)/loy; if n >= 4, Yr = Yr + date(4)/(24*loy); end if n >= 5, Yr = Yr + date(5)/(60*24*loy); end if n >= 6, Yr = Yr + date(6)/(60*60*24*loy); end maci/tsview/datevecl.m 0000644 0000765 0000765 00000012074 11155777336 013424 0 ustar tah tah function [y,mo,d,h,mi,s] = datevecl(t,pivotyear) %DATEVEC Date components. % C = DATEVEC(T) separates the components of date strings and date % numbers into date vectors containing [year month date hour mins % secs] as columns. If T is a date string, it must be in one of the % date formats 0,1,2,6,13,14, 15,16 (as defined by DATESTR). Date % strings with 2 character years are interpreted to be within the 100 years % centered around the current year. % % [Y,M,D,H,MI,S] = DATEVEC(T) returns the components of the date % vector as individual variables. % % [...] = DAVEVEC(T,PIVOTYEAR) uses the specified pivot year instead of the % current year minus 50 years. % % Examples % d = '12/24/1984'; % t = 725000.00; % c = datevec(d) or c = datevec(t) produce c = [1984 12 24 0 0 0]. % [y,m,d,h,mi,s] = datevec(d) returns y=1984, m=12, d=24, h=0, mi=0, s=0. % c = datevec('5/6/03') produces c = [2003 5 6 0 0 0] until 2054. % c = datevec('5/6/03',1900) produces c = [1903 5 6 0 0 0]. % % See also DATENUM, DATESTR, CLOCK. % Copyright (c) 1984-98 by The MathWorks, Inc. % $Revision: 1.18 $ $Date: 1998/06/09 16:33:32 $ if nargin < 1 error('Not enough input arguments.'); end if nargin<2 clk = clock; pivotyear = clk(1)-50; % Pivot year is (current year - 50 years) end [date_row,date_col] = size(t); if isstr(t) pm = -1; % means am or pm is not in datestr dts = zeros(date_row,6); siz = [date_row 1]; for count = 1:date_row % Convert date input to date vector % Initially, the six fields are all unknown. c(1,1:6) = NaN; d = [' ' lower(t(count,:)) ' ']; % Replace 'a ', 'am', 'p ' or 'pm' with ': '. p = max(find(d == 'a' | d == 'p')); if ~isempty(p) if (d(p+1) == 'm' | d(p+1) == ' ') & d(p-1) ~= lower('e') pm = (d(p) == 'p'); if d(p-1) == ' ' d(p-1:p+1) = ': '; else d(p:p+1) = ': '; end end end % Any remaining letters must be in the month field; interpret and delete them. p = find(isletter(d)); if ~isempty(p) k = min(p); if d(k+3) == '.', d(k+3) = ' '; end M = ['jan'; 'feb'; 'mar'; 'apr'; 'may'; 'jun'; ... 'jul'; 'aug'; 'sep'; 'oct'; 'nov'; 'dec']; c(2) = find(all((M == d(ones(12,1),k:k+2))')); d(p) = setstr(' '*ones(size(p))); end % Find all nonnumbers. p = find((d < '0' | d > '9') & (d ~= '.')); % Pick off and classify numeric fields, one by one. % Colons delinate hour, minutes and seconds. k = 1; while k < length(p) if d(p(k)) ~= ' ' & d(p(k)+1) == '-' f = str2double(d(p(k)+1:p(k+2)-1)); k = k+1; else f = str2double(d(p(k)+1:p(k+1)-1)); end if ~isnan(f) if d(p(k))==':' | d(p(k+1))==':' if isnan(c(4)) c(4) = f; % hour if pm == 1 & f ~= 12 % Add 12 if pm specified and hour isn't 12 c(4) = f+12; elseif pm == 0 & f == 12 c(4) = 0; end elseif isnan(c(5)) c(5) = f; % minutes elseif isnan(c(6)) c(6) = f; % seconds else error(['Too many time fields in ' t]) end elseif isnan(c(2)) if f > 12 error([num2str(f) ' is too large to be a month.']) end c(2) = f; % month elseif isnan(c(3)) c(3) = f; % date elseif isnan(c(1)) if (f >= 0) & (p(k+1)-p(k) == 3) % two char year % Moving 100 year window centered around current year c(1) = pivotyear + rem(f + 100 - rem(pivotyear,100),100); else c(1) = f; % year end else error(['Too many date fields in ' t]) end end k = k+1; end if sum(isnan(c)) >= 5 error(['Cannot parse date ' t]) end % If the any of the day fields have been set, set an unspecified % year to the current year if isnan(c(1)) & any(~isnan(c(2:3))), clk = clock; c(1) = clk(1); end % If any field has not been specified, set it to zero. p = find(isnan(c)); if ~isempty(p) c(p) = zeros(1,length(p)); end dts(count,:) = c; end c = dts; else siz = size(t); c = dvcorel(86400*t); end % Make sure time part is properly rounded, the day number is within % range, and the first five fields are integers. maxc = ones(size(c,1),1)*[24 60 60]; [e,col] = find(any((c(:,4:6) >= maxc)')' | ... any((c(:,3) > eomday(c(:,1),c(:,2)))')' | ... any((c(:,1:5) ~= floor(c(:,1:5)))')'); if ~isempty(e), dn = datenum(c(e,1),c(e,2),c(e,3),c(e,4),c(e,5),c(e,6)); t = datevec(dn); if dn < 1, % Time only c(e,4:6) = t(:,4:6); else c(e,:) = t; end end if nargout <= 1 y = c; else y = reshape(c(:,1),siz); mo = reshape(c(:,2),siz); d = reshape(c(:,3),siz); h = reshape(c(:,4),siz); mi = reshape(c(:,5),siz); s = reshape(c(:,6),siz); end maci/tsview/dvcorel.m 0000644 0000765 0000765 00000002671 11155777336 013275 0 ustar tah tah function c = dvcorel(t) %DVCORE Core routine for DATEVEC. % D = DVCORE(T) converts the time in seconds into a date vector. % % Called by DATEVEC and DATESTR. % Copyright (c) 1984-98 by The MathWorks, Inc. % $Revision: 1.4 $ $Date: 1998/04/07 18:36:00 $ t = t(:); % Make t a column vector len = length(t); if any(~isfinite(t)), error('Date numbers must be finite.'); end % If necessary, make t positive. % Shift by an integer multiple of 400 years, which is p seconds. p = (400*365+97)*86400; shift = max(0,ceil(-t ./ p)); t = t + shift*p; % Convert seconds into a date vector components. t is in seconds. c(:,6) = rem(t,60); t = floor(t); % Truncate to the nearest second t = floor(t ./ 60); % minutes. c(:,5) = rem(t,60); t = floor(t ./ 60); % hours. c(:,4) = rem(t,24); t = floor(t ./ 24); % days. a = 365+97/400; y = floor(t ./ a); % year i = find(t <= 365*y + ceil(0.25*y)-ceil(0.01*y)+ceil(0.0025*y)); if ~isempty(i), y(i) = y(i) - 1; end c(:,1) = y - shift*400; dpm = [31 28 31 30 31 30 31 31 30 31 30 31]; dpm = dpm(ones(len,1),:); t = t - (365*y + ceil(0.25*y)-ceil(0.01*y)+ceil(0.0025*y)); i = find((rem(y,4) == 0 & rem(y,100) ~= 0) | rem(y,400) == 0); dpm(i,2) = 29+zeros(size(i)); cdm = cumsum(dpm')'; % dpm rows may now differ cdm = [zeros(len,1) cdm(:,1:11)]; c(:,2) = sum((t(:,ones(1,12)) - cdm > 0)')'; c(:,3) = t - cdm(len*(c(:,2)-1)+(1:len)'); maci/tsview/EarthScope_BoiseUAL0509.pdf 0000644 0000765 0000765 00000216465 11161162173 016202 0 ustar tah tah %PDF-1.3 % 4 0 obj << /Length 5 0 R /Filter /FlateDecode >> stream x\YoV~2/W( 8^:$qj+M<(2k"K$'?9wbh %Ӽ=v⋑(OϢh/Ók_ԹtS? 9~1$ˆ! Gx>~0^y;,giGH\_28-01Go?oNzW' `?rY 7w~ܿ I'7oydO'C F^3t-1mxc+d\ذP ̪ C<Ej\"F:w #bK|m"@*Z3)r[ᶲp W;pGQ28e/$