function DayPlot(x,y,c);	%Instead of plot(x,y,c);
%   Plots data grouped by day, or rather, by integral values of x. This grouping does not
%plot a line joining each individual (x,y) point as in the crude plot, which with many points
%closely spaced in the x-direction merely results in a smear from the minimum to maximum
%y-values at each resolved x-position, instead it uses summary statistics for each day's
%collection of y-values and draws lines correspondingly.
%         plot(x,y,'c');                  might have been what you used for a crude plot,
%      DayPlot(x,y,'c');                  will produce a more stylish plot, or better,
%      DayPlot(x,y,'c1,c2,c3,c4,c5,c6');  using different colours.
%   For each day's worth of data (identified by changes in the integer value of x),
%at the position corresponding to the middle of the x-variable's integer value,
%lines will be drawn with colours according to this schedule:
%      c1 for a vertical line from the minimum to the lower quartile,
%      c2 for a short horizontal dash marking the median value,
%      c3 for a vertical line from the upper quartile to the maximum,
%      c4 for a line joining the successive daily averages,
%      c5 for a stepwise line showing the successive weekly averages,
%      c6 for a horizantal line showing the overall average, and possibly, stepwise annual averages.
%   If you do not supply six separate entries (that is, a text string with five commas)
%then the necessary collection will be developed by repeating the supplied text
%a sufficient number of times. If you do not want a particular sequence to be shown,
%supply a blank for the unwanted sequence's entry.
%   Note that each day's average is the average of the "good" values only.
%The week's average is the average of the contained good daily averages,
%and is not necessarily the same as the average of all the individual values.
%Similarly, the overall average is the average of the weekly averages,
%not of the contained values. These averages would be the same only if
%every day had the same number of values and every week had the same number of days.
%   Matlab's date numbering scheme is not usefully aligned. 1/1/2001, a Monday
%is given daynumber 730852 and the remainder modulo 7 is 3, not 1. Thus the usage
%of WeekShift to arrange that a week runs from Sunday to Saturday, as per Genesis.
%   Concocted by R.N.McLean (whom God preserve), 29'th June MMIX.
WeekShift = mod(datenum(2001,1,0),7);	%1/1/2001 is a Monday, so 0/1/2001 is a Sunday.
ccc = c;			%Start by inspecting the colour code collection.
if size(ccc,2) < 1		%If we have been given nothing,
 ccc = 'g,r,m,c,b,k';		%This will do.
end;				%And prevent flops.
while size(findstr(ccc,','),2) < 5	%Enough yet?
 ccc = strcat(ccc,',',c);		%No, concoct more.
end;				%Perhaps not one at a time.
cc = textscan(ccc,'%s%s%s%s%s%s','Delimiter',',');	%I require six plot instruction codes.
%disp([char(cc{1}),',',char(cc{2}),',',char(cc{3}),',',char(cc{4}),',',char(cc{5}),',',char(cc{6})]);

ug = find(~isnan(y),1,'last');	%Disregard trailing non-values.
if size(ug,1) == 0, disp('All y-values bad!'); return; end;	%Surely not.
hop = find(fix(x(2:end)) - fix(x(1:end - 1)));	%Indices of the last element before a step in the day number.
if size(hop,1) <= 0	%Were there any day number changes?
  nh = 1;			%No, so there is only one segment.
  clear hop; hop(1) = ug;	%And this fingers the last value worth plotting.
 else			%Otherwise, there is at least one day step, with last values of the day fingered..
  nh = find(hop < ug,1,'last');	%The last day has no step to the next day so its end is unfingered.
  if size(nh,1) == 0		%Was there a step before the last good value? Which may not be in the last day.
    clear nh; nh = 1		%No, so one segment only.
   else				%But if there are hops, index ug is somewhere in the last day.
    nh = nh + 1;		%So add another segment ender, as the last hop marks the end of the penultimate day.
  end;			%So much for the hop count.
  hop(nh) = ug;		%End the last segment.
