% This code creates the endogenous clock figure for the cyst-paper 
% It is a sub-section of the more extensive code: endo_plot_all.m, which
% includes plots of each individual data and more extensive analysis
% Julian days are calculated for a common year base (1995)
% to enable plotting against the same gregorian axis.

clf
clear all

addpath timeplt

% read and plot Matrai Data
[mat_date,st1,st2,st3] = textread('Matrai_NaN.txt','%s %d %d %d');
mat_date = char(mat_date);
mat_date = mat_date;

% loop through each date and isolate the day and the month based on the
% position of the /

for i = 1:size(mat_date,1);
    aa = find(mat_date(i,:) == '/');
    mon_mat(i) = str2num(mat_date(i,1:(aa(1)-1)));
    day_mat(i) = str2num(mat_date(i,(aa(1)+1):(aa(2)-1)));
    year_mat(i) = str2num(mat_date(i,(aa(2)+1:end)));
    jday_mat(i) = julian([year_mat(i),mon_mat(i),day_mat(i),12,0,0]);
end

% read 84-88 data, each of three isolates
[old_date,iso1,iso2,iso3] = textread('clock84_88.txt','%s %d %d %d');
old_date = char(old_date);

for i = 1:size(old_date,1);
    aa = find(old_date(i,:) == '/');
    mon_old(i) = str2num(old_date(i,1:(aa(1)-1)));
    day_old(i) = str2num(old_date(i,(aa(1)+1):(aa(2)-1)));
    year_old(i) = str2num(old_date(i,(aa(2)+1:end)));
    jday_old(i) = julian([year_old(i),mon_old(i),day_old(i),12,0,0]);
end

% loop through to isolate day and month based on position of /

% read in 1998 data, 1 isolate, data points 1-22
[amy_date,amy,dum] = textread('clock_1998.txt','%s %f %f');

amy_date = char(amy_date);

for i = 1:size(amy_date,1);
    aa = find(amy_date(i,:) == '/');
    mon_amy(i) = str2num(amy_date(i,1:(aa(1)-1)));
    day_amy(i) = str2num(amy_date(i,(aa(1)+1):(aa(2)-1)));
    year_amy(i) = str2num(amy_date(i,(aa(2)+1:end)));
    jday_amy(i) = julian([year_amy(i)+1900,mon_amy(i),day_amy(i),12,0,0]);
end

aa = find(isfinite(st1));
mod61(:,1) = st1(aa(1:22));
mod61(:,2) = mon_mat(aa(1:22))';
mod61(:,3) = day_mat(aa(1:22))';

aa = find(isfinite(st2));
mod94(:,1) = st2(aa(1:22));
mod94(:,2) = mon_mat(aa(1:22))';
mod94(:,3) = day_mat(aa(1:22))';

aa = find(isfinite(st3));
mod4(:,1) = st3(aa(1:19));
mod4(:,2) = mon_mat(aa(1:19))';
mod4(:,3) = day_mat(aa(1:19))';

aa = find(isfinite(iso1));
mod29_84(:,1) = iso1(aa(2:7));
mod29_84(:,2) = mon_old(aa(2:7))';
mod29_84(:,3) = day_old(aa(2:7))';

aa = find(isfinite(iso2));
mod29_85(:,1) = iso2(aa(1:11));
mod29_85(:,2) = mon_old(aa(1:11))';
mod29_85(:,3) = day_old(aa(1:11))';

aa = find(isfinite(iso3));
mod29_87(:,1) = iso3(aa(1:3));
mod29_87(:,2) = mon_old(aa(1:3))';
mod29_87(:,3) = day_old(aa(1:3))';

aa = find(isfinite(amy));
mod38(:,1) = amy(aa(2:27));
mod38(:,2) = mon_amy(aa(2:27))';
mod38(:,3) = day_amy(aa(2:27))';

% plot on uniform, 1995 axis - do conversion to year day using function
% "yearday" which considers january 1st yearday 1.
% neglect the influence of the leap-year - we are trying to map 
% onto the same 365 day axis - we need to either treat them all
% as leap-years or non - leap-years...either way there is only
% a slight inaccuracy.

bb = size(mod61,1);
for i = 1:bb
  temp = yearday([1995 mod61(i,2:3)]) + [0 21];
%  temp = yearday([1995 mod61(i,2:3)]);
  yd61(i) = temp(1,2);
end
aa = find(yd61 > 365);
yd61(aa) = yd61(aa) - 365;

clear temp;
temp = mod61;
clear mod61;
[yd61 sort_ind] = sort(yd61);
mod61 = temp(sort_ind);
clear temp;

bb = size(mod94,1);
for i = 1:bb
  temp = yearday([1995 mod94(i,2:3)]) + [0 21];
%  temp = yearday([1995 mod94(i,2:3)]);
  yd94(i) = temp(1,2);
end
aa = find(yd94 > 365);
yd94(aa) = yd94(aa) - 365;

clear temp;
temp = mod94;
clear mod94;
[yd94 sort_ind] = sort(yd94);
mod94 = temp(sort_ind);
clear temp;

bb = size(mod4,1);
for i = 1:bb
  temp = yearday([1995 mod4(i,2:3)]) + [0 21];
%  temp = yearday([1995 mod4(i,2:3)]);
  yd4(i) = temp(1,2);
end
aa = find(yd4 > 365);
yd4(aa) = yd4(aa) - 365;

