function h_v=loadsegyheaders(filename,field,value1,value2)
% LOADSEGYHEADERS(filename,field,value1,value2)
%
% reads header values from a SEG-Y file into a structure
% FUNCTION ARGUMENTS: 
% filename (string): Mandatory. The path and filename of the SEG-Y file.
% field (string), value1 (number), value2 (number): Optionals.
% field must be the exact name of a header field from headers.m
% Use value1 and value2 if you want to load only the headers for which the
% value of a certain field is between (or equal to) value1 and value2. If the
% value2 argument is not present, only headers for which the values of 'field'
% is EQUAL to value1 will be read. If only the filename argument is present,
% every header will be read until the end of the file. Examples:
% H=loadsegyheaders('/nick/tape1/file23') reads all headers
% H=loadsegyheaders('/nick/tape1/file23','channel',1) reads only the vertical
%			 component headers (field 'channel' has value 1)
% H=loadsegyheaders('/nick/tape1/file23','offset',100,1000) reads only the
% traces whose offset is between 100 and 1000.
% FUNCTION OUTPUT: a structure of 'horisontal' vectors holding the header
% values. The last field of the structure ('S.locator') holds the number that
% tells us to which trace, numbered from the begginning of the file, the header
% belongs. For more help type 'type loadsegyheaders' at the Matlab prompt.
% See also LOADSEGYTRACES, HEADERS

% Beware, if you deal with reel files, the variable reel should be 3600.
% The headers are written as a cell array in headers.m
% that must be somewhere in the Matlab path.
% The cell describing the headers need not contain all the possible 
% information, nor does it need to be in a certain order. You can thus build 
% 'reduced headers' that contain only the fields that you need.
% The program assumes that, according to the SEG-Y standard, the number of
% samples in each trace is a 2-byte integer written on the bytes 115 and 116 of
% the trace header, and  that the trace values are 4 bytes long each.
% For the whole thing to run faster, compile with 'mcc -x loadsegyheaders.m'
% Last revised: Nick Vlad, Jan. 2001

reel=3600; % length of reel header (skipped)
ns_pos=115; % position in the trace header of the number of samples
ns_format='int16'; % position in the trace header of the number of samples
samples_format='float32'; % assumption about the format
sample_space=4; % assumption about the number of bytes a sample occupies
		% on disk (depends on sample_format)
h_l=240; % trace header length assumption
t_h_s=headers; % loads the header structure from headers.m
h_v=cell2struct(cell(size(t_h_s(:,1))),t_h_s(:,1),1); % future header values
h_v.locator=[]; % will hold positions of traces.
                % Assumption no header field is called 'locator'
if nargin==1
	f_o_c=cell2mat(t_h_s(1,1)); % field of choice
	f_o_c_pos=cell2mat(t_h_s(1,2)); % f_o_c position in the header
	f_o_c_frmt=cell2mat(t_h_s(1,3)); % f_o_c format in the header
else
	f_o_c=field;
	for i=1:length(t_h_s(:,1))
		if size(cell2mat(t_h_s(i,1)))==size(f_o_c)
			if cell2mat(t_h_s(i,1))==f_o_c
				f_o_c_pos=cell2mat(t_h_s(i,2));
				f_o_c_frmt=cell2mat(t_h_s(i,3));
				break
			end
		end
	end
end

FID = fopen(filename,'r','ieee-be'); % The SEG-Y file is now OPEN for reading
if FID==-1 error(cat(2,'Cannot open ',filename));end
status = fseek(FID,0,'eof');filesize=ftell(FID); % determine file size
status = fseek(FID,reel+ns_pos-1,'bof'); % go to number of samples position
ns=fread(FID,1,ns_format); % read the  number of samples in the trace
count=0;

% we will select the headers that we want to read
for i=1:floor(filesize/(h_l+ns*sample_space))
	if reel+i*(h_l+ns*sample_space)>filesize
		break
	end
	status=fseek(FID,reel+(i-1)*(h_l+ns*sample_space)+f_o_c_pos-1,'bof');
	temp=fread(FID,1,f_o_c_frmt);
	if nargin==1
		h_v.locator(i)=i; % mark all headers for loading
	elseif (nargin==3)&(temp==value1)
		count=count+1; h_v.locator(count)=i;
	elseif (nargin==4)&(value1<=temp)&(temp<=value2)
		count=count+1; h_v.locator(count)=i;
	end
end

% Now reading the headers: 
for i=1:length(h_v.locator)
	for j=1:length(t_h_s(:,1)) % 'locator' must be the last field in h_v
		status=fseek(FID,reel+(h_v.locator(i)-1)*(h_l+ns*sample_space)+...
					cell2mat(t_h_s(j,2))-1,'bof');
		eval(strcat('h_v.',cell2mat(t_h_s(j,1)),'(',num2str(i),...
			')=',num2str(fread(FID,1,cell2mat(t_h_s(j,3)))),';'));
	end
end
fclose(FID);




