function [st,ed,rec] = geowhen(w, geostages)
%[start, end, record] = geowhen('stage name' OR time [Ma], geologicstages)
%May be called with any of 0-3 outputs.
%geologicstages argument is optional.
%
%This function, when called with a numeric input, prints to the screen the periods,
%epochs, stages, etc. that were occuring on that date in millions of years
%ago (Ma).
%
%When called with a string input it responds differently depending on whether or
%not outputs are requested.  If no outputs are requested, information describing
%the stage is displayed to the screen.  If one are more outputs are specified
%it returns the start date, end date and database record to fill those outputs
%upto the number requested.  The database record is a structure of the type
%defined in geologicstages.m.  When outputs are requested, no information is 
%displayed to the screen.
%
%If called with a "geologicstages" input parameter, then this must be a structure
%of the type created by geologicstages.m.  This record will be used instead of
%calling geologicstages.m.  This parameter is provided for backward compatibility
%and can improve access time if geowhen is to be called repeatedly.
%
%If this program can find no stage name that exactly matches the one given, it will
%look to see if any stages start with the specified string.  (For instance returning
%Cretaceous when given 'Cret').  If multiple names match it gives a list of such
%stages.  If no names start with the specified string then the program lists
%those stages which are closest in alphabetical order.  In all cases in which the
%string fails to match exactly, the output is directed to the screen and specified
%output arguments are ignored.  The matching functions are performed by calling
%geolist.m

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

%load the database
if nargin < 2
	stages = geologicstages;
else
   stages = geostages;
end

%defines categories of time higher than stages for use with ordering heirarchies.
divisions = {'Time','Eon','Era','Period','Epoch','Sub-Epoch'};

