% ReadSegy : Reads a SEG Y rev 1 fromatted file
%
% Call :
% [Data,TraceHeaderInfo,SegyTraceHeaders,SegyHeader]=ReadSegy(filename);
%
% Implemented using the syntax of the SEGY Y revised format :
% SEGY Y rev 1 DRAFT 6
% http://seg.org/publications/tech-stand/
%
% (C) 2001, Thomas Mejer Hansen, thomas@cultpenguin.com
%
% Extended Textural Header is not yet tested
% If you would like it implemented, please send me an SEGY file with
% that sort of information, as well as a description of the segy file
%
% 
%
function [Data,TraceHeaderInfo,SegyTraceHeaders,SegyHeader]=ReadSegy(filename);
segyid = fopen(filename,'r','b');   % ALL DISK FILES ARE IN BIG
                                    % ENDIAN FORMAT, ACCORDING TO SEG
                                    % Y rev 1
% ans
                    
if ~(exist(filename)==2),
  disp([mfilename,' : ', filename,' does not exist !'])
  Data=[];SegyTraceHeaders=[];SegyHeader=[];
  return
end
% ans
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BINARY HEADERS
SegyHeader=GetSegyHeader(segyid);

% ans
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ TEXTURAL FILE HEADER EXTENSION IF NEEDED

                                   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% READ DATA
%Segy=fread(segyid,4000,'float32');
fseek(segyid,0,'eof'); DataEnd=ftell(segyid);
ans
DataStart=3600+3200*SegyHeader.NumberOfTexturalHeaders;
fseek(segyid,DataStart,'bof');       % Go to the beginning of the file

ntraces=(DataEnd-DataStart)./(240+(SegyHeader.ns)*(SegyHeader.BytesPerSample/8));

disp(['Number of Samples Per Trace=',num2str(SegyHeader.ns)])
disp(['Number of Traces=',num2str(ntraces)])
ans
t=0;




oldclock=clock
while (~(ftell(segyid)>=DataEnd))
  t=t+1;
  if t==1
    
        TraceStart=ftell(segyid);
        [SegyDataTest,SegyTraceHeadersTest]=GetSegyTrace(segyid,TraceStart,SegyHeader.DataFormat,SegyHeader.ns);
        
      
        SegyData(1:ntraces)=SegyDataTest;
        SegyTraceHeaders(1:ntraces)=SegyTraceHeadersTest;
        disp('Inializing space to get our groove into!')
    else
  
        if ((t/200)==round(t/200)), 
            disp(['Reading trace ',num2str(t),' of ',num2str(ntraces)]),
            newclock=clock;
            disp(['taking' num2str(etime(newclock,oldclock)) ' per 200 traces'])    
            oldclock=clock;
        end
        TraceStart=ftell(segyid);
        [SegyData(t),SegyTraceHeaders(t)]=GetSegyTrace(segyid,TraceStart,SegyHeader.DataFormat,SegyHeader.ns);
  
        
        
    end
  
end



% MOVE DATA from SegyData.data to a regular variable
Data=zeros(SegyHeader.ns,t);
for i=1:t
    Data(:,i)=SegyData(i).data;
end

% MOVE common data from structure to regular variables 
S=size(SegyTraceHeaders,2);
cdp=zeros(1,S);
offset=zeros(1,S);
for i=1:S
  if ((i/200)==round(i/200)), disp(['header=',num2str(i),' of ',num2str(S)]);end
    cdp(i)=SegyTraceHeaders(i).cdp;
    offset(i)=SegyTraceHeaders(i).offset;
    Inline3D(i)=SegyTraceHeaders(i).Inline3D;
    Crossline3D(i)=SegyTraceHeaders(i).Crossline3D;
    SourceX(i)=SegyTraceHeaders(i).SourceX;
    SourceY(i)=SegyTraceHeaders(i).SourceY;
    GroupX(i)=SegyTraceHeaders(i).GroupX;
    GroupY(i)=SegyTraceHeaders(i).GroupY;    
    cdpX(i)=SegyTraceHeaders(i).cdpX;
    cdpY(i)=SegyTraceHeaders(i).cdpY;
end

% Create a strcture of common header values if requested
%
if nargout==4,
  TraceHeaderInfo.cdp=cdp;
  TraceHeaderInfo.offset=offset;
  TraceHeaderInfo.Inline3D=Inline3D;
  TraceHeaderInfo.Crossline3D=Crossline3D;
  TraceHeaderInfo.SourceX=SourceX;
  TraceHeaderInfo.SourceY=SourceY;
  TraceHeaderInfo.GroupX=GroupX;
  TraceHeaderInfo.GroupY=GroupY;
  TraceHeaderInfo.cdpX=cdpX;
  TraceHeaderInfo.cdpY=cdpY;
end