%  runNikoCoordModelSopac.m


%  Operates on *.neu files (format from createSeries) to either
%  clean or model the timeseries.
 
%  'clean' - read *r.neu (raw) data, estimate model parameters, 
%            clean resids using IQR criteria, output *c.neu files
%            (parameter output files and plots optional)
 
%  'model' - read *c.neu (clean) data, estimate model parameters,
%            compute model vector and model position at epoch, 
%            rotate model position at epoch and rates to xyz 
%            and output to file (parameter output files and 
%            plots optional)

%  calls clean.m, coord_model.m, runNikoCoordModelSopac_plot.m

%  by Rosanne Nikolaidis, April 2001

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set defaults and get user inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; clear; 

run_type = input('run type (in single quotes, clean/model) :   ');
if run_type == 'clean',
   criteria = input('edit criteria for IQR cleaning:    ');
   window = input('window size (days):    ');
end
site_list = input('site or site list (in single quotes):   ');
stack_flag = input('subtract stack file? [0 1]:   ');
offset_list = input('offset list (in single quotes):    ');
rate_inp = input('solve for rate [1] OR enter apr xyz rate file (in single quotes):    ');
ann_inp = input('annual flag [0 1] OR annual list (in single quotes):   ');
semi_inp = input('semi-annual flag [0 1] OR semi-annual list (in single quotes):   ');
tau_inp = input('decay list (in single quotes) or 0:   ');
vbreak_inp = input('vbreak list (in single quotes) or 0:   ');
output_flag = input('output model parameters? [0 1]:    ');
outputdata_flag = input('output residual timeseries and filtered data? [0 1]:    ');
plot_flag = input('make plots? [0 1]:    ');
noise_flag = input('use data unc. [0], white noise model [1] or white+flicker noise model? [2]:     ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read in lists
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(site_list) == 4
  sites = site_list;
else
  [sites] = textread(site_list,'%4c','commentstyle','shell');
  array_name = site_list(1:4);
end
if ischar(offset_list)
   [sites_offset,dir_offset,yr_offset] = textread(offset_list,'%4c %3c %f %*d','commentstyle','shell');
end
if ischar(ann_inp)
   [sites_annual,dir_annual] = textread(ann_inp,'%s %s','commentstyle','shell');
   sites_annual=char(sites_annual);  dir_annual=char(dir_annual);
end
if ischar(semi_inp)
   [sites_semi,dir_semi] = textread(semi_inp,'%s %s','commentstyle','shell');
   sites_semi=char(sites_semi);  dir_semi=char(dir_semi);
end
if ischar(tau_inp)
   [sites_tau,dir_tau,time_tau,tau] = textread(tau_inp,'%s %s %f %f','commentstyle','shell');
   sites_tau=char(sites_tau);  dir_tau=char(dir_tau);
end
if ischar(vbreak_inp)
   [sites_vbreak,dir_vbreak,vbreak] = textread(vbreak_inp,'%s %s %f','commentstyle','shell');
   sites_vbreak=char(sites_vbreak);  dir_vbreak=char(dir_vbreak);
end
if rate_inp~=1 
  [sites_xyzrate,xrates,yrates,zrates] = textread(rate_inp,'%4c %f %f %f','commentstyle','shell');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read in data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for s = 1:length(sites(:,1))  % BEGIN BIG LOOP OVER SITES
site = sites(s,:)

if run_type == 'clean',
   file=[site 'r.neu'];
else 
   file=[site 'c.neu'];
end

clear neu_flt neusig_flt
clear tneusig time neu neusig resids resids_c wrms_pf chi2_pf rms_pf
clear neu_epoch neu_rate dneu_epoch dneu_rate

if exist(file,'file') ~= 0      % Check to see if the file exists, if not, end

[header] = textread(file,'%s',12,'whitespace','\b\r\n\t');
[t,n,e,u,ns,es,us] = textread(file,'%f %f %f %f %f %f %f','commentstyle','shell');
tneusig = [t n e u ns es us];
clear t n e u ns es us

sizeneusig = size(tneusig);

% pj, 07/29/2004: model fails if file has 4 or less points.
% after increasing this by one each time it fails, increase to 10 now
if sizeneusig(1) > 4   % Check length of file, if too short, end 

% units are mm and years
time=tneusig(:,1);
time1=time(1);
neu=[tneusig(:,2), tneusig(:,3), tneusig(:,4)].*1000;
neusig=[tneusig(:,5), tneusig(:,6), tneusig(:,7)].*1000;
ndata = length(neu)

rate_inp = 1;
if  (max(time) - min(time)) < 0.1    % If less than 3 months long, set rate = 0
    rate_inp = 0;
end

if noise_flag == 2
% construct flicker noise matrix using Hosking method (by Hadley and Simon)
  [trash,j1] = power_law_cov(time,1,-1);
end

flt_flag = 'unf';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if input xyz rates are to be assumed, rotate to neu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if rate_inp~=1
%  k = strmatch(site,sites_xyzrate);
%  if length(k) ==0       
%    [site ' not found in xyzrate file'] 
%  else
%    xrate_apr=xrates(k)*1000;
%    yrate_apr=yrates(k)*1000;
%    zrate_apr=zrates(k)*1000;
%  end
    xrate_apr=0.0;
    yrate_apr=0.0;
    zrate_apr=0.0;
end

% use refXYZ from neu file header to compute jac 
clear iX iY iZ refX refY refZ jac neurate_apr 
[refname, refvalue] = textread(file,'%*s %s %*s %f',5,'headerlines',3);
iX = strmatch('Reference_X',refname);
iY = strmatch('Reference_Y',refname);
iZ = strmatch('Reference_Z',refname);
refX=refvalue(iX);
refY=refvalue(iY);
refZ=refvalue(iZ);
jac = xyzJacobian(refX,refY,refZ);

iLat = strmatch('Latitude',refname);
iLon = strmatch('Longitude',refname);
Lat=refvalue(iLat);
Lon=refvalue(iLon);

if rate_inp~=1
   neurate_apr = jac * [xrate_apr;yrate_apr;zrate_apr];
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for comp = 1:3              % BEGIN BIG LOOP OVER NEU 

if comp == 1
   neu_flag = 'n'; neu_txt = 'north'
elseif comp == 2
   neu_flag = 'e'; neu_txt='east'
elseif comp == 3
   neu_flag = 'u'; neu_txt='up'
end   

clear data sig
data=neu(:,comp);
sig = neusig(:,comp);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if input rates assumed, detrend data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if rate_inp==1
   rate_flag = 1;   % solve for rates
else 
   rate_flag = 0;   % assume input rates
   data = data - neurate_apr(comp) .*time;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subtract the stack data to create filtered data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if stack_flag ~= 0 
  clear idata
  clear istack
  flt_flag = 'flt';
% read in stack and find intersection times
  clear idata istack stack
  file=['stack_' array_name '_' neu_flag '.dat'];
  eval(['load ' file])
  array=['stack_' array_name '_' neu_flag];
  stack=[eval(array)];
  [cc,idata,istack] = intersect(time,stack(:,1));
% subtract stack from data

  data_flt = data;
  data_flt(idata) = data(idata) - stack(istack,2).*1000;
  sig_flt = sig;
  sig_flt(idata) = sqrt(sig(idata).^2 + (stack(istack,3).*1000).^2);
  data=data_flt;
  sig = sig_flt;

  neu_flt(:,comp) = data_flt;
  neusig_flt(:,comp) = sig_flt;

end   %END stack_flag section

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compile model parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jumps=[];
i = strmatch(site,sites_offset);
for a = 1:length(i) 
   if findstr(neu_flag,dir_offset(i(a),:)) & (yr_offset(i(a)) > min(time)) & (yr_offset(i(a)) < max(time))      
       jumps = [jumps  yr_offset(i(a))];
   end
end
njumps=length(jumps);

if ischar(ann_inp)
  j = strmatch(site,sites_annual);           %check for site match
  if findstr(neu_flag,dir_annual(j,:))       %check for dir match
     ann_flag = 1;
  else
     ann_flag = 0;
  end
else
  if ((max(time) - min(time)) >= 2.0 ) & (ndata >= 360)  % if more than 2 years of data, annual 
     ann_flag = ann_inp;                                 % term is estimated 
  else
    ann_flag = 0;
  end
end

if ischar(semi_inp)
  j = strmatch(site,sites_semi);           %check for site match
  if findstr(neu_flag,dir_semi(j,:))       %check for dir match
     semi_flag = 1;
  else
     semi_flag = 0;
  end
else
  if ((max(time) - min(time)) >= 1.0) & (ndata >= 180) % if more than 1 year of data, semi-annual 
     semi_flag = semi_inp;                             % term is estimated 
  else
     semi_flag = 0;
  end
end

TAUs = [];
epoch_tau = [];
if ischar(tau_inp)
  q = strmatch(site,sites_tau);               %check for site match, q is vector of indices in decay list
  for a = 1:length(q)                         % a is number of decays to model
    if findstr(neu_flag,dir_tau(q(a),:))        %check for dir match
      TAUs = [TAUs tau(q(a))];
      epoch_tau = [epoch_tau time_tau(q(a))];
    end
  end
end
 
VBREAKs = [];
if ischar(vbreak_inp)
  j = strmatch(site,sites_vbreak);           %check for site match
  for a = 1:length(j) 
    if findstr(neu_flag,dir_vbreak(j(a),:))       %check for dir match
       VBREAKs = [VBREAKs vbreak(j(a))];
    end
  end
end

% compute degrees-of-freedom
nvb=length(VBREAKs);
ntau=length(TAUs);
dof=1+rate_flag+2*ann_flag+2*semi_flag+njumps+nvb+ntau;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% construct data noise covariance from noise model or data unc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear c 
no_noiseamp = 0;
if noise_flag == 0           %use data uncertainties for white noise covariance
  c = diag(sig.^2);
elseif noise_flag == 1       %use white noise model for covariance
  wh_noise_file=['WH_' flt_flag '_fn.dat'];
  [sites_noise,wh] = textread(wh_noise_file,'%4c %*13c %*4c %*s %*s %*s %f %*s %*f','commentstyle','shell');
  m = strmatch(site,sites_noise);           %check for site match
  wh_amp = wh(m(comp));
  c = wh_amp^2 .* diag(ones(ndata,1));
elseif noise_flag == 2      %use white+flicker noise model for covariance
  wh_noise_file=['WH_' flt_flag '_fn.dat'];
  [sites_noise,wh] = textread(wh_noise_file,'%4c %*13c %*4c %*s %*s %*s %f %*s %*f','commentstyle','shell');
  if strmatch(site,sites_noise)                %check for site match
     m = strmatch(site,sites_noise);
     wh_amp = wh(m(comp));
     pl_noise_file=['PL_' flt_flag '_fn.dat'];
     [sites_noise,pl] = textread(pl_noise_file,'%4c %*13c %*4c %*s %*s %*s %f %*s %*f','commentstyle','shell');
     m = strmatch(site,sites_noise);           %check for site match
     pl_amp = pl(m(comp));
     c = wh_amp^2 .* diag(ones(ndata,1)) + pl_amp^2 .* j1;
   else
     no_noiseamp = 1;
     c = diag(sig.^2);
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% estimate model parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear fit cov_fit data_hat fit_sig
[fit,cov_fit,data_hat,chi_2_misfit]=coord_model(time,data,c,rate_flag,ann_flag,semi_flag,[jumps],[epoch_tau],[TAUs]./365,[VBREAKs]);

%compute postfit statistics
resids(:,comp)=data-data_hat;
if noise_flag == 2
  chi2_pf(comp) = 1.000;
else
  [wrms_pf(comp),chi2_pf(comp)]=wstat2(resids(:,comp),sqrt(diag(c)),dof);
end
rms_pf(comp) = std(resids(:,comp));
% scale fit covariance by chi2
cov_fit = cov_fit.*chi2_pf(comp);
fit_sig=sqrt(diag(cov_fit));

fit=[fit fit_sig];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% extract fit parameters for output 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear yint slope dslope sine cosine offset decay

yint=fit(1,:); 
skip=1;
if rate_flag==1
  slope(1,:)=fit(skip+1,:); 
  skip=skip+1;
else
  slope(1,:)=[neurate_apr(comp) 0];
end
for k=1:nvb
    slope(k+1,:)=fit(skip+k,:);               
    dslope(k,1)=slope(k+1,1) - slope(k,1);
    dslope(k,2)= sqrt(cov_fit(skip+k,skip+k) + cov_fit(skip+k-1,skip+k-1) - 2.* cov_fit(skip+k,skip+k-1));
end
skip=skip+nvb;
if ann_flag 
   sine=fit(skip+1,:); 
   cosine=fit(skip+2,:); 
   skip=skip+2;
%   phase = atan2(-sine(1),cosine(1));
   phase = -(atan2(-sine(1),cosine(1)));
   if phase<0
      phase=phase+2*pi;
   end
   annual(1) =  sqrt(sine(1).^2 + cosine(1).^2);
%   annual(1) =  cosine(1)./cos(phase);    % also, annual = sqrt(sine(1).^2 + cosine(1).^2)
   annual(2) =  sqrt(sine(2).^2 + cosine(2).^2);
else 
   sine=[0 0];   
   cosine=[0 0];
end

if semi_flag 
   sine_semi=fit(skip+1,:); 
   cosine_semi=fit(skip+2,:); 
   skip=skip+2;
   phase_semi = -(atan2(-sine_semi(1),cosine_semi(1)));
   if phase_semi<0
      phase_semi=phase_semi+2*pi;
   end
   semi(1) =  sqrt(sine_semi(1).^2 + cosine_semi(1).^2);
   semi(2) =  sqrt(sine_semi(2).^2 + cosine_semi(2).^2);
else 
   sine_semi=[0 0];   
   cosine_semi=[0 0];
end

if jumps
  for k=1:njumps
     offset(k,:)=fit(skip+k,:);
  end
else
  offset=[0 0];
end
skip=skip+njumps;

if TAUs
  for k=1:ntau
       decay(k,:)=fit(skip+k,:);
  end
else
  decay=[0 0];
end

% subtract decay amplitude from total offset at eqs to compute coseismic offset
[cdecay,ij,id]=intersect(jumps,epoch_tau);
offset(ij,1) = offset(ij,1) + decay(id,1);
offset(ij,2) = sqrt(offset(ij,2).^2 + decay(id,2).^2);

% find offsets at earthquakes
eqs = [2002.8425 2004.9850 2005.2383 2005.2720 2005.2890 2005.3000 2005.3658 2005.5082 2007.1795 2007.6986 2007.7013];
[ceqs,ij,ieqs] = intersect(jumps,eqs);
neqs = length(ij);
offset_eq = offset(ij,:);
jumps_eq = eqs(ieqs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if run_type == 'model', compute model vector and epoch position
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if run_type == 'model',

clear time_m model model_rate_jumps model_norate_nojumps

% pj, 03/31/2003: use 1/365 interval
time_m = time1:0.0001:max(time);
%time_m = time1:0.00274:max(time);

model = yint(1) + slope(1,1).*time_m + sine(1).*sin(time_m*2*pi) + cosine(1).*cos(time_m*2*pi);
model = model +  sine_semi(1).*sin(time_m*4*pi) + cosine_semi(1).*cos(time_m*4*pi);
model_rate_jumps = yint(1) + slope(1,1).*time;
model_norate_nojumps = sine(1).*sin(time_m*2*pi) + cosine(1).*cos(time_m*2*pi);
model_norate_nojumps = model_norate_nojumps + sine_semi(1).*sin(time_m*4*pi) + cosine_semi(1).*cos(time_m*4*pi);
for k=1:nvb
  model = model + (time_m>=VBREAKs(k)) .* dslope(k,1) .* (time_m-VBREAKs(k)); 
  model_norate_nojumps = model_norate_nojumps + (time_m>=VBREAKs(k)) .* dslope(k,1) .* (time_m-VBREAKs(k)); 
end

if nvb < 1
   numdays(1) = ndata;
else
   mask = ones(size(time));
  [C,numdays(1),IA] = intersect( time<=VBREAKs(1),mask);
  for i=1:nvb-1,
      [C,IB,IA] = intersect(time<=VBREAKs(i+1), mask);
      [C,ID,IA] = intersect(time<=VBREAKs(i), mask);
      numdays(i+1) = IB - ID + 1;
  end
      [C,IB,IA] = intersect( time<VBREAKs(nvb),mask);
      numdays(nvb+1) = length(time) - IB;
end

for k=1:njumps
  model = model + (time_m>=jumps(k)) .* offset(k,1);
  model_rate_jumps = model_rate_jumps + (time>=jumps(k)) .* offset(k,1);
end
for k=1:ntau
  model = model + (time_m>=epoch_tau(k)) .* (decay(k,1).* exp(-(time_m-epoch_tau(k))./(TAUs(k)./365)) - decay(k,1));
  model_norate_nojumps = model_norate_nojumps + (time_m>=epoch_tau(k)) .* (decay(k,1).* exp(-(time_m-epoch_tau(k))./(TAUs(k)./365)) - decay(k,1));
end
  

if outputdata_flag == 1,

% pj, 03/31/2003: write model data (includes jumps and trend) to output file
   xml_outfile=[site '_data_model_' neu_flag '_' flt_flag '.xml'];
   xml_out = fopen(xml_outfile,'w');

%   header = ['<?xml version="1.0" standalone="no"?><!DOCTYPE plot PUBLIC "-//UC Berkeley//DTD PlotML 1//EN" "http://sopac.ucsd.edu/plotml/plotML_1.dtd"><plot><size width="590" height="210"><title>' site ' ' neu_txt ' ' flt_flag ' data v model</title><yLabel>position (mm)</yLabel><xLabel>epoch (yr)</xLabel><noGrid/>'];

   printTitle = '';
   VBREAKs_out=[time1 VBREAKs max(time)];
   for t=1:(1+nvb)
     txt = sprintf('%3.1f +- %2.1f (%6.2f:%6.2f)',slope(t,:),VBREAKs_out(t), VBREAKs_out(t+1));
     printTitle = [printTitle ' ' txt];  
   end

   rms = sprintf('rms: %3.1f mm', rms_pf(comp));
   printTitle = [printTitle '; ' rms];

   header = ['<?xml version="1.0" standalone="no"?><!DOCTYPE plot PUBLIC "-//UC Berkeley//DTD PlotML 1//EN" "http://sopac.ucsd.edu/plotml/plotML_1.dtd"><plot><size width="590" height="210"><title>' site ' ' neu_txt ' slope(s) (mm/yr):' printTitle '</title><yLabel>position (mm)</yLabel><xLabel>epoch (yr)</xLabel><noGrid/>'];

   modelHeader = '<default connected="no" marks="points"/><dataset connected="yes" marks="none">';

   datasetFooter = '</dataset>';

   fprintf (xml_out, '%s\n', header);
   fprintf (xml_out, '%s\n', modelHeader);

   for i = 1:length(time_m)
     fprintf(xml_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',time_m(i), model(i));
   end

   fprintf (xml_out, '%s\n', datasetFooter);

% pj: we need to close this file here since it requires re-opening later 
% (open by component-specific name in data section below)
   fclose(xml_out);


% pj, 03/31/2003: write norate nojumps model data to output file
   xml_nr_nj_outfile=[site '_data_model_nr_nj_' neu_flag '_' flt_flag '.xml'];
   xml_nr_nj_out = fopen(xml_nr_nj_outfile,'w');

   printTitle = '';
   VBREAKs_out=[time1 VBREAKs max(time)];
   for t=1:(1+nvb)
     txt = sprintf('%3.1f +- %2.1f (%6.2f:%6.2f)',slope(t,:),VBREAKs_out(t), VBREAKs_out(t+1));
     printTitle = [printTitle ' ' txt];  
   end

   rms = sprintf('rms: %3.1f mm', rms_pf(comp));
   printTitle = [printTitle '; ' rms];

   header = ['<?xml version="1.0" standalone="no"?><!DOCTYPE plot PUBLIC "-//UC Berkeley//DTD PlotML 1//EN" "http://sopac.ucsd.edu/plotml/plotML_1.dtd"><plot><size width="590" height="210"><title>' site ' ' neu_txt ' slope(s) (mm/yr):' printTitle '</title><yLabel>position (mm)</yLabel><xLabel>epoch (yr)</xLabel><noGrid/>'];

   modelHeader = '<default connected="no" marks="points"/><dataset connected="yes" marks="none">';

   datasetFooter = '</dataset>';

   fprintf (xml_nr_nj_out, '%s\n', header);
   fprintf (xml_nr_nj_out, '%s\n', modelHeader);

% need to start connected plots with a disconnected point, so multiple plots
% do not connect
   for i = 1:1
     fprintf(xml_nr_nj_out,'<m x=\"%9.4f\" y=\"%7.3f\"/>\n',time_m(i), model_norate_nojumps(i));
   end

   for i = 1:length(time_m)
     fprintf(xml_nr_nj_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',time_m(i), model_norate_nojumps(i));
   end

   fprintf (xml_nr_nj_out, '%s\n', datasetFooter);
end
% stop pj new


% These variables will be used to compute the output xyz coordinates at the LAST epoch.

neu_epoch(comp,1) = model(length(model));
neu_rate(comp,1) = slope(nvb+1,1);
%dneu_epoch(comp,1) = yint(2);  
dneu_epoch(comp,1) = rms_pf(comp);  
dneu_rate(comp,1) = slope(nvb+1,2);

end       % end 'model' section

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% put trend back in data if detrended, detrend/dejump data for plots 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if rate_flag == 0
   data = data + neurate_apr(comp) .*time;
end

if run_type == 'model',
   data_norate_nojumps = data - model_rate_jumps;

% pj, 03/31/2003: write norate nojumps data to output file

  %dataHeader = '<dataset connected="no">';
  dataHeader = ['<dataset connected="no" name="' site '">'];

   nr_nj_file = [site '_nr_nj_' neu_flag '_' flt_flag ];   
   fid_nr_nj = fopen(nr_nj_file,'w');
   fprintf (xml_nr_nj_out, '%s\n', dataHeader);
   for i = 1:length(time)
     fprintf(xml_nr_nj_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',time(i), data_norate_nojumps(i));
     fprintf(fid_nr_nj,'%9.4f %9.4f\n',time(i),data_norate_nojumps(i));
   end

   fprintf (xml_nr_nj_out, '%s\n', datasetFooter);  
   fclose(fid_nr_nj);

% stop pj new

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write output file 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if output_flag==1,

outfile=[site '_' neu_flag '_' flt_flag '.out']; 
%outsumfile=['SugarSeriesSummary' date '.txt'];   %Linette
fid_out = fopen(outfile,'w');
%fid_sum_out = fopen(outsumfile,'a');   %Linette
%fprintf(fid_sum_out,'%-14s' \n',outfile);  %Linette

 
VBREAKs_out=[time1 VBREAKs max(time)];
fprintf(fid_out,'%-12s  %9.3f\n','postfit chi2', chi2_pf(comp));
fprintf(fid_out,'%-12s  %9.3f\n','postfit rms', rms_pf(comp));
if (noise_flag < 2) | (no_noiseamp == 1)
  fprintf(fid_out,'%-17s \n','white noise model');
else
  fprintf(fid_out,'%-25s \n','white+flicker noise model');
  fprintf(fid_out,'%-12s  %9.3f\n','white noise amp', wh_amp);
  fprintf(fid_out,'%-12s  %9.3f\n','flicker noise amp', pl_amp);
end
fprintf(fid_out,'%12s  %16.6f \n','Reference_X', refX);
fprintf(fid_out,'%12s  %16.6f \n','Reference_Y', refY);
fprintf(fid_out,'%12s  %16.6f \n','Reference_Z', refZ);
fprintf(fid_out,'%12s  %8.4f  \n','start_epoch', min(time));
fprintf(fid_out,'%12s  %8.4f  \n','end_epoch', max(time));
fprintf(fid_out,'%12s  %6i  \n','num_days', ndata);
fprintf(fid_out,'%12s  %9.3f  %9.3f\n','y-intercept', yint);
for t=1:(1+nvb)
   fprintf(fid_out,'%11s%1i %10.3f  %9.3f  %8.4f %1s %8.4f\n','slope_',t, slope(t,:), VBREAKs_out(t), '-', VBREAKs_out(t+1));
   fprintf(fid_out,'%11s%1i %6i  \n','num_days_',t, numdays(t));
end
for t=1:nvb
   fprintf(fid_out,'%11s%1i %10.3f  %9.3f  %8.4f \n','dslope_',t, dslope(t,:), VBREAKs_out(t+1));
end
if ann_flag
  fprintf(fid_out,'%12s  %9.3f  %9.3f  %8.4f %1s %8.4f\n','sine_ann', sine, min(time), '-',max(time));
  fprintf(fid_out,'%12s  %9.3f  %9.3f  %8.4f %1s %8.4f\n','cosine_ann', cosine, min(time),'-', max(time));
  fprintf(fid_out,'%12s  %9.3f  \n','  phase_ann', phase);
  fprintf(fid_out,'%12s  %9.3f  %9.3f  %8.4f %1s %8.4f\n','annual', annual, min(time),'-', max(time));
end
if semi_flag
  fprintf(fid_out,'%12s  %9.3f  %9.3f  %8.4f %1s %8.4f\n','sine_semi', sine_semi, min(time), '-',max(time));
  fprintf(fid_out,'%12s  %9.3f  %9.3f  %8.4f %1s %8.4f\n','cosine_semi', cosine_semi, min(time),'-', max(time));
  fprintf(fid_out,'%12s  %9.3f  \n','  phase_semi', phase_semi);
  fprintf(fid_out,'%12s  %9.3f  %9.3f  %8.4f %1s %8.4f\n','semi', semi, min(time),'-', max(time));
end

  % pj: write offsets to xml file
offsetHeader = '<dataset connected="no" marks="none" stems="yes">';
fprintf (xml_nr_nj_out, '%s\n', offsetHeader);

for t=1:njumps
   if intersect(jumps(t),eqs)    % if offset is coseismic, tag it
      fprintf(fid_out,'%11s%1i  %9.3f  %9.3f  %8.4f %9s\n','offset_',t, offset(t,:), jumps(t), 'coseismic');
   else
      fprintf(fid_out,'%11s%1i  %9.3f  %9.3f  %8.4f\n','offset_',t, offset(t,:), jumps(t));
   end

% pj: write offsets to xml file
% using stems function here, need to get max value of data 
   fprintf(xml_nr_nj_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',jumps(t),max(data_norate_nojumps));
end


% pj: close xml file

fprintf (xml_nr_nj_out, '%s\n', datasetFooter);
footer = '</size></plot>';
fprintf (xml_nr_nj_out, '%s\n', footer);
fclose(xml_nr_nj_out);

for t=1:ntau
   fprintf(fid_out,'%11s%1i  %9.3f  %9.3f  %8.4f  %4s%1i%2s  %-4i\n','decay_',t, decay(t,:), epoch_tau(t), 'tau_',t,' =', TAUs(t));
end

fclose(fid_out);

end    % end parameter output section

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  plots       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if plot_flag==1 & run_type ~= 'clean',
   runNikoCoordModelPlotSopac
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end            %  end big loop over neu 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if run_type == 'clean', then clean resids and output sitec.neu file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if run_type == 'clean',

%clean residuals array of outliers
clear imaj nedit resids_c tneusig_c tneusig_e
[imaj, resids_c]=var_clean(criteria,window,resids);
nedit = length(imaj)

tneusig_e = tneusig(imaj,:);
tneusig_c = tneusig; tneusig_c(imaj,:)=[];

%write cleaned neu file and file of outlier indices
itype = strmatch('# Type ',header);
header(itype)=cellstr(sprintf('%s %i %s', '# Type of file : clean ( nedits:',nedit,')'));
outfile=[site 'c.neu'];
fid_neu = fopen(outfile,'w');
for l = 1:length(header)
   fprintf(fid_neu,'%s\n', char(header(l)));
end
fprintf(fid_neu,'%9.4f %10.4f %10.4f %10.4f %7.4f %7.4f %7.4f\n', tneusig_c');
fclose(fid_neu);

%for i = 1:length(tneusig_c)
%  fprintf('%9.4f %9.4f\n',tneusig_c(i,1),tneusig_c(i,2));
%end

ioutfile = [site '_imaj.dat'];
eval(['save ' ioutfile ' imaj -ascii'])

outlierfile = [site '_outlier_' flt_flag '.dat'];
eval(['save ' outlierfile ' tneusig_e -ascii'])

% pj: "c" file written here contains output for data section of unfiltered
% model_data (with jumps and rate) xml file

% pj, 06/10/2003: commented out.  we don't use sitec.ps files,
% rather we convert the site_flt.ps and unf.ps files to gif
%if plot_flag==1,
%   runNikoCoordModelPlotSopac
%end

end            % end 'clean' section  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write out residual data and filtered neu data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output units are meters
if outputdata_flag==1,

   % pj, 05/19/2003: do we need this file?  comment out for now
   % pj, 07/14/2003: yes we do, for stacking
   resids_output = [time,resids./1000,neusig./1000];
   resids_file = [site '_resids_' flt_flag '.dat']; 
   eval(['save ' resids_file ' resids_output -ascii'])


% pj need to loop over components here
   for comp = 1:3
      if comp == 1
         neu_flag = 'n'; neu_txt = 'north'
      elseif comp == 2
         neu_flag = 'e'; neu_txt='east'
      elseif comp == 3
         neu_flag = 'u'; neu_txt='up'
      end   


      % pj output resids to xml file
      xml_resid_outfile=[site '_resids_' neu_flag '_' flt_flag '.xml'];
      xml_resid_out = fopen(xml_resid_outfile,'w');

      rms = sprintf('rms = %3.1f mm', rms_pf(comp));       
      header = ['<?xml version="1.0" standalone="no"?><!DOCTYPE plot PUBLIC "-//UC Berkeley//DTD PlotML 1//EN" "http://sopac.ucsd.edu/plotml/plotML_1.dtd"><plot><size width="590" height="210"><title>' site ' ' neu_txt ' data-model residuals: ' rms '</title><yLabel>residual (mm)</yLabel><xLabel>epoch (yr)</xLabel><noGrid/>'];
      fprintf (xml_resid_out, '%s\n', header);

      % include name here so legend is shown in ptplot
      dataHeader = ['<dataset connected="no" marks="points" name="' site '">'];
      fprintf (xml_resid_out, '%s\n', dataHeader);

      for i = 1:length(time);
         fprintf(xml_resid_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',time(i), resids(i,comp));
      end

      datasetFooter = '</dataset>';
      fprintf (xml_resid_out, '%s\n', datasetFooter);

      footer = '</size></plot>';
      fprintf (xml_resid_out, '%s\n', footer);

      fclose(xml_resid_out);


      % filtered ( is stack flag not equal to zero - interesting logic)
      if stack_flag~=0


         % pj: "f" file here contains data for data section of data of filtered
         % model_data (with jumps and rate) xml file
	
         % pj, 06/12/2003: this is the file used as input to the 
         % siteCleanFlt.neu files we create tar files of.  uncomment
         % this code

         data_flt_output = [time,neu_flt./1000,neusig_flt./1000];
         data_flt_file = [site 'f.neu']; 
         eval(['save ' data_flt_file ' data_flt_output -ascii'])

         % append to the respective file, model data should already be written
         xml_outfile=[site '_data_model_' neu_flag '_' flt_flag '.xml'];
         xml_out = fopen(xml_outfile,'a');

         % include name here so legend is shown
         dataHeader = ['<dataset connected="no" name="' site '">'];
         fprintf (xml_out, '%s\n', dataHeader);

         for i = 1:length(time);
           fprintf(xml_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',time(i), neu_flt(i,comp));
         end

         datasetFooter = '</dataset>';
         fprintf (xml_out, '%s\n', datasetFooter);

         % we will not be writing offsets later.  close file
         footer = '</size></plot>';
         fprintf (xml_out, '%s\n', footer);
         fclose(xml_out);

     % unfiltered
     else

     % stack flag is not set - this is the unfiltered run.  output cleaned data
     % we run the clean step first, but these unfiltered data must appear
     % after model data in file.

       xml_outfile=[site '_data_model_' neu_flag '_' flt_flag '.xml'];
       xml_out = fopen(xml_outfile,'a');
  
       %dataHeader = '<dataset connected="no">';
       dataHeader = ['<dataset connected="no" name="' site '">'];
       fprintf (xml_out, '%s\n', dataHeader);

       % re-create cleaned data arrays and write them now
       %clean residuals array of outliers
       clear imaj nedit resids_c tneusig_c tneusig_e

       %PJ NOTE HARDCODING
       criteria = 3;
       window = 365;
      
       [imaj, resids_c]=var_clean(criteria,window,resids);
       nedit = length(imaj)

       % pj output resids to xml file
        xml_resid_outfile=[site '_resids_' neu_flag '_' flt_flag '.xml'];
        xml_resid_out = fopen(xml_resid_outfile,'w');

        rms = sprintf('rms = %3.1f mm', rms_pf(comp));        
        header = ['<?xml version="1.0" standalone="no"?><!DOCTYPE plot PUBLIC "-//UC Berkeley//DTD PlotML 1//EN" "http://sopac.ucsd.edu/plotml/plotML_1.dtd"><plot><size width="590" height="210"><title>' site ' ' neu_txt ' data-model residuals: ' rms '</title><yLabel>residual (mm)</yLabel><xLabel>epoch (yr)</xLabel><noGrid/>'];        fprintf (xml_resid_out, '%s\n', header);
     
        % include name here so legend is shown in ptplot
        dataHeader = ['<dataset connected="no" marks="points" name="' site '">'];
        fprintf (xml_resid_out, '%s\n', dataHeader);

        time_c = time; time_c(imaj,:)=[];
        for i = 1:length(time_c);
          fprintf(xml_resid_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',time_c(i), resids_c(i,comp));
        end

        datasetFooter = '</dataset>';
        fprintf (xml_resid_out, '%s\n', datasetFooter);

        footer = '</size></plot>';
        fprintf (xml_resid_out, '%s\n', footer);

        fclose(xml_resid_out);

       % this is not used
       %tneusig_e = tneusig(imaj,:);

       tneusig_c = tneusig; tneusig_c(imaj,:)=[];
  
       useCol = comp+1;

	% pj, 09/22/2003: 
        % here linette changed from length to size.  this was giving number
	% of columns (7) instead of rows (e.g., 4) if rows were less than
	% columns.  script would fail here.
       for i = 1:size(tneusig_c,1);
         fprintf(xml_out,'<p x=\"%9.4f\" y=\"%7.3f\"/>\n',tneusig_c(i,1), tneusig_c(i,useCol).*1000.);
       end

       datasetFooter = '</dataset>';
       fprintf (xml_out, '%s\n', datasetFooter);

       % we will not be writing offsets later.  close file
       footer = '</size></plot>';
       fprintf (xml_out, '%s\n', footer);
       fclose(xml_out);

     % end stack_flag conditional
     end

  % end neu loop
  end

% end output_dataflag conditional
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if run_type == 'model', then output site_{filter_type}.xyz file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if run_type == 'model',

% compute XYZ position/velocity at epoch  (units in m and m/yr)            
% using refXYZ to compute jac and add refXYZ to xyz to get total XYZ; 
% (undoes xyz to neu conversion in /usr/local/postprocessing/src/converttoneu.c)
clear xyz_epoch XYZ_epoch XYZ_rate 
clear cov_neu_epoch cov_neu_rate cov_xyz_epoch cov_xyz_rate dxyz_epoch dxyz_rate

xyz_epoch = (jac' * neu_epoch) ./1000;
XYZ_epoch = xyz_epoch + [refX; refY; refZ]; 
XYZ_rate = (jac' * neu_rate) ./1000;

cov_neu_epoch = diag(dneu_epoch.^2);
cov_neu_rate = diag(dneu_rate.^2);
cov_xyz_epoch = jac' * cov_neu_epoch * jac;
cov_xyz_rate = jac' * cov_neu_rate * jac;
dxyz_epoch = sqrt(diag(cov_xyz_epoch))./1000;
dxyz_rate = sqrt(diag(cov_xyz_rate))./1000;

% write out xyz file at epoch
SITE=strcat(upper(site), '_GPS');
xyzfile=[site '_' flt_flag '.xyz'];
fid_xyz = fopen(xyzfile,'w');
format = '%9s %14.4f %14.4f %14.4f %7.4f %7.4f %7.4f %8.3f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\n';
fprintf(fid_xyz,format,SITE,XYZ_epoch(1),XYZ_epoch(2),XYZ_epoch(3),XYZ_rate(1),XYZ_rate(2),XYZ_rate(3),max(time_m),dxyz_epoch(1),dxyz_epoch(2),dxyz_epoch(3),dxyz_rate(1),dxyz_rate(2),dxyz_rate(3));
fclose(fid_xyz);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write output gmt file (coords and vels in gmt format)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if stack_flag
   gmtoutfile=['vel_flt.gmt']; 
   fid_gmtout = fopen(gmtoutfile,'a');
   fprintf(fid_gmtout,'%10.4f %10.4f %9.3f %9.3f %9.3f %9.3f %3.1f %4s\n', Lon,Lat,neu_rate(2),neu_rate(1),dneu_rate(2),dneu_rate(1),0.0,site);
   fclose(fid_gmtout);
else
   gmtoutfile=['vel_unf.gmt']; 
   fid_gmtout = fopen(gmtoutfile,'a');
   fprintf(fid_gmtout,'%10.4f %10.4f %9.3f %9.3f %9.3f %9.3f %3.1f %4s\n', Lon,Lat,neu_rate(2),neu_rate(1),dneu_rate(2),dneu_rate(1),0.0,site);
   fclose(fid_gmtout);
end

% End output gmt file

end       % end 'model' section

end       % end if check for whether timeseries is long enough
end       % end checking for file's existence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end            %  end big loop over sites