%Execute this block if input is a string
if ischar(w)
	%Indicates the record of interest, set to 0 until it is found.
	fnd = 0;
	%search for the record by name ...
	for a = 1:length(stages)
		%... in the primary name
		if (strcmp(lower(stages(a).Name),lower(w)))
			fnd = a;
			break;
		end
		%... in the aliases
		if ~isempty(stages(a).AKA)
			for tc = 1:length(stages(a).AKA)
				if (strcmp(lower(stages(a).AKA{tc}),lower(w)))
					fnd = a;
					a = length(stages);
					break;
				end
			end
			if (fnd~=0)
				break;
			end
		end
	end

	%If the stage was not found, run the listing routine to show close matches.
	if (fnd==0)
		disp(' ');
		disp('-- Stage Name Not Recognized.  Did you mean? --');
		disp(' ');
		geolist(w);
		disp(' ');
		return
	end

	%if output is requested, set the outputs and return.
	if (nargout > 0)
		st = stages(fnd).Start;
		ed = stages(fnd).End;
		rec = stages(fnd);
		return;
		%else, display to the screen all the information.
	else
		disp(' ');
		%set the width of screen display size.
		ln = 80;
		%store the record we found in "stage"
		stage = stages(fnd);
		
		%The title line
		str = [stage.Name ' ' stage.Category];
		disp([blanks(floor(ln/2 - length(str)/2)), str]);
		
		%the aliases
		str = [];
		for a = 1:length(stage.AKA)
			str = [str, stage.AKA{a}, ', '];
		end
		if ~isempty(str)
			str = str(1:(end-2));
			str = ['(' str ')'];
			disp([blanks(floor(ln/2 - length(str)/2)), str]);
		end
	
		%date information
		starts = stage.Start;
		ends = stage.End;
		dst = stage.dStart;
		dend = stage.dEnd;
		
		%process the different kinds of error formats into strings
		if isempty(dend)
			dend = '';
		end
		if isempty(dst)
			dst = '';
		end
		if isnumeric(starts)
			starts = num2str(starts);
		end
		if isnumeric(ends)
			ends = num2str(ends);
		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 iscell(dend)
			d1 = dend{1};
			d2 = dend{2};
			if isnumeric(d1)
				d1 = num2str(abs(d1));
			end
			if isnumeric(d2)
				d2 = num2str(abs(d2));
			end
			dend = ['+' d1 '/-' d2];
		end
		if isnumeric(dend)
			dend = ['+/- ' num2str(dend)];
		end
	
		%write out dates and errors	  
		str = ['From ' starts];
		if ~isempty(dst)
			str = [str ' ' dst];
		end
		str = [str ' To ' ends];
		if ~isempty(dend)
			str = [str ' ' dend];
		end
		disp(' ');
		str = [str, ' Ma'];
		disp([blanks(floor(ln/2 - length(str)/2)), str]);
		
		%Describe Definitions
		if (~isempty(stage.Start_Defined_By)) | (~isempty(stage.End_Defined_By))
			disp(' ');
		end
		if ~isempty(stage.Start_Defined_By)	  
			wrapline(['Start Defined By:  ' stage.Start_Defined_By],ln);
		end
		if ~isempty(stage.End_Defined_By)	  
			wrapline(['End Defined By:  ' stage.End_Defined_By],ln);
		end
	
		%Describe Sources
		if (~isempty(stage.Start_Source)) | (~isempty(stage.End_Source))
			disp(' ');
		end
		if ~isempty(stage.Start_Source)	  
			wrapline(['Start Based On:  ' stage.Start_Source],ln);
		end           
		if ~isempty(stage.End_Source)	  
			wrapline(['End Based On:  ' stage.End_Source],ln);
		end
	
		%Describe GSSP	  
		if ~isempty(stage.GSSP)	  
			disp(' ');
			wrapline(['Global Stratotype Section and Point:  ' stage.GSSP],ln);	  
		end
	
		%Add Comment
		if ~isempty(stage.Comment)  
			disp(' ');
			wrapline(['Comment:  ' stage.Comment],ln);
		end
      
      %Determine Current Stage's Division
      sdiv = length(divisions) + 1;
		for b = 1:(sdiv-1)
			if strcmp(stage.Category,divisions{b})
				sdiv = b;
			end
		end

		%This section searches for related stages.
		prior = [];
		after = [];
		overlap = [];
		superset = [];
		subset = [];
		%temp variables for the start and end of this stage
		sst = stages(fnd).Start;
		sed = stages(fnd).End;
		%loop over all stages other than the current one, and find 
		%what categories they fall in.
		for b = [1:(fnd-1),(fnd+1):length(stages)]
			if (stages(b).Start >= sst) & (stages(b).End <= sed)
				subset = [subset,b];
			end
			if (stages(b).Start <= sst) & (stages(b).End >= sed)
				superset = [superset,b];
			end
			if ((stages(b).Start > sst) & (stages(b).End < sst)) | ((stages(b).Start > sed) & (stages(b).End < sed))
				if (~ismember(b,subset))&(~ismember(b,superset));
					overlap = [overlap,b];
				end
			end
			if (stages(b).Start == sed)
				after = [after,b];
			end
			if (stages(b).End == sst)
				prior = [prior,b];
			end		 
		end
		%use the standard divisions to select only those stages of the next
		%finest ordering when listing subdivisions of the current stage.
      %First remove any equal length stages of higher order...
      j = 1;
      while j <= length(superset);
         if stages(superset(j)).Start - stages(superset(j)).End == stage.Start - stage.End
            for b = 1:(sdiv-1)
					if strcmp(stages(superset(j)).Category,divisions{b})
                  superset = superset([1:j-1,j+1:end]);
                  break;
					end
            end
         end
         j = j + 1;
		end
 		maxd = length(divisions) + 1;
		%... then find the subdivision of highest order ...        
      for a = 1:length(superset);
			for b = 1:(maxd-1)
				if strcmp(stages(superset(a)).Category,divisions{b})
               maxd = b;
				end
			end
		end
		%... then include only divisions of the same order 
		if maxd <= length(divisions)
			ls = [];
			for a = 1:length(superset);
				if strcmp(stages(superset(a)).Category,divisions{maxd})
					ls = [ls,superset(a)];
				end
			end
			superset = ls;
		end
      
      j = 1;
      while j <= length(subset);
         if stages(subset(j)).Start - stages(subset(j)).End == stage.Start - stage.End
            good = -1;
            for b = 1:min([sdiv,length(divisions)])
               if strcmp(stages(subset(j)).Category,divisions{b})
                  good = 1;
               end
            end
            if good == -1
               subset = subset([1:j-1,j+1:end]);
            end
         end
         j = j + 1;
		end      
      
		%if any overlaps found, then display them.
		if ~isempty(overlap) 
			disp(' '); 
			disp(' -- Overlaps With --');
			displaylist(stages, overlap);
		end
	
		%if any stages are superior to this one, then display them.
		if ~isempty(subset) 
			disp(' '); 
			disp(' -- A Substage of --');
			displaylist(stages, subset);
		end
	
		%if any stages preceed this one, display them.
		if ~isempty(prior) 
			disp(' ');
			disp(' -- Preceeded By --');
			displaylist(stages, prior); 
		end
	
		%if any stages follow this one, display them. 
		if ~isempty(after) 
			disp(' ');
			disp(' -- Followed By --');
			displaylist(stages, after);
		end
	
		%if this stage has any subdivisions, display them.	 
		if ~isempty(superset) 
			disp(' ');
			if maxd <= length(divisions)
				disp([' -- Contains the Following ' divisions{maxd} 's --']);
			else		 
				disp(' -- Contains the Following Stages --');
			end
			displaylist(stages, superset);
		end
		disp(' ');
	end

	%ELSE, IF the function is called with a numeric parameter,
	%then search for stages occuring at that time and display them.
