% clamp_public_plot.m
% time series plotting program for clamp aws data (4 per page)
% based on aws_model_plot.m - Ian Renfrew Jan 2000/March 2001

load C:/public/data/clamp/matlab_data/clamp_aws.mat

%set up time period 
sdmin=datenum('01-may-1999 00:00:00'); 
sdmax=datenum('14-may-1999 24:00:00');

sertick = round((sdmax-sdmin)/7);
e1 = find((sd > sdmin)&(sd < sdmax)&(loc==1));
e2 = find((sd > sdmin)&(sd < sdmax)&(loc==2));
e3 = find((sd > sdmin)&(sd < sdmax)&(loc==3));
e4 = find((sd > sdmin)&(sd < sdmax)&(loc==4));
e5 = find((sd > sdmin)&(sd < sdmax)&(loc==5));

% set up some plotting arrays: 
% locations (clamp 1 to 4); plus 5 = halley,
stations = char('C1','C2','C3','C4','Halley');
cl_colour = char('c-','b-','r--','g-.','k-');

% variables: sd; all variables are 1 by n vectors, where loc gives location
% i=1 mslp; i=2 ws; i=3 wd; i=4 t0; i=5 t1; i=6 q0; i=7 q1; i=8 rh0; i=9 rh1
% i=10 rhi0; i=11 rhi1; i=12 theta0; i=13 theta1; i=14 p; i=15 t1-t0;
% i=16 (t1-t0)*v i=17 u; i=18 v; i=19 accumulation; i=20 voltage;
% i=21 cloudh; i=22 ppert; i=23 h_c3_ppert; i=24 difference theta1; 

%--------------------------------------------------------------
% loop over panels per page, where variables plotted via i subscript
% eg. 2 3 22 21; 13 15 11 7;  
for j=1:4
 if (j==1) i=2; end;
 if (j==2) i=3; end;
 if (j==3) i=22; end;
 if (j==4) i=21; end;

%reset plotting variables
 c1=NaN*ones(size(e1)); c2=NaN*ones(size(e2)); c3=NaN*ones(size(e3)); 
 c4=NaN*ones(size(e4)); h=NaN*ones(size(e5));

 if (i==1)
  c1=mslp(e1); c2=mslp(e2); c3=mslp(e3); c4=mslp(e4); h=mslp(e5);
  maxax=1020; minax =960;
  yl='sea level pressure (mb)'
 elseif (i==2)
  c1=ws(e1); c2=ws(e2); c3=ws(e3); c4=ws(e4); h=ws(e5);
  maxax=30; minax=0;
  yl='wind speed (m s^{-1})'
 elseif (i==3)
  c1=wd(e1); c2=wd(e2); c3=wd(e3); c4=wd(e4); h=wd(e5);
  maxax=360; minax=0;
  yl='wind direction'
 elseif (i==4)
  c1=t0(e1); c2=t0(e2); c3=t0(e3); c4=t0(e4);
  maxax=0; minax=-40;
  yl='t0 (^o C)' 
 elseif (i==5)
  c1=t1(e1); c2=t1(e2); c3=t1(e3); c4=t1(e4); h=t1(e5);
  maxax=0; minax=-40;
  yl='t_a (^o C)' 
 elseif (i==6)
  c1=q0(e1); c2=q0(e2); c3=q0(e3); c4=q0(e4); 
  maxax=2; minax=0;
  yl='q0 (g kg^{-1})'
 elseif (i==7)
  c1=q1(e1); c2=q1(e2); c3=q1(e3); c4=q1(e4); h=q1(e5);
  maxax=2; minax=0;
  yl='specific humidity (g kg^{-1})'
 elseif (i==8)
  c1=rh0(e1); c2=rh0(e2); c3=rh0(e3); c4=rh0(e4);
  maxax=110; minax=0;
  yl='relative humidity [0]'
 elseif (i==9)
  c1=rh1(e1); c2=rh1(e2); c3=rh1(e3); c4=rh1(e4); h=rh1(e5); 
  maxax=110; minax=0;
  yl='RH'
 elseif (i==10)
  c1=rhi0(e1);  c2=rhi0(e2); c3=rhi0(e3); c4=rhi0(e4); 
  maxax=110; minax=0;
  yl='relative humidity w.r.t ice [0]'
 elseif (i==11)
  c1=rhi1(e1);  c2=rhi1(e2); c3=rhi1(e3); c4=rhi1(e4); h=rhi1(e5);
  maxax=110; minax=0;
  yl='RH_{ice}'
 elseif (i==12)
  c1=theta0(e1);  c2=theta0(e2); c3=theta0(e3); c4=theta0(e4); 
  maxax=280; minax=240;
  yl='potential temperature [0] (K)'
 elseif (i==13)
  c1=theta1(e1);  c2=theta1(e2); c3=theta1(e3); c4=theta1(e4); h=theta1(e5);
  maxax=287; minax=237;
  yl='\theta (K)'
 elseif (i==14)
  c1=p(e1);  c2=p(e2); c3=p(e3); c4=p(e4); h=p(e5); 
  maxax=1050; minax=750;
  yl='station pressure (mb)'
 elseif (i==15)
  c1=(theta1(e1)-theta0(e1))/1.5; c2=(theta1(e2)-theta0(e2))/1.5; 
  c3=(theta1(e3)-theta0(e3))/1.5; c4=(theta1(e4)-theta0(e4))/1.5; 
  h=(theta1(e5)-theta0(e5))/2.0; 
  maxax=4.5; minax=-4.5;
  yl=' \Delta \theta / \Delta z (K m^{-1})'
 elseif (i==16)
  c1=(t1(e1)-t0(e1)).*v(e1);  c2=(t1(e2)-t0(e2)).*v(e2);
  c3=(t1(e3)-t0(e3)).*v(e3);  c4=(t1(e4)-t0(e4)).*v(e4); 
  maxax=20; minax=-10;
  yl='(t1 - t0) x v'
 elseif (i==17)
  c1=u(e1); c2=u(e2); c3=u(e3); c4=u(e4); h=u(e5); 
  maxax=15; minax=-15;
  yl='u (m s^{-1})'
 elseif (i==18)
  c1=v(e1); c2=v(e2); c3=v(e3); c4=v(e4); h=v(e5); 
  maxax=15; minax=-15;
  yl='v (m s^{-1})'
 elseif (i==19)
  disp('no accumulation data')
 elseif (i==20)
  c1=volt(e1); c2=volt(e2); c3=volt(e3); c4=volt(e4);
  maxax=16; minax=5;
  yl='battery voltage'
 elseif (i==21)
  h1=cloud(e5);
  sdcloud1 = ~isnan(h1).*sd(e5);
  sdcloud = sdcloud1(find(sdcloud1 >1));
  h = filternan(h1);
  maxax=8.5; minax=-0.5;
  yl='total cloud'
 elseif (i==22)
  c1=ppert(e1);  c2=ppert(e2); c3=ppert(e3); c4=ppert(e4); 
  h=ppert(e5); 
  maxax=30; minax=-30;
  yl='perturbation pressure (mb)'
 elseif (i==23)
  h=h_c3_ppert; 
  maxax=10; minax=-10;
  yl='difference pert. p (mb)'
