function geotimeline(varargin)
%geotimeline(start date,end date,'stages' OR 'all') OR geotimeline(stage name,'stages' OR 'all')
%
%This function produces a schematic timeline showing the ordering and
%relative hierarchy of stages in the geologic timescale.  This function
%can be called either with start and end dates so that it shows the
%divisions in that time interval, or with the name of a particular stage
%in which case it display that stage and the surrounding ones.  If called
%with no arguments, it displays the whole length of geologic time.
%
%If the last argument is the word 'stages', this routine will include the
%lowest level ICS stages/ages.  If the last argument is 'all', the routine
%will include all stages including regional and archaic usages.  If both 
%words are omitted, the timeline displayed will only include divisions down
%to Sub-Epochs and higher.
%
%Note: the columns only indicate relative heirarchy, they do not conform to
%particular categories of time.  This is because different part of the record
%are divided to differing degrees of detail (e.g. sub-era, sub-epochs, etc).

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

%load the geologic database
stages = geologicstages;
%define the heirarchy of standard categories
divisions = {'Eon','Era','Period','Epoch','Sub-Epoch','ICS Stage'};

%faunal is a indicator of whether to display the finest level of stages
%0 = no, 1 = yes.  Defaults to no.
faunal = 0;
%stageprint is a indicator of whether to display the ICS stages
%0 = no, 1 = yes.  Defaults to no.
stageprint = 0;

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

%if there is input, check for the word 'stages','faunal', or 'all'
if nargs > 0
	if (strcmp(lower(va{nargs}),'faunal'))
		%if true, set indicator to yes.
		faunal = 1;
		%set to ignore final argument.
		nargs = nargs - 1;
	elseif (strcmp(lower(va{nargs}),'all'))
		%if true, set indicator to yes.
		faunal = 1;
		%set to ignore final argument.
		nargs = nargs - 1;
   elseif (strcmp(lower(va{nargs}),'stages'))
		%if true, set indicator to yes.
		stageprint = 1;
		%set to ignore final argument.
		nargs = nargs - 1;
   end
end
%this is an indicator of whether to limit the display to an interval of time
%0 = no, 1 = yes.  No is default.
interval = 0;
%if there are 2 arguments, use them to set the start and end time and set
%interval to true
if nargs > 1
	starttime = va{1};
	endtime = va{2};
	%Type Checking
	if ~isnumeric(starttime)
		if (faunal == 1)
			error('First Argument is Not Numeric');
		else
			error('Two or more arguments which are not numeric, and/or word ''stages'' not recognized');
		end
	end
	if ~isnumeric(endtime)
		error('Second argument is neither numeric, nor the word ''stages''.');
	end
	%if out of order, switch them
	if starttime < endtime
	t = starttime;
	starttime = endtime;
	endtime = t;
end
%set interval to true.
interval = 1;
end
%Else If there is one argument, look for matching stage name and set interval around it. 
if nargs == 1
	%if the one argument is numeric, use geowhen to display stages during that date.
	if isnumeric(va{1})
		geowhen(va{1});
		return;
		%else look for the stage matching the given name
	elseif ischar(va{1})
		%the location in the stages array of the record we want.  Set to 0 until found.
		fnd = 0;
		%loop over all stages looking for a match
		for a = 1:length(stages)
			if strcmp(stages(a).Name,va{1})
				fnd = a;
				break;
			else
				%loop over the aliases
				for b = 1:length(stages(a).AKA)
					if strcmp(stages(a).AKA{b},va{1})
						fnd = a;
						break;
					end
				end
				%If stated name was an alias, display message indicating primary name.
				if fnd~=0
					disp(' ');
					disp(['  ---  ' va{1} ' Also Known As ' stages(a).Name '  ---']);
					break;
				end
			end
		end
		%if not found, give an error.
		if fnd==0
			error(['Stage Not Recognized: ' va{1}]);
		end
		%Type of category for this stage.
		CatName = stages(fnd).Category;
		%Find this category name in the list of divisions and record it's number
		%if not in list, indicate as length of list + 1.
		CatNum = length(divisions) + 1;
		for a = 1:length(divisions)
			if strcmp(lower(CatName),lower(divisions{a}))
				CatNum = a;
			end
		end
	
		%if division is Epoch or lower, force it to use the stage list.
		if CatNum >= length(divisions)-2
			stageprint = 1;
		end
	
		%indicate the index numbers in the data structure of the stages
		%before and after the current stage.  0 until found.
		before = 0;
		follow = 0;
		%if category was in the divisions list, we want the stages of same category
		%that come before and after.
		if CatNum <= length(divisions)
			%loop over stages looking for same category
			for b = 1:length(stages)
				if strcmp(stages(b).Category,CatName)
					if stages(b).Start == stages(fnd).End
						follow = b;
					elseif stages(b).End == stages(fnd).Start
						before = b;
					end
				end
			end
			%Duration of stage, used if before or after not found.
			td = stages(fnd).Start - stages(fnd).End;
			%set start of interval
			if before ~= 0
				starttime = stages(before).Start;
			else
				starttime = stages(fnd).Start + td;
			end
			%set end of interval
			if follow ~= 0
				endtime = stages(follow).End + 0.0001;
			else
				endtime = stages(fnd).End - td;
			end
			%set indicator of interval to true.
			interval = 1;
			%if category was not on divisions list then simply take it's length
			%and add to time before and after   
		else
			td = stages(fnd).Start - stages(fnd).End;
			starttime = stages(fnd).Start + td;
			endtime = stages(fnd).End - td;
			%set indicator of interval to true.
			interval = 1;
		end
	end				
