% called by run_coord_model.m to make plots

if run_type == 'clean',

figure
subplot(3,1,1)
plot(tneusig_c(:,1),tneusig_c(:,2),'k.')
hold on
%plot(tneusig_e(:,1),tneusig_e(:,2),'go')
plot(tneusig_e(:,1),tneusig_e(:,2),'g.')
title([site ',' blanks(2) 'north'])
ylabel 'position (m)'
subplot(3,1,2)
plot(tneusig_c(:,1),tneusig_c(:,3),'k.')
hold on
%plot(tneusig_e(:,1),tneusig_e(:,3),'go')
plot(tneusig_e(:,1),tneusig_e(:,3),'g.')
title([site ',' blanks(2) 'east'])
ylabel 'position (m)'
subplot(3,1,3)
plot(tneusig_c(:,1),tneusig_c(:,4),'k.')
hold on
%plot(tneusig_e(:,1),tneusig_e(:,4),'go')
plot(tneusig_e(:,1),tneusig_e(:,4),'g.')
title([site ',' blanks(2) 'up'])
xlabel 'time (years)'
ylabel 'position (m)'
eval(['print -dpsc ' site  'c.ps'])
close all

elseif run_type == 'model',

vrange1 = max(data)-min(data);
vmin1 = min(data) - vrange1.*0.12*(neqs);
vmax1 = max(data) + vrange1.*0.12*(nvb+1);
vrange1 = vmax1-vmin1;
vrange2 = max(data_norate_nojumps)-min(data_norate_nojumps);
vmin2 = min(data_norate_nojumps) - vrange2.*0.12*(ntau+1);
vmax2 = max(data_norate_nojumps) + vrange2.*0.12*2;
vrange2 = vmax2-vmin2;
vrange3 = max(resids(:,comp))-min(resids(:,comp));
vmin3 = min(resids(:,comp)) - vrange3.*0.12*2;
vmax3 = max(resids(:,comp)) + vrange3.*0.12*2;
vrange3 = vmax3-vmin3;
hrange = max(time)-min(time);
hmin = min(time) - hrange/40.0;
hmax = max(time) + hrange/40.0;
vshift = 0.065;
hshift = 0.2;

figure
orient tall
subplot(3,1,1)
plot(time,data,'.k')
hold on
plot(time_m,model,'r')
title(['Data with Model:' blanks(3) site ',' blanks(2) neu_txt])
kk=1;
for t=1:(1+nvb)
  txt=sprintf('slope = %6.2f +- %6.2f mm/yr (%8.4f - %8.4f)', slope(t,1), slope(t,2),VBREAKs_out(t), VBREAKs_out(t+1));
  text(hmin+hshift*hrange,vmax1-kk*vshift*vrange1,txt)
  kk=kk+1;
end
kk=1;
for i=1:neqs
  txt=sprintf('coseismic offset(%8.4f) = %6.2f +- %6.2f mm',jumps_eq(i), offset_eq(i,1), offset_eq(i,2));
  text(hmin+hshift*hrange,vmin1+((neqs+1)-kk)*vshift*vrange1,txt)
  kk=kk+1;
end
ylabel 'position (mm)'
axis([hmin hmax vmin1 vmax1])
subplot(3,1,2)
plot(time,data_norate_nojumps,'.k')
hold on
plot(time_m,model_norate_nojumps,'r')
for i=1:njumps
   line([jumps(i),jumps(i)],[-200,200])
end
kk = 1;
for i=1:ntau
   txt=sprintf('postseismic decay(%8.4f) = %6.2f +- %6.2f mm, tau = %d days',epoch_tau(i),decay(i,1), decay(i,2), TAUs(i));
   text(hmin+hshift*hrange,vmin2+((ntau+1)-kk)*vshift*vrange2,txt)
   kk = kk+1;
end
if ann_flag
  txt=sprintf('annual = %6.2f +- %6.2f mm, phase = %6.2f',annual(1,1), annual(1,2), phase);
  text(hmin+hshift*hrange,vmax2-vshift*vrange2,txt)
end
if semi_flag
  txt=sprintf('semi = %6.2f +- %6.2f mm, phase = %6.2f',semi(1,1), semi(1,2), phase_semi);
  text(hmin+hshift*hrange,vmax2-2*vshift*vrange2,txt)
end
title(['Data (trend and jumps removed), with Model:' blanks(3) site ',' blanks(2) neu_txt])
ylabel 'position (mm)'
axis([hmin hmax vmin2 vmax2])
subplot(3,1,3)
plot(time,resids(:,comp),'.k')
hold on
for i=1:njumps
   line([jumps(i),jumps(i)],[-200,200])
end
title(['Data - Model:' blanks(3) site ',' blanks(2) neu_txt])
xlabel 'time (years)'
ylabel 'residual (mm)'
%txt=sprintf('wrms = %6.2f mm',wrms_pf(comp));
txt=sprintf('rms = %6.2f mm',rms_pf(comp));
text(hmin+hshift*hrange,vmax3-vshift*vrange3,txt)
axis([hmin hmax vmin3 vmax3])
eval(['print -dpsc ' site '_' neu_flag '_' flt_flag '.ps'])
%eval(['print -djpeg ' site '_' neu_flag '_' flt_flag '.jpg'])
close all

end
