clear all

data = csvread('sounding_dec2013_08-21.csv');

pot_snd = data(:,9:11:end);
hgt_snd = data(:,2:11:end);
dir_snd = data(:,7:11:end);
spd_snd = data(:,8:11:end);

temp   = importdata('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/dec2013/time_series/SLC.d03.TH',' ',1);
temp_k = temp.data(:,2:end);
hgt    = importdata('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/dec2013/time_series/SLC.d03.PH',' ',1);
hgt_m  = hgt.data(:,2:end);
uws    = importdata('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/dec2013/time_series/SLC.d03.UU',' ',1);
uws_ms = uws.data(:,2:end);
vws    = importdata('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/dec2013/time_series/SLC.d03.VV',' ',1);
vws_ms = vws.data(:,2:end);

wind_abs = (uws_ms.^2 + vws_ms.^2).^(1/2);
wind_dir_trig_to = atan2(uws_ms./wind_abs, vws_ms./wind_abs);
wind_dir_trig_to_degrees = wind_dir_trig_to .* (180/pi);

wind_dir_trig_from_degrees = wind_dir_trig_to_degrees + 180;

wind_dir_cardinal = 90 - wind_dir_trig_from_degrees;

%%

figure('Position',[100, 100, 1200, 600])

subplot(1,2,1)
plot(pot_snd(:,1),hgt_snd(:,1),'k-','LineWidth',2)
hold on
plot(temp_k(1,:),hgt_m(1,:),'k--','LineWidth',2)
axis([280 290 hgt_snd(1,1) 3000])
xlabel('Potential Temperature (K)')
ylabel('Height A.S.L. (m)')
legend('Station','Model','Location','northwest')

subplot(1,2,2)
plot(dir_snd(:,1),hgt_snd(:,1),'k.','MarkerSize',20)
hold on
plot(wind_dir_trig_from_degrees(1,:),hgt_m(1,:),'k*','MarkerSize',10)
legend('Station','Model','Location','northwest')

ax1 = gca;
set(ax1,'box','off')
xlabel('Wind Direction (^o)');
axis(ax1,[0 360 hgt_snd(1,1) 3000])

ax1_pos = ax1.Position; % position of first axes
ax2 = axes('Position',ax1_pos,...
    'XAxisLocation','top',...
    'YAxisLocation','right',...
    'Color','none');
hold on

ax2.XColor = 'b';
ax2.YColor = 'b';
xlabel('Wind Speed');

plot(spd_snd(:,1),hgt_snd(:,1),'Parent',ax2,'Color','b','LineWidth',2)
hold on
plot(wind_abs(1,:),hgt_m(1,:),'b--','LineWidth',2)
axis(ax2,[0 25 hgt_snd(1,1) 3000])


%%

ds = 8640;

day =  {'08' '08' '09' '09' '10' '10' '11' '11' '12' '12' '13' '13' '14' '14' '15' '15' '16' '16' '17' '17' '18' '18' '19' '19' '20' '20' '21'};
hour = {'00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00' '12' '00'};

for i = 1:26
   figure('Position',[100, 100, 1200, 600])

   subplot(1,2,1)
   plot(pot_snd(:,i),hgt_snd(:,i),'k-','LineWidth',2)
   hold on
   plot(temp_k(i*ds,:),hgt_m(i*ds,:),'k--','LineWidth',2)
   axis([260 300 hgt_snd(1,i) 3000])
   xlabel('Potential Temperature (K)')
   ylabel('Height A.S.L. (m)')
   legend('Station','Model','Location','northwest')

   subplot(1,2,2)
   plot(dir_snd(:,i),hgt_snd(:,i),'k.','MarkerSize',20)
   hold on
   plot(wind_dir_trig_from_degrees(i*ds,:),hgt_m(i*ds,:),'k*','MarkerSize',10)
   legend('Station','Model','Location','northwest')

   ax1 = gca;
   set(ax1,'box','off')
   xlabel('Wind Direction (^o)');
   axis(ax1,[0 360 hgt_snd(1,i) 3000])

   ax1_pos = ax1.Position; % position of first axes
   ax2 = axes('Position',ax1_pos,...
       'XAxisLocation','top',...
       'YAxisLocation','right',...
       'Color','none');
   hold on

   ax2.XColor = 'b';
   ax2.YColor = 'b';
   xlabel('Wind Speed');

   plot(spd_snd(:,i),hgt_snd(:,i),'Parent',ax2,'Color','b','LineWidth',2)
   hold on
   plot(wind_abs(i*ds,:),hgt_m(i*ds,:),'b--','LineWidth',2)
   axis(ax2,[0 25 hgt_snd(1,i) 3000])
   
   set(gcf,'PaperUnits','inches','PaperPosition',[0 0 12 6])
   saveas(gcf,strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/daq_runs/dec2013/images/soundings/sounding_',day{i+1},hour{i+1},'.png'))

   close all
end