end

%Figures out the category number for all of the stages in the database.
%Used for sorting.
cats = lower(strvcat(stages(:).Category));
[stages.CatNum] = deal(length(divisions)+1);
for a = 1:length(divisions);
	stn = strmatch(lower(divisions{a}),cats);
	[stages(stn).CatNum] = deal(a);
end

%if we don't want the lowest stages, then remove them from the record.
if ~faunal
	ct = cat(1,stages(:).CatNum);
	f = find(ct <= length(divisions));
	stages = stages(f);
end
if ~stageprint & ~faunal
	ct = cat(1,stages(:).CatNum);
	f = find(ct < length(divisions));
	stages = stages(f);
end

%Create lists of the start and end times of each stage
tslist = cat(1,stages(:).Start);
telist = cat(1,stages(:).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 = stages(re(end:-1:1));

%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 = {};
%loop over all stages
for a = 1:length(stages)
	%if we haven't already added this time...
	if stages(a).Start ~= tlast
		%if set to an interval, check bounds
		if interval
			if ((stages(a).Start > starttime)&(stages(a).End < endtime))|(stages(a).End < starttime)
				times{end+1,1} = stages(a).Start;
				times{end,2} = stages(a).dStart;
				tlast = stages(a).Start;
			end
			%if beyond the endtime, set endtime to the start of first following stage and
			%exit the loop.
			if stages(a).Start < endtime
			endtime = stages(a).Start;
			break;
		end
		%else, add every time.
	else
		times{end+1,1} = stages(a).Start;
		times{end,2} = stages(a).dStart;
		tlast = stages(a).Start;
	end			
end
end

%if set to only a specific interval, remove all stages that
%don't overlap on the interval
if interval
	a = 1;
	while a <= length(stages)
		if ~((stages(a).End <= starttime) & (stages(a).Start >= endtime))
			stages = stages([1:(a-1),(a+1):end]);
			a = a - 1;
		end
		a = a +1;
	end	
end	

%FINALLY... we start the output.
disp(' ');
%loop over all the transition times.
for a = 1:length(times)
	%b is an index to the position in the stages array while seraching
	b = 1;
	%current is an array of stages that are present at the current time.
	%c is an index in the current array.
	current = stages(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(b).Start >= times{a,1})
		%if stage ends after the current time then add it to the list.
		if (stages(b).End < times{a,1})
			current(c) = stages(b);
			c = c + 1;
			%else remove it from the stages list as no longer neccesary.
		else
			stages = stages([1:(b-1),(b+1):end]);
			b = b - 1;
		end
		b = b + 1;
		%if out of stages to sort, then break.
		if (b > length(stages))
			break;
		end
	end
	%loop over the current array and sort it by position
	%in the heirarchy of stages.
	for j = 1:length(current)
		t1 = current(j);
		for b = 1:length(current)
			t2 = current(b);
			if t1.CatNum < t2.CatNum
				t = current(j);
				current(j) = current(b);
				current(b) = t;
				t1 = t2;
			end
		end
	end
	%str is the string to be printed.
	str = '';
	%loop over the current array and add the names to str.
	for c = 1:length(current)
		%if name was previously printed and it is a defined heirarchy stage
		%then we print only blanks in it's place.
		if current(c).CatNum <= length(divisions) & current(c).Start~=times{a,1}
			str = [str, blanks(length(current(c).Name)),' - '];
			%else, just add the name   
		else
			str = [str, current(c).Name,' - '];
		end
	end
	%dst is the error associated with the current time
	dst = times{a,2};
	%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

	%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,1});
	%if there is an error estimate, add a seperator space.
	if ~isempty(dst)
		dst = [dst,' '];
	end
	%format the final string.
	stra = [st ' ' dst];
	str = [blanks(20 - length(stra)) stra 'Ma: ' str];
	%display it.
	disp(str);
end
disp(' ');