elseif (i==24)
  c2=c2_c4_theta1; c3=c3_c4_theta1;
  maxax=10; minax=-10;
  yl='difference theta (K)'
 end;

%----------------------------------------------------------------------------
% plot the time series
%----------------------------------------------------------------------------
 if (i==21) 
  subplot(4,1,j), plot(sdcloud, h, cl_colour(5,:),'linewidth',1.5);
 else
  subplot(4,1,j), plot(sd(e5), h, cl_colour(5,:),'linewidth',1.5);
 end;
hold on;
 h=gca;
 subplot(4,1,j), plot(sd(e1), c1, cl_colour(1,:));
 subplot(4,1,j), plot(sd(e2), c2, cl_colour(2,:));
 subplot(4,1,j), plot(sd(e3), c3, cl_colour(3,:));
 subplot(4,1,j), plot(sd(e4), c4, cl_colour(4,:));
 if (i==21) subplot(4,1,j), plot(sdcloud, h, cl_colour(5,:)); end;
 set(h, 'fontsize', 11);
 set(h, 'xtick',[sdmin:sertick:sdmax])
 datetick('x','mm/dd')
 ylabel(yl);
 axis([sdmin sdmax minax maxax]); 
 if (i==3) set(h, 'ydir', 'reverse'); end;
 orient tall

% add legend
 legend_plot=1;
 if (j == 1) 
  if (legend_plot == 1)
   dy=(maxax-minax)/12;
   dx=(sdmax-sdmin)/14;
   for i=2:4
    subplot(4,1,1), plot([sdmin+0.5*dx sdmin+2*dx],...
			[maxax-(i-1)*dy maxax-(i-1)*dy], cl_colour(i,:));
    tt=text(sdmin+2.5*dx, maxax-(i-1)*dy, stations(i,:));
   end;
   i=5;
   subplot(4,1,1), plot([sdmin+0.5*dx sdmin+2*dx], ...
	[maxax-(i-1)*dy maxax-(i-1)*dy], cl_colour(i,:),'linewidth',1.5);
   tt=text(sdmin+2.5*dx, maxax-(i-1)*dy, stations(i,:));
  end; 
 end;

end; %j loop

