function georegion(varargin)
%georegion(name1, name2, ...)
%
%This function produces a timeline comparing regional stratigraphic systems
%to the ICS 2004 geologic timescale.  The input to the function are names of 
%regional stratigraphic systems, where the recognized names are defined in
%stratigraphicregions.m.  The stratigraphicregions.m and geologicstages.m files
%must be available for this function to work.

%Created by Robert A. Rohde (geowhen@robertrohde.com)
%Last Modified January 18, 2005

%load regional information
regions = stratigraphicregions;

%store the inputs in editable variables.
nargs = nargin;
va = varargin;

%processes the inputs as either system names or system id numbers.
if ischar(va{1})
	for a=1:length(va)
	   system(a) = 0;
      for b = 1:length(regions)
         %compare input argument to system name
	   	if strcmp(lower(regions(b).name),lower(va{a}))
            system(a) = b;
         end
         %also compare it to aliases.
	   	for c = 1:length(regions(b).alias)
		   	if strcmp(lower(regions(b).alias{c}),lower(va{a}))
   		   	system(a) = b;
   			end
   		end               
      end      
      %if name is not recognized, throw an error.
   	if system(a)==0
      	error(['Regional System "' va{a} '" not found.']);
   	end
   end
elseif isnumeric(va{1}) 
   system= [];
   %if numeric, simply concatenate list
   %no error checking is performed.
   for a =1:length(va)
      system = [system,va{a}];
   end
end

%The first and second systems are reserved for ICS general and specific
%timescales.  First we load these.
stages{1} = loadsystem(regions(1));
stages{2} = loadsystem(regions(2));
regiontitle{1} = regions(1).name;
regiontitle{2} = regions(2).name;
%now load regional systems.
for a = 1:length(system)
	stages{a+2} = loadsystem(regions(system(a)));
	regiontitle{a+2} = regions(system(a)).name;
   %figure out the start and end time for this system.
   starttime(a) = max(regions(system(a)).dates(:,1));   
   endtime(a) = min(regions(system(a)).dates(:,2));   
end   

%figure out the start and end times for the whole comparison.
starttime = max(starttime);
endtime = min(endtime);

%Create lists of the start and end times of each stage
for a = 1:length(stages)
   tslist = cat(1,stages{a}(:).Start);
   telist = cat(1,stages{a}(:).End);
   
   %create a superlist that emphasizes the start date but
	%also gives weight to the end date so that longer stages
	%appear before short stages.
	tlist = tslist + telist*1e-8;
	%use this list to sort the stage list.
	[tlist2,re] = sort(tlist);
	stages{a} = stages{a}(re(end:-1:1));
end

%create a list of the start times, for stages that we want to include
%tlast is the last time added to the list, to avoid repetitions.
%initial value set to impossibly long ago.
tlast = 1e6;
%cell array of start times and their errors
times = [];
terr = {};
%loop over all stage systems
for j = 1:length(stages)
   %loop over all stages in a system
	for a = 1:length(stages{j})
		%if we haven't already added this time...
		if stages{j}(a).Start ~= tlast
			if ((stages{j}(a).Start > starttime)&(stages{j}(a).End < endtime))|(stages{j}(a).End < starttime)
				times(end+1) = stages{j}(a).Start;
				terr(end+1) = stages{j}(a).dStart;
				tlast = stages{j}(a).Start;
			end
			%if beyond the endtime, set endtime to the start of first following stage and
			%exit the loop.
			if stages{j}(a).Start < endtime
            if (j>2) 
               endtime = stages{j}(a).Start;
            else
               times(end) = [];
            	terr(end) = [];   
            end
				break;
			end
   	end
	end
end

%sort out the times and errors
[times,I] = sort(times);
terr = terr(I);
%remove duplicate times, while preserving meaningful errors.
c = 2;
while c<=length(times)
   if (times(c) == times(c-1))
      times(c) = [];
      if ~isempty(terr(c)) 
         terr(c) = [];
      else
         terr(c-1) =[];
      end
   else
      c = c + 1;
   end
end
%place in oldest first order.
times = times(end:-1:1);
terr = terr(end:-1:1);

%Remove all stages that don't overlap on the interval
for j = 1:length(stages);
   a = 1;
   %loop over stages within each system
   while a <= length(stages{j})
      if ~((stages{j}(a).End <= starttime) & (stages{j}(a).Start >= endtime))
         stages{j} = stages{j}([1:(a-1),(a+1):end]);
         a = a - 1;
      end
      a = a +1;
   end	
end

%Determine the max size needed for each column.
for j = 1:length(stages);
   col = 0;
   for a = 1:length(stages{j})
      if length(stages{j}(a).Name) > col & ~strcmp(stages{j}(a).Name(1:3),'End')
         col = length(stages{j}(a).Name);
      end
   end
	colwidth(j) = col;   
end
      
%FINALLY... we start the output.
disp(' ');
str = blanks(24);
for j = 1:length(stages)
   rt = regiontitle{j};
   %update column width if it can't accomodate region title.
   if colwidth(j) < length(rt)
      colwidth(j) = length(rt);
   end
   %write the region title.
   str = [str,blanks(floor((colwidth(j)-length(rt))/2)),rt,blanks(ceil((colwidth(j)-length(rt))/2)),blanks(3)];
end
disp(str);