end;			%That was painful.
YT = 0; Yn = 0; Y = -Inf; lY = 0;	%Set up for some annual averages.
nd = 0; ld = fix(x(end)) - fix(x(1)) + 1;	%Number of days in the time span; there may be gaps in the data.
%disp(['nhop=',int2str(nh),',hop(nh)=',int2str(hop(nh)),',ug=',int2str(ug),',dayspan=',int2str(ld)]);
Avx(1:ld) = 0; Avy(1:ld) = 0;	%Some days may prove to have no data, though.
nb = 0; lb = ld*3;		%Number of bar pieces, supposing all days have data.
X(1:lb) = 0; QL(1:lb) = 0; QM(1:lb) = 0; QU(1:lb) = 0;	%Prepare some arrays, lb known.
ws = 1; lw = ceil(ld/7) + 1;	%Week step count. Add one for the tail end.
Wxx(1:lw) = 0; Wav(1:lw) = 0;	%Reserve some space.
Wsum = 0; WsumN = 0;		%The starting week awaits day averages.
Av = 0; AvN = 0;		%And the overall average awaits week averages.
xx = fix(x(1)) - 0.5;		%The middle of the previous day.
l1 = 1;				%Now scan the day chunks.
for i = 1:nh			%One by one.
 l2 = hop(i);			%The last of this day's chunk.
 if l2 < l1			%This shouldn't happen if my inspections worked...
  disp(['Hic! l1,l2=',int2str(l1),',',int2str(l2),',i=',int2str(i),', nh=',int2str(nh)]);	%But alas,
 end;				%There may be some confusion, especially over "bad" values.
 clear live; live = find(~isnan(y(l1:l2))) + l1 - 1;	%Some values might be bad.
 nn = size(live,1);		%The number of good values in this day.
 if nn > 0			%Surely positive. Unless data are missing.
   nd = nd + 1;			%Another day is complete.
   pxx = xx;			%Recall where its predecessor was.
   xx = (x(l1) + x(l2))/2;	%Its centre.
   if xx - pxx > 3		%An excessive lurch?
     Avx(nd) = NaN;		%Yes. Break the chain.
     Avy(nd) = NaN;		%So there will be a gap in the plot.
     nd = nd + 1;		%Onwards.
   end;				%So perhaps not all day points will be used.
   Avx(nd) = xx;		%Save for the plot.
   clear snarl; snarl = sort(y(live));	%MatLab doesn't know about quartiles, nor Medianof(up to median)).
   Avy(nd) = sum(snarl)/nn;		%Calculate the average.
   clear idiot; idiot = datevec(x(l1)); Year = idiot(1);	%The year of the first value of the day.
   if Year > Y			%A new year yet?
     if Yn >=  222		%A properly populous year?
       lY = lY + 1;		%Yes. Another one.
       Yav(lY) = YT/Yn;		%Its average of its numerous live days.
       YY(lY) = Y;		%I think there may be gaps.
     end;			%So much for the old year.
     Y = Year;			%This is the year it is.
     YT = 0; Yn = 0;		%Refresh the counts for the new year.
   end;				%So much for year crossings.
   YT = YT + Avy(nd); Yn = Yn + 1;	%Count this day into its year average.
   tw = fix(fix(x(l2) - WeekShift)/7);	%This day's week number.
   if WsumN <= 0		%This catches the start of the first week with live data.
     Wxx(ws) = fix(x(l1));	%The start of its first day.
     wn = tw;			%This is the week we're in.
   end;				%And I don't know when it finishes.
   if wn ~= tw			%Starting a new week number?
    Wav(ws) = Wsum/WsumN;		%Yes. Save the average of the day averages.
   %disp(['Wav(ws=',int2str(ws),')=',num2str(Wav(ws)),',WsumN=',int2str(WsumN)]);
    Av = Av + Wav(ws); AvN = AvN + 1;	%And augment the overall average.
    Wsum = 0; WsumN = 0;		%And prepare for the next.
    ws = ws + 1;			%Another week completed. Wxx(ws) was set at its start.
    if tw - wn > 1		%A lurch?
      Wxx(ws) = fix(pxx) + 0.9999;	%Yes. The end of the last day that had been admitted.
      Wav(ws) = Wav(ws - 1);		%Thus form the level part, not extending into the absent data..
      ws = ws + 1;%			%Advance another.
      Wxx(ws) = fix(pxx) + 1;		%Start of the next day.
      Wav(ws) = NaN;			%Break the chain.
      ws = ws + 1;			%Get ready to start the resumed step.
    end;			%Enough of gaps. Perhaps extra week points will be needed.
    Wxx(ws) = fix(x(l1));	%The week boundary. (Time zero of the incoming day)
    wn = tw;			%Recall the week number to spot changes.
   end;				%So much for a new week. Wav(ws) will be placed later.
   Wsum = Wsum + Avy(nd);	%Augment the sum of the daily averages.
   WsumN = WsumN + 1;		%Counting them.
   nb = nb + 1;			%Now prepare the vertical lines.
   X(nb:nb + 2) = xx;		%The x-position for the bars.
   QL(nb) = snarl(1);		%The minimum.
   if nn <= 1			%Trouble with insufficient points?
     QL(nb + 1) = QL(nb);	%Yes. One point means no second point to interpolate with.
     QM(nb:nb + 1) = QL(nb);	%So just make everything the same.
     QU(nb:nb + 1) = QL(nb);	%Splotching atop each other.
    else			%But ordinarily, there will be many points.
     q = (nn - 1)/4;		%Points 1:nn, compute the offset on 0:(nn - 1) for the quartile position.
     QL(nb + 1) = (fix(q + 1) - q)*snarl(1 + fix(q)) + (q - fix(q))*snarl(1 + fix(q + 1));	%Interpolate.
     q = (nn - 1)/2;		%Find the median.
     QM(nb) = (fix(q + 1) - q)*snarl(1 + fix(q)) + (q - fix(q))*snarl(1 + fix(q + 1));		%Interpolate.
     QM(nb + 1) = QM(nb);		%This will have a horizontal bar.
     QU(nb) = snarl(nn);		%The maximum.
     q = (nn - 1)*3/4;		%Interpolation via weighted average based on the distance from adjacent indices.
     QU(nb + 1) = (fix(q + 1) - q)*snarl(1 + fix(q)) + (q - fix(q))*snarl(1 + fix(q + 1));	%The right way around, I hope.
   end;				%So much for those strokes.
   nb = nb + 2;			%Advance the finger, knowing another is to come.
   QL(nb) = NaN; QM(nb) = NaN; QU(nb) = NaN;	%Split one bar from the next.
 end;				%So much for that day.
 l1 = l2 + 1;			%Advance.