elseif (isnumeric(w))
	cs = [];
	%Find stages that start at or before the given date,
	%and who end after it.
	for a = 1:length(stages)
		if ((stages(a).Start >= w)&(stages(a).End < w))
			cs = [cs,a];
		end      
	end

	%figure out the durations of all such stages.
	for b = 1:length(cs)
		df(b) = stages(cs(b)).Start - stages(cs(b)).End;
	end

	%sort them with longest stage first
	[s,re] = sort(df);
	df = [cs(re(end:-1:1));df(re(end:-1:1))]';
	
	%figure out how much room we need for the output
	disp(' ');
	cmax = 0;
	nmax = 0;
	for a = 1:length(df(:,1))
		if (length(stages(df(a,1)).Category)>cmax)
			cmax = length(stages(df(a,1)).Category);
		end
		if (length(stages(df(a,1)).Name)>nmax)
			nmax = length(stages(df(a,1)).Name);
		end
	end 

	%for each stage occurring at this time display the output.
	for a = 1:length(df(:,1))
		%temp variables
		er = stages(df(a,1)).Category;
		nm = stages(df(a,1)).Name;
		starts = stages(df(a,1)).Start;
		ends = stages(df(a,1)).End;
		if (isnumeric(starts))
			starts = num2str(starts);
		end
		if (isnumeric(ends))
			ends = num2str(ends);
		end
		%construct the string and display it
		str = blanks(cmax + 2 - length(er));
		str = [str,er, ': ' nm];
		str = [str,blanks(nmax + 2 - length(nm)),'(',ends,' - ',starts, ')'];
		disp(str);
	end
	disp(' ');
end


%sub function used to write lists of stages.
function displaylist(stages,lst)
%stages is the GeologicStages data structure, and lst is a list
%of indices in that data structure that one wants output.

divisions = {'Time','Eon','Era','Period','Epoch','Sub-Epoch','ICS Stage'};
%sort list so that they are ordered first by division of geologic time and 
%then by oldest first.  Parts of the named heirachy always precede regional
%designations.

%for each stage dsignate which division it belongs to.
for a = 1:length(lst)
   divis(lst(a)) = length(divisions) + 1;
   for j = 1:length(divisions)
		if strcmp(stages(lst(a)).Category,divisions{j})
         divis(lst(a)) = j;
		end
   end
   %if it doesn't match to any division on the primary
   %heirachy, number it based on it's name.
   if (divis(lst(a)) == length(divisions) + 1)
      cat = lower(stages(lst(a)).Category);
      for k = 1:length(cat)
         divis(lst(a)) = divis(lst(a)) + (cat(k)-'a')/26^k;
      end
   end
end
%sort the stages
for a = 1:length(lst)
	t1 = stages(lst(a));
	for b = 1:length(lst)
		t2 = stages(lst(b));
		%first by division rank
      if(divis(lst(a))<divis(lst(b)))
			t = lst(a);
			lst(a) = lst(b);
         lst(b) = t;
			t1 = t2;
      elseif (divis(lst(a))>divis(lst(b))) 
         t1=t2;
      end
      %then by start order
      if (t2.Start < t1.Start)
			t = lst(a);
			lst(a) = lst(b);
			lst(b) = t;
			t1 = t2;
		end
		%then by length
		if (t2.Start == t1.Start) & (t2.End > t1.End)
			t = lst(a);
			lst(a) = lst(b);
			lst(b) = t;
			t1 = t2;
      end
	end
end

%figure out the width of display needed to write the lines
cmax = 0;
nmax = 0;
for a = 1:length(lst)
	if (length(stages(lst(a)).Category) > cmax)
		cmax = length(stages(lst(a)).Category);
	end
	if (length(stages(lst(a)).Name) > nmax)
		nmax = length(stages(lst(a)).Name);
	end
end	

%display the list
for a = 1:length(lst)
	tmp = stages(lst(a));
	disp(['  ' tmp.Category ':' blanks(cmax-length(tmp.Category)) ' ' tmp.Name blanks(nmax-length(tmp.Name)) ' (' num2str(tmp.Start) ' - ' num2str(tmp.End) ')']);
end


%subfunction used to display lines that may be too
%long for the specified width of screen.
function wrapline(text,ln)
%text is the string to display and ln is the screen width limit.

%fst is a indicator designating the first iteration of the loop.
fst = 1;

%while there is still text in the string to display
while (~isempty(text));
	%if first loop, indent less.
	if fst
		str = '  ';
	else
		str = '     ';
	end
	%not "first" loop from here on.
	fst = 0;
	%If text is shorter than a full line... 
	if length(text) < (ln - length(str));
		%set ln to the number of characters left to print
		ln = length(text) + length(str);
		%st is number of character to write on this pass
		%in this case, all of them.
		st = ln-length(str);
	else
		%else, find the best breaking point.
		st = ln-length(str);
		%characters considered good breaking points.
		while (~ismember(text(st),' /,:?'))&(st>1)
			st = st - 1;
		end
	end
	%add to str the text we want on this pass.
	str = [str,text(1:st)];
	%drop the bit that is about to be written.
	text = text((st+1):end);
	%write it.
	disp(str);
end