%indicates when a column has started printing output, all set to false (0)
%bt default.
good(1:length(stages)) = 0;

%To avoid spurious output, remove times corresponding to the present day.
f = find(times>0);
times=times(f);

%loop over all the transition times.
for a = 1:length(times)
   str = '';
   %find the stage to print in each system.
   for j = 1:length(stages)
      %b is an index to the position in the stages array while searching
      b = 1;
      %if we have already exhausted all of the stages in the current system,
      %print a placeholder.
      if isempty(stages{j}) 
         if (j>2)
            sta = [blanks(floor((colwidth(j)-7)/2)) 'use ICS'];
         else
            sta = [blanks(floor((colwidth(j)-8)/2)) 'no stage'];
         end   
         str = [str,sta,blanks(colwidth(j)-length(sta)),' - '];
      else
         %current is an array of stages that are present at the current time.
         %c is an index in the current array.
         current = stages{j}(1);
         c = 1;
         %loop over stages to find those active during the current time.
         %loop exits when stages are starting after the current time
         %since stages list was time sorted on start times.
         while (stages{j}(b).Start >= times(a))
            %if stage ends after the current time then add it to the list.
            if (stages{j}(b).End < times(a))
               current(c) = stages{j}(b);
               c = c + 1;
            elseif (stages{j}(b).End == stages{j}(b).Start) & (stages{j}(b).End==times(a))
               current(c) = stages{j}(b);
               c = c + 1;
               %else remove it from the stages list as no longer neccesary.         
            else         
               stages{j} = stages{j}([1:(b-1),(b+1):end]);
               if (c==1) & ~isempty(stages{j})
                  current(c) = stages{j}(1);
               else
                  %if isempty, print a placeholder.
                  if isempty(stages{j}) 
                     if (j>2)
                        sta = [blanks(floor((colwidth(j)-7)/2)) 'use ICS'];
                     else
                        sta = [blanks(floor((colwidth(j)-8)/2)) 'no stage'];
                     end   
                     str = [str,sta,blanks(colwidth(j)-length(sta)),' - '];
                     current(c) = [];
                     break;
                  end
               end
               b = b - 1;
            end
            b = b + 1;
            %if out of stages to sort, then break.
            if (b > length(stages{j}))
               break;
            end
         end
         %str is the string to be printed.
         %loop over the current array and add the names to str.
         for c = 1:length(current)
            %if name was previously printed
            %then we print only blanks in it's place.
            if current(c).Start>=times(a) & current(c).End~=times(a)
               if current(c).Start~=times(a)
                  str = [str, blanks(colwidth(j)),' - '];
                  %else, just add the name   
               else
                  str = [str, current(c).Name,blanks(colwidth(j)-length(current(c).Name)),' - '];
               end
               good(j) = 1;
            else
               %else, if no stage in current system
               if (good(j) == 1)
                  %if system has already started, print place holders.
                  if (j>2)
                     sta = [blanks(floor((colwidth(j)-7)/2)) 'use ICS'];
                  else
                     sta = [blanks(floor((colwidth(j)-8)/2)) 'no stage'];
                  end   
                  %else, if not started, print blanks.
               else
                  sta = [blanks(colwidth(j))];
               end
               str = [str,sta,blanks(colwidth(j)-length(sta)),' - '];
            end   
         end
      end
   end
   
   dst = terr{a};
	%process the possible formats for the error.
	if isempty(dst)
		dst = '';
	end
	if iscell(dst)
		d1 = dst{1};
		d2 = dst{2};
		if isnumeric(d1)
			d1 = num2str(abs(d1));
		end
		if isnumeric(d2)
			d2 = num2str(abs(d2));
		end
		dst = ['+' d1 '/-' d2];
	end
	if isnumeric(dst)
		dst = ['+/- ' num2str(dst)];
	end
   
   if ~isempty(dst)
		dst = [dst,' '];
	end
   
   %drop the last three characters ' - ' which aren't needed after the last name.
	str = str(1:(end-3));
	%st equals the string version of the time.
	st = num2str(times(a));
	%format the final string.
	stra = [st ' ' dst];
	str = [blanks(20 - length(stra)) stra 'Ma: ' str];
	%display it.
	disp(str);
end
disp(' ');

%helper function that load information about the regional system
function stagelist = loadsystem(region)

%fix start and end times.
starttime = max(region.dates(:,1));
endtime = min(region.dates(:,2));

%number of stages in system.
sti = 1;
eti = length(region.member);

%place information in stages array
for a = sti:eti
   stages(a).Name = region.member{a-sti+1};
   stages(a).Start = region.dates(a-sti+1,1);
   stages(a).dStart = region.error(a-sti+1,1);
	stages(a).End = region.dates(a-sti+1,2);
end
%add placeholders for end dates that don't mark the beginning of new
%stage in the same system.
slist = [stages.Start];
elist = [stages.End];
%first find the unmatched ends.
[mlist,I] = setdiff(elist,slist);
st2 = eti+1;
et2 = eti + length(mlist);
%now add placeholders
for a = st2:et2;
   stages(a).Name = ['End' stages(I(a-st2+1)).Name];
   stages(a).Start = mlist(a-st2+1);
   stages(a).dStart = region.error(I(a-st2+1),2);
	stages(a).End = mlist(a-st2+1);
end   

%return array.
stagelist = stages;