end;			%On to the next day.
if Yn >=  222			%Tail end. Some year has finished.
 lY = lY + 1;			%And if populous enough,
 Yav(lY) = YT/Yn;		%Record its average of its days.
 YY(lY) = Y;			%Not forgetting what year it is.
end;			%Righto, enough calculation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Plotting commences.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on;			%In multiple assaults.
expose = 'on';			%No plot action has been exposed to the Legend List.
if size(char(cc{1}),2) > 0	%Show the min to lower quartile?
 line(X(1:nb),QL(1:nb),'Color',char(cc{1}),'HandleVisibility',expose);	%Yes.
 expose = 'off';		%It has been sent.
end;
if size(char(cc{3}),2) > 0	%Show the upper quartile to max?
 line(X(1:nb),QU(1:nb),'Color',char(cc{3}),'HandleVisibility',expose);	%Yes.
 expose = 'off';		%The first appearance will do for the legend.
end;
if size(char(cc{2}),2) > 0	%Show the median values?
 X(1:2:nb) = X(1:2:nb) - 0.25;	%Yes. These have horizontal lines.
 X(2:2:nb) = X(2:2:nb) + 0.25;	%So, corrupt the X-array (last usage). Half a day long.
 line(X(1:nb),QM(1:nb),'Color',char(cc{2}),'HandleVisibility',expose);
 expose = 'off';		%As I have chosen the colours to match.
end;

if size(char(cc{4}),2) > 0	%Show the daily averages?
 plot(Avx(1:nd),Avy(1:nd),char(cc{4}),'HandleVisibility',expose);
 expose = 'off';		%Each component having a similar colour.
end;

Wav(ws) = Wsum/WsumN;		%Late placement for the last average's beginning.
Av = (Av + Wav(ws))/(AvN + 1);	%Thus an overall average, before the tail end is attached..
ws = ws + 1;			%Advance to mark its end, at the last live day.
Wxx(ws) = fix(xx) + 0.9999;	%Add 1, and damn matlab regards a new day as starting.
Wav(ws) = Wav(ws - 1);		%Thus level for the time interval.
if size(char(cc{5}),2) > 0	%Show the weekly averages?
 stairs(Wxx(1:ws),Wav(1:ws),char(cc{5}),'HandleVisibility',expose);	%I hope this is aligned correctly.
 expose = 'off';		%Any one component colour will have to do for the legend.
end;

if size(char(cc{6}),2) > 0 && ws > 6	%Show the overall average and annual averages?
  if lY > 3			%Yes. Worth showing annual averages?
    ys = 0;			%There seems to be enough of them.
    pY = YY(1) - 1;		%Syncopation to spot year lurches.
    for i = 1:lY		%Step through the collection of annual averages.
      ys = ys + 1;		%Another step starts.
      if YY(i) - pY > 1		%A lurch in the years?
        Yy(ys) = Yy(ys - 1);		%Yes. End the previous step with its value.
        Yx(ys) = datenum(pY,12,31,59,59,59.99);	%At the end of its year.
        ys = ys + 1;		%Prepare for a new step.
        Yy(ys) = NaN;		%But with a break.
        Yx(ys) = datenum(pY,1,1,0,0,0);	%Starting with the next year.
        ys = ys + 1;		%Advance for the standard step placement to come.
      end;			%Enough of splits.
      Yy(ys) = Yav(i);		%Place the annual average.
      pY = YY(i);		%Remember this year's number.
      Yx(ys) = datenum(pY,1,1,0,0,.999);	%Avoid encouraging matlab's span extensuion.
    end;			%On to the next annual average start time.
    ys = ys + 1;		%Now end the last annual average.
    Yy(ys) = Yy(ys - 1);	%Its value lasts all through its year.
    Yx(ys) = datenum(YY(lY),12,31,23,59,59.999);	%End time, again discouraging matlab's extensions.
    Yx(ys) = min(Yx(ys),x(end));%Likewise, the year may end after the last datum.
    Yx(1)  = max(Yx(1), x(1));	%And start before. So don't extend the plot area.
    stairs(Yx(1:ys),Yy(1:ys),char(cc{6}),'HandleVisibility',expose);	%All that struggle to suit this.
    expose = 'off';		%It might be easier to just plot the points.
  end;				%Now for the overall average..
 plot([x(1) Wxx(1) Wxx(ws) x(end)], [NaN Av Av NaN],char(cc{6}),'HandleVisibility',expose);	%Use the span of the weekly averages.
 expose = 'off';		%Might as well be complete.
end;				%That was a struggle.

set(gca,'TickDir','out');	%Do not deface the data area.
set(gca,'TickLength',get(gca,'TickLength')/4);	%And be quiet about it.
%END. Just fall off.