clear temp;
temp = mod4;
clear mod4;
[yd4 sort_ind] = sort(yd4);
mod4 = temp(sort_ind);
clear temp;

bb = size(mod29_84,1);
for i = 1:bb
  temp = yearday([1995 mod29_84(i,2:3)]) + [0 14];
%  temp = yearday([1995 mod29_84(i,2:3)]);
  yd29_84(i) = temp(1,2);
end
aa = find(yd29_84 > 365);
yd29_84(aa) = yd29_84(aa) - 365;

clear temp;
temp = mod29_84;
clear mod29_84;
[yd29_84 sort_ind] = sort(yd29_84);
mod29_84 = temp(sort_ind);
clear temp;

bb = size(mod29_85,1);
for i = 1:bb
  temp = yearday([1995 mod29_85(i,2:3)]) + [0 14];
%  temp = yearday([1995 mod29_85(i,2:3)]);
  yd29_85(i) = temp(1,2);
end
aa = find(yd29_85 > 365);
yd29_85(aa) = yd29_85(aa) - 365;

clear temp;
temp = mod29_85;
clear mod29_85;
[yd29_85 sort_ind] = sort(yd29_85);
mod29_85 = temp(sort_ind);
clear temp;

bb = size(mod29_87,1);
for i = 1:bb
  temp = yearday([1995 mod29_87(i,2:3)]) + [0 14];
%  temp = yearday([1995 mod29_87(i,2:3)]);
  yd29_87(i) = temp(1,2);
end
aa = find(yd29_87 > 365);
yd29_87(aa) = yd29_87(aa) - 365;

clear temp;
temp = mod29_87;
clear mod29_87;
[yd29_87 sort_ind] = sort(yd29_87);
mod29_87 = temp(sort_ind);
clear temp;

bb = size(mod38,1);
for i = 1:bb
  temp = yearday([1995 mod38(i,2:3)]) + [0 14];
%  temp = yearday([1995 mod38(i,2:3)]);
  yd38(i) = temp(1,2);
end
aa = find(yd38 > 365);
yd38(aa) = yd38(aa) - 365;

clear temp;
temp = mod38;
clear mod38;
[yd38 sort_ind] = sort(yd38);
mod38 = temp(sort_ind);
clear temp;

% plot all together with different color codings

figure(4)

hold on
lines = 'no '
if lines == 'no '
plot(yd61,mod61,'ko','MarkerSize',5,'MarkerFaceColor',[0 0 0]);
plot(yd94,mod94,'ks','MarkerSize',5,'MarkerFaceColor',[0 0 0]);
plot(yd4,mod4,'k^','MarkerSize',5,'MarkerFaceColor',[0 0 0]);
plot(yd38,mod38,'kd','MarkerSize',5,'MarkerFaceColor',[0 0 0]);
plot(yd29_84,mod29_84,'ko','MarkerSize',5,'MarkerFaceColor',[1 1 1]);
plot(yd29_85,mod29_85,'ks','MarkerSize',5,'MarkerFaceColor',[1 1 1]);
plot(yd29_87,mod29_87,'k^','MarkerSize',5,'MarkerFaceColor',[1 1 1]);
else
plot(yd61,mod61,'ko-','MarkerSize',5,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'Color',[0.6 0.6 0.6],'LineWidth',0.3);
plot(yd94,mod94,'ks-','MarkerSize',5,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'Color',[0.6 0.6 0.6],'LineWidth',0.3);
plot(yd4,mod4,'k^-','MarkerSize',5,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'Color',[0.6 0.6 0.6],'LineWidth',0.3);
plot(yd38,mod38,'kd-','MarkerSize',5,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],'Color',[0.6 0.6 0.6],'LineWidth',0.3);
plot(yd29_84,mod29_84,'ko-','MarkerSize',5,'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0],'Color',[0.6 0.6 0.6],'LineWidth',0.3);
plot(yd29_85,mod29_85,'ks-','MarkerSize',5,'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0],'Color',[0.6 0.6 0.6],'LineWidth',0.3);
plot(yd29_87,mod29_87,'k^-','MarkerSize',5,'MarkerFaceColor',[1 1 1],'Color',[0.6 0.6 0.6],'LineWidth',0.3);
end

ylabel('germination potential(%)')
set(gca,'Xtick',[0:30:360]);

% lump all data together

col1 = [yd61 yd94 yd4 yd29_84 yd29_85 yd29_87 yd38]';
col2 = [mod61 mod94 mod4 mod29_84 mod29_85 mod29_87 mod38]';
mod_data(:,1) = col1;
mod_data(:,2) = col2;

% monthly periods
mon_div = [0 31 59 90 120 151 181 212 243 273 304 334 365];

for i = 1:12
  i
  aa = find((mod_data(:,1) > mon_div(i)) & (mod_data(:,1) <= mon_div(i+1)));
  mod_data(aa,2)
  monthly_median(i) = median(mod_data(aa,2));
  monthly_size(i) = size(aa,1);
  mon_xplot(i) = mon_div(i) + (mon_div(i+1) - mon_div(i))/2; 
end

plot(mon_xplot,monthly_median,'k-','Marker','o','MarkerSize',10,'MarkerFaceColor',[0 0 0],'LineWidth',3)

hold on
xlabel('year-day');
set(gca,'Xtick',[0:30:360]);
axis([0 370 0 100]);
