clear all;close all;
%load ab1202_db
%load rv1203_db
% extract the data that you want from the structure
%  this routine assumes that your structure is in a very specific format
%  could be made more general to search btl.names to find the appropriate
%  column for the data you want.
station = [];
theta_primary = [];
t_primary = [];
theta_secondary = [];
t_secondary = [];
pr = [];

for II = 1: length(btl)
    station = [station; II*ones(length(btl(II).data(:,3)),1)];
    theta_primary = [theta_primary; btl(II).data(:,3)];
    t_primary = [t_primary; btl(II).data(:,13)];
    theta_secondary = [theta_secondary; btl(II).data(:,4)];
    t_secondary = [t_secondary; btl(II).data(:,14)];
    pr = [pr; btl(II).data(:,12)];
end

% 
% extract only the data below 500 dbar
%
II = find(pr >= 200.0);
deep_station = station(II);
deep_theta_primary = theta_primary(II);
deep_t_primary = t_primary(II);
deep_theta_secondary = theta_secondary(II);
deep_t_secondary = t_secondary(II);
deep_pr = pr(II);

mean_theta_diff = mean(deep_theta_primary-deep_theta_secondary);
stddev_theta_diff = std(deep_theta_primary-deep_theta_secondary);
median_theta_diff = median(deep_theta_primary-deep_theta_secondary);
median_t_diff = median(deep_t_primary-deep_t_secondary);
stddev_t_diff = std(deep_t_primary-deep_t_secondary);
iqr_theta_diff = iqr(deep_theta_primary-deep_theta_secondary)/1.349;
var_median_theta_diff = var_median(deep_theta_primary-deep_theta_secondary)/1.349;

%KK = find (deep_station >= 29 & deep_station <= 33);

%
% Now make the figure that you want
%

figure;
  plot(station, theta_primary-theta_secondary,'+')
  
figure;
  subplot(2,1,1)
  plot(deep_station, deep_theta_primary-deep_theta_secondary,'+')
  h=gca;
  set(h,'ylim',[-0.0032 0.0032])
  xlimits = get(h,'xlim');
  xlabel('Station number','fonts',14)
  ylabel('Primary-Secondary','fonts',14)
  title('Potential Temperature (Prs > 200)','fonts',14)
  font(14)
  hold on;
  plot(xlimits, [median_theta_diff median_theta_diff],'r')
  plot(xlimits, [median_theta_diff-iqr_theta_diff median_theta_diff-iqr_theta_diff],'r--')
  plot(xlimits, [median_theta_diff+iqr_theta_diff median_theta_diff+iqr_theta_diff],'r--')
  
  plot(xlimits, [mean_theta_diff mean_theta_diff],'g')
  plot(xlimits, [mean_theta_diff-stddev_theta_diff mean_theta_diff-stddev_theta_diff],'g--')
  plot(xlimits, [mean_theta_diff+stddev_theta_diff mean_theta_diff+stddev_theta_diff],'g--')
  
  

  %plot(deep_station(KK), deep_theta_primary(KK)-deep_theta_secondary(KK),'k+')
  
  subplot(2,1,2)
  plot(deep_pr, deep_theta_primary-deep_theta_secondary,'+')
  h=gca;
  set(h,'ylim',[-0.0032 0.0032])
  xlimits = get(h,'xlim');
  xlabel('Pressure, dbar','fonts',14)
  ylabel('Primary-Secondary','fonts',14)
  font(14)
  %title('Potential Temperature')
  hold on;
  plot(xlimits, [median_theta_diff median_theta_diff],'r')
  plot(xlimits, [median_theta_diff-iqr_theta_diff median_theta_diff-iqr_theta_diff],'r--')
  plot(xlimits, [median_theta_diff+iqr_theta_diff median_theta_diff+iqr_theta_diff],'r--')
  %plot(deep_pr(KK), deep_theta_primary(KK)-deep_theta_secondary(KK),'k+')
  
  print('-depsc','figures/t_sensor_diff_gom')
  
  figure;
  plot(station,t_primary-t_secondary,'+')
  ylabel('T - Difference {\circ}C')
  xlabel('Station')
  
  
  figure;
  font(14)
  plot(station, t_primary-t_secondary,'g+')
  hold on;
  plot(deep_station,deep_t_primary-deep_t_secondary,'+')
  ylabel('T1 - T2 ({\circ}C)','fonts',16)
  xlabel('Station ','fonts',16)
  title('Temperature Differences ','fonts',16)
  hold on;
  h=gca;
  xlimits = get(h,'xlim');
  plot(xlimits, [median_t_diff median_t_diff],'r')
  plot(xlimits, [median_t_diff-stddev_t_diff median_theta_diff-stddev_t_diff],'r--')
  plot(xlimits, [median_t_diff+stddev_t_diff median_theta_diff+stddev_t_diff],'r--')
  
  
  figure;
  font(14)
  plot(pr, t_primary-t_secondary,'g+')
  hold on;
  plot(deep_pr,deep_t_primary-deep_t_secondary,'+')
  ylabel('T1 - T2 ({\circ}C)','fonts',16)
  xlabel('Pressure ','fonts',16)
  title('Temperature Differences ','fonts',16)
  hold on;
  h=gca;
  xlimits = get(h,'xlim');
  plot(xlimits, [median_t_diff median_t_diff],'r')
  plot(xlimits, [median_t_diff-stddev_t_diff median_theta_diff-stddev_t_diff],'r--')
  plot(xlimits, [median_t_diff+stddev_t_diff median_theta_diff+stddev_t_diff],'r--')