% run_coord_model_psd.m

% Compute power spectral density of resid files output by run_coord_model.m
% Rosanne Nikolaidis January 2002

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get user inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; clear;
site_list = input('site or site list (in single quotes):   ');
nwin = input('number of windows to use for psd?   ');
flt_flag = input('include psd of flt data? [0 1]:   ');
output_flag = input('write output freq/power files? [0 1]:   ');

if length(site_list) == 4
  sites = site_list;
else
  [sites] = textread(site_list,'%4c','commentstyle','shell');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read in data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for s = 1:length(sites(:,1))  % BEGIN BIG LOOP OVER SITES
close all;
site = sites(s,:)
file=[site '_resids_unf.dat'];
clear tneusig time neu neusig wrms chi2 data sig f pout pout_flt
[time,n,e,u,ns,es,us] = textread(file,'%f %f %f %f %f %f %f','commentstyle','shell');
tneusig = [time n e u ns es us];
clear n e u ns es us
%neu=[tneusig(:,2), tneusig(:,3), tneusig(:,4)].*1000;
%neusig=[tneusig(:,5), tneusig(:,6), tneusig(:,7)].*1000;
neu=[tneusig(:,2), tneusig(:,3), tneusig(:,4)];
neusig=[tneusig(:,5), tneusig(:,6), tneusig(:,7)];

if flt_flag == 1
  file_flt=[site '_resids_flt.dat'];
  clear tneusig_flt time_flt neu_flt neusig_flt wrms_flt chi2_flt data_flt sig_flt
  [time_flt,n,e,u,ns,es,us] = textread(file_flt,'%f %f %f %f %f %f %f','commentstyle','shell');
  tneusig_flt = [time_flt n e u ns es us];
  clear n e u ns es us
%  neu_flt=[tneusig_flt(:,2), tneusig_flt(:,3), tneusig_flt(:,4)].*1000;
%  neusig_flt=[tneusig_flt(:,5), tneusig_flt(:,6), tneusig_flt(:,7)].*1000;
  neu_flt=[tneusig_flt(:,2), tneusig_flt(:,3), tneusig_flt(:,4)];
  neusig_flt=[tneusig_flt(:,5), tneusig_flt(:,6), tneusig_flt(:,7)];
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);
data_filled = fill_data(time,data,sig);

if flt_flag == 1
  clear data_flt sig_flt
  data_flt=neu_flt(:,comp);
  sig_flt = neusig_flt(:,comp);
  data_filled_flt = fill_data(time_flt,data_flt,sig_flt);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot psd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
nfft = length(data_filled);
window=round(nfft./nwin);
%overlap = floor(window/4);
overlap = 0;
fsamp = 1./(86400);           %sampling freq is once per day
%fref = 1./[172800 604800 1209600 2628003  7884009 15768017 31536000];
%textref = ['2d ';'1wk';'2wk';'1m ';'3m ';'6m ';'1yr'];
fref = 1./[172800 604800 1209600 2628003  7884009 15768017 31536000  63115200];
textref = ['2d ';'1wk';'2wk';'1m ';'3m ';'6m ';'1yr';'2yr'];


[p,f]=psd(data_filled,window,fsamp,window,overlap,'linear');
pout(:,comp) = p;
semilogx(f.*1000,10*log10(2*p/fsamp),'b-')
hold on 

% fit slopes for f<15d
i=find(f~=0 & f<=1/1296000);
fs=f(i(2:length(i))); ps=p(i(2:length(i)));
lfs=log10(fs); lps=10*log10(2*ps/fsamp);
zpoly=polyfit(lfs,lps,1);
slope_unf=zpoly(1)./10;
zcomp=polyval(zpoly,lfs);
semilogx(fs*1000,zcomp,'-k');
%assume stdev of data unknown and equal, so estimate it from fit
res=lps-zcomp;
chi2 = sum(res.^2)/(length(res)-2);
h=[lfs ones(length(lfs),1)];
cov=chi2.*inv(h'*h);
slope_sig_unf=sqrt(cov(1,1))./10;

if flt_flag == 1
  [p_flt,f_flt]=psd(data_filled_flt,window,fsamp,window,'linear');
  pout_flt(:,comp) = p_flt;
  semilogx(f_flt.*1000,10*log10(2*p_flt/fsamp),'r-')
  i=find(f_flt~=0 & f_flt<=1/1296000);
  fs=f_flt(i(2:length(i))); ps=p_flt(i(2:length(i)));
  lfs=log10(fs); lps=10*log10(2*ps/fsamp);
  zpoly=polyfit(lfs,lps,1);
  slope_flt=zpoly(1)./10;
  zcomp=polyval(zpoly,lfs);
  semilogx(fs*1000,zcomp,'-k');
  res=lps-zcomp;
  chi2 = sum(res.^2)/(length(res)-2);
  h=[lfs ones(length(lfs),1)];
  cov=chi2.*inv(h'*h);
  slope_sig_flt=sqrt(cov(1,1))./10;
end

xlabel('frequency (mHz)')
ylabel('power (dB)')
if comp <3
  pmin = -25;
  pmax = 30;
else
  pmin = -15;
  pmax = 40;
end
fmin = 0.000006;
fmax = 0.01;
axis([fmin,fmax,pmin,pmax])

for r=1:length(fref)
   semilogx([fref(r) fref(r)].*1000,[pmax-5 pmax-3],'k-')
   text(fref(r).*1000,pmax-2,textref(r,:))
end
txt=sprintf('%s    %s     nwin =%3.0f', site, neu_txt, nwin);
text(fmin*1.5,pmin+7,txt)
txt=sprintf('slope(unf) = %6.2f +- %6.2f mm/yr', slope_unf, slope_sig_unf);
text(fmin*1.5,pmin+4,txt);
if flt_flag == 1
  txt=sprintf('slope(flt) = %6.2f +- %6.2f mm/yr', slope_flt, slope_sig_flt);
  text(fmin*1.5,pmin+2,txt);
end

eval(['print -dpsc ' site '_' neu_flag '_psd.ps'])

end          %end dir comp loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write output freq/power files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output units are hz and db
if output_flag == 1
  fpw_output = [f pout];
  fpw_file = [site '_psd_unf.dat'];
  eval(['save ' fpw_file ' fpw_output -ascii'])
  if flt_flag == 1
    fpw_flt_output = [f pout_flt];
    fpw_flt_file = [site '_psd_flt.dat'];
    eval(['save ' fpw_flt_file ' fpw_flt_output -ascii'])
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end            %end site loop
