
function [seis,hdr,reel]=read_segy(filename,shotno,tlim,vred,rlim,rpass,decfac,chan)
%function [seis,hdr,reel]=read_segy(filename,shotno,tlim,vred,rlim,rpass,decfac,chan);
%                                              *      *    *    *
%Reads standard SEGY in file FILENAME
%SEIS is the seismic records and HDR and REEL are structures containing 
%trace and reel header data.
%Note: A non-standard millisecond field is read from bytes 207-208 of the trace
%      header and is read into structure field hdr.msec and is NOW added to
%      hdr.bigtime (Changed 1/24/05)
%
%Optional input arguments can be used to return subsets of data
%SHOTNO  - Shot (field record) numbers to read if positive
%          Trace number within the file if negative 
%TLIM    - Two element vector containing minimum and maximum times (s) relative
%          to the shot time
%VRED    - A reduction velocity (km/s) to apply to traces before applyin TLIM
%RLIM    - Two element vector containing minimum and maximum ranges to return
%          SEGY convention is +ve range when range increases with shot number
%RPASS   - If file contains more than on section of data within the range limits
%          return on the the RPASS'th if RPASS>0
%Values of [], 0/[0 0], or no input mean the arguments are inoperative
%DECFAC  - Integer Decimation factor to reduce trace length using Matlab
%          function resample.m.  A value of 2 would apply an anti-alias filter
%          and return every other sample
%
% There are also 3 global scalars which affect how READ_SEGY works.  Each of 
% these has a default value (in square brackets) when not set externally
% MONOTONIC - 1 if shot numbers increase monotonically throughout the file [0]
% RANGESIGNED - 1 if ranges in SEGY headers are properly signed [1]
%               WHOI OBH ranges are signed but OBS and ORB data are not
%               If ~RANGESIGNED this program will assign ranges based on the
%               sign of the differene in range between adjacent shots.
%               Changes in sign at large ranges are eliminated for the first 
%               and last 5% of shots (sometimes the ship turns near line ends).
%               If the range sign still changes more than once (i.e., complex
%               geometry), all ranges are made positive
% DELTARANGEMAX - Maximum spacing in metres of shots (0 - not known) [0]
% Note: global variables are capitalized
%
% SEGY files should be big-endian but this routine will attempt to read 
% little-endian segy aswell
% 
%@EEEH  modified to read only the selected channel

global MONOTONIC RANGESIGNED DELTARANGEMAX

%----Check inputs & GLOBALS
hdr=[]; seis=[]; reel=[];
if nargin<1 | nargin>8
  error('READSEGY needs 1-8 input arguments')
end
if nargin<2; shotno=0; end;
if nargin<3; tlim=[0 0]; end;
if nargin<4; vred=0;; end;
if nargin<5; rlim=[0 0]; end;
if nargin<6; rpass=0; end;
if nargin<7; decfac=1; end;
decfac=max(1,round(decfac));
if decfac>1
  disp(['  READ_SEGY: Decimating w/ antialias filter - Returning every ',int2str(decfac),' sample'])
  disp(' '); drawnow
end
if ~length(shotno); shotno=0; end;
if length(tlim)~=2 ; tlim=[0 0]; end;
if tlim(1)>tlim(2)
  tlim=[tlim(2) tlim(1)];
  disp('  READ_SEGY: Swapping time limits')
  disp(' '); drawnow
end
if length(vred)~=1; vred=0; end;
if length(rlim)~=2 ; rlim=[0 0]; end;
if rlim(1)>rlim(2)
  rlim=[rlim(2) rlim(1)];
  disp('  READ_SEGY: Swapping range limits')
  disp(' '); drawnow
end
if length(rpass)~=1; rpass=0; end;
if ~length(MONOTONIC); MONOTONIC=0; end
if ~length(RANGESIGNED); RANGESIGNED=1; end
if ~length(DELTARANGEMAX); DELTARANGEMAX=0; end
% EEEH - start
if ~length(chan) ; chan==0;  
elseif (length(chan)>=1 & chan~=0); 
  MONOTONIC=0; 
end
% EEEH - end

%----Check for the existence of the file and get its length
dirinfo=dir(filename);
if ~length(dirinfo)
  disp(['  Problems finding SEGY file ',filename]);
  disp(' '); drawnow
  return
end
filesize=dirinfo.bytes;

%----Open file
little_endian=0;   % Little endian is non-standard
reread=1;
while reread
  reread=0;
  if ~little_endian
    fid=fopen(filename,'r','ieee-be');
  else
    fid=fopen(filename,'r','ieee-le');
  end
  if fid<0
    disp(['  Failed to open SEGY file ',filename]);
    disp(' '); drawnow
    return
  end

%----Skip EBCDIC portion of reel header
  if fseek(fid,3200,'bof')
    disp('  Failed to skip EBCDIC portion of reel header');
    disp(' '); drawnow
    fclose(fid);
    return
  end

%----Read binary portion of reel header
  [in,count(1)]=fread(fid,3,'int32');
  [in2,count(2)]=fread(fid,24,'int16');
  in=[in; in2];
  if sum(count)~=27
    disp('  Failed to read binary portion of reel header');
    disp(' '); drawnow
    fclose(fid);
    return
  end
  reel.job=in(1);
  reel.line=in(2);
  reel.no=in(3);
  reel.trace_per_record=in(4);
  reel.aux=in(5);
  reel.sintr=in(6);
  reel.nsamp=in(8);
  reel.format=in(10);
  if reel.format~=1 & reel.format~=2 & reel.format~=3 & reel.format~=5 & reel.format~=0,
    if reel.format==1280 ,
      little_endian=1;
      reread=1;
      fclose(fid);
      disp(['  READ_SEGY will attempt to read this file as non-standard little endian'])   
      disp(' '); drawnow
    else
      disp(['  READ_SEGY currently cannot read reel format ' int2str(in(10))])
      disp(['  This can be fixed with a little work!'])
      disp(' '); drawnow
      keyboard
    end
  end
end
reel.cdp_fold=in(11);
reel.sort=in(12);
reel.units=in(25);
reel.polarity=in(26);
if reel.format==1 | reel.format==2| reel.format==5 | reel.format ==0
  bytes_per_sample=4;
elseif reel.format==3
  bytes_per_sample=2;
end
bytes_per_trace=240+bytes_per_sample*reel.nsamp;

%----Determine the number of traces
ntrace=(filesize-3600)/bytes_per_trace;
if nargin<2
  shotno=-(1:ntrace);
end

% EEEH - start
%---Determine shot numbers. Look at all traces if not monotonic
if ( any(shotno==0)  & any(chan==0) );
  trace=1:ntrace;
elseif ( any(shotno<0) & any(chan==0) ) ;
  trace=abs(shotno);
else
  if ( any(shotno==0)  & any(chan~=0) );
    shotno=1:ntrace;
  end
  if ( MONOTONIC & any(chan==0) )
    if fseek(fid,3608,'bof');
      disp('  Failed to skip to shot number in trace header 1')
      disp(' '); drawnow
      fclose(fid);
      return
    end
    [shot,count]=fread(fid,1,'int32');
    if sum(count)~=1
      disp('  Failed to read field record number for trace 1');
      disp(' '); drawnow
      fclose(fid);
      return
    end
    shot=(1:ntrace)+shot-1;
  else
    for itrace=1:ntrace
    if fseek(fid,3600+bytes_per_trace*(itrace-1),'bof');
      disp(['  Failed to skip to start of trace header ',itrace])
      disp(' '); drawnow
      fclose(fid);
      return
    end
      [in,count(1)]=fread(fid,7,'int32');
      [in2,count(2)]=fread(fid,4,'int16');
      [in3,count(3)]=fread(fid,8,'int32');
      [in4,count(4)]=fread(fid,2,'int16');
      [in5,count(5)]=fread(fid,4,'int32');
      [in6,count(6)]=fread(fid,76,'int16');
      in=[in; in2; in3; in4; in5; in6];
      if sum(count)~=101
        disp(['  Failed to read header ',int2str(sum(count)),' for trace ',itrace]);
        disp(' '); drawnow
        fclose(fid);
        return
      end
      shot(itrace) = in(3);
%       % For Newberry eeeh start
%       shot(itrace) = in(68);
%       if itrace == 1, disp( ' For Newberry setting shot = hdr.recid'), end
%       % For Newberry eeeh end

%      [shot(itrace),count]=fread(fid,1,'int32');
%       if sum(count)~=1
%         disp(['  Failed to read field record number for trace ',itrace]);
%         disp(' ')
%       disp(' '); drawnow
%         fclose(fid);
%         return
%       end
      if any(chan~=0)
        channel(itrace)=in(4);
 %       [channel(itrace),countchan]=fread(fid,1,'int32');
%         if sum(countchan)~=1
%           disp(['  Failed to read channel field record number for trace ',itrace]);
%           disp(' ')
%           fclose(fid);
%           return
%	        end
      end
    end
  end
  
  %### KLUGE (DTD) #######
  trace=[];  
  % Determine if SIO or WHOI instrument.
  insttype = strfind(filename,'SIO_segy'); % DTD add-on
 % for i=shotno                 % EEH changed for vector of shots  
  for i=1:length(shotno)                    
%     if any(chan~=0)
    if any(chan~=0) && isempty(insttype)==1 % DTD add-on
      for j=chan
    %    trace=[trace find(shot==i & channel==j)];  EEH changed to:
        trace=[trace find(shot==shotno(i) & channel==j)];
      end
    else
    %  trace=[trace find(shot==i)];      % EEEH cahnged to:
      trace=[trace find(shot==shotno(i))];
    end
  end
  %#######################  
end
if ~length(trace)
  disp('  READ_SEGY found no traces with required channels and shot numbers')
  disp(' '); drawnow
  hdr=[];
  seis=[];
  fclose(fid);
  return
end

%----Find traces with ranges within limits (or absolute limits if the 
%                                           ranges are unsigned)
if any(rlim)
  if RANGESIGNED
    rlimuse=rlim*1000;
  else
    rlimuse=[~diff(sign(rlim))*min(abs(rlim)) max(abs(rlim))]*1000;
  end
  if all(diff(trace)==1)
    it1=trace(1); it2=trace(end);
  else
    it1=1; it2=ntrace;
  end
  it=1; ipass=0; lastwithin=0; i1=[]; i2=[];
  while it<=it2 & (rpass==0 | ipass<rpass | (ipass==rpass & lastwithin))
    if fseek(fid,3636+bytes_per_trace*(it-1),'bof');
      disp(['  Failed to skip to range in trace header ',itrace])
      disp(' '); drawnow
      fclose(fid);
      return
    end
    [r,count]=fread(fid,1,'int32');
    if sum(count)~=1
      disp(['  Failed to read range for trace ',itrace]);
      disp(' '); drawnow
      fclose(fid);
      return
    end
    if ~RANGESIGNED; r=abs(r); end 
    within=r>=rlimuse(1) & r<=rlimuse(2);
    if within & ~lastwithin
      ipass=ipass+1;
      i1(ipass)=it;
    elseif ~within & lastwithin
      i2(ipass)=itlast;
    end
    itlast=it;
    lastwithin=within;
    iskip=floor(min(abs(rlimuse-r))/abs(DELTARANGEMAX))*~(~DELTARANGEMAX)+1;
    it=min(it2,it+iskip)+(it==it2);
  end
  if within & length(i1)>length(i2)
    i2(ipass)=itlast;
  end
  index=zeros(size(trace));
  if ~rpass
    for i=1:ipass
      index=index+(trace>=i1(i) & trace<=i2(i));
    end
  elseif rpass<=ipass
    index=(trace>=i1(rpass) & trace<=i2(rpass));
  end
  trace=trace(find(index));
  if ~length(trace)
    disp('  READ_SEGY found no traces with required shot numbers and ranges')
    disp(' '); drawnow
    hdr=[];
    seis=[];
    fclose(fid);
    return
  end
end

%----Loop through traces
i=0; 
fprintf('  Creating SEIS Matrix\n  ')
nread=ceil(reel.nsamp/decfac);
if any(tlim); nread=min(nread,round(diff(tlim)*1e6/reel.sintr/decfac)); end
seis=zeros(nread,length(trace));
for itrace=trace;
  %----Read trace header
  if fseek(fid,3600+bytes_per_trace*(itrace-1),'bof');
    disp(['  Failed to skip to start of trace header ',itrace])
    disp(' '); drawnow
    fclose(fid);
    return
  end
  [in,count(1)]=fread(fid,7,'int32');
  [in2,count(2)]=fread(fid,4,'int16');
  [in3,count(3)]=fread(fid,8,'int32');
  [in4,count(4)]=fread(fid,2,'int16');
  [in5,count(5)]=fread(fid,4,'int32');
  [in6,count(6)]=fread(fid,76,'int16');
  in=[in; in2; in3; in4; in5; in6];
  if sum(count)~=101
    disp(['  Failed to read header ',int2str(sum(count)),' for trace ',itrace]);
    disp(' '); drawnow
    fclose(fid);
    return
  end
  i=i+1;
  %----Assign header values (add your own favorites)
  IN(:,i) = in;
  hdr(i).trace=in(1);
  hdr(i).shot=in(3);
  hdr(i).channel=in(4);
  hdr(i).nstack=in(9);
  hdr(i).range_segyhdr=in(12);
  hdr(i).rec_elev=in(13);
  hdr(i).source_elev=in(14);
  temp=1; if in(20)>0; temp=in(20); elseif in(20)<0; temp=-1/in(20); end
  hdr(i).source_depth=in(15)*temp;
  hdr(i).source_water_depth=in(18)*temp;
  hdr(i).rec_water_depth=in(19)*temp;
  temp=1; if in(21)>0; temp=in(21); elseif in(21)<0; temp=(-1/in(21)); end
  hdr(i).source_long=in(22)*temp;
  hdr(i).source_lat=in(23)*temp;
  hdr(i).rec_long=in(24)*temp;
  hdr(i).rec_lat=in(25)*temp;
  hdr(i).units=in(26);
  hdr(i).delay=in(36);
  hdr(i).nsamp=in(39);
  if hdr(i).nsamp~=reel.nsamp
    error('Reel and trace headers - NSAMP inconsistent')
  end
  nsampdec=ceil(in(39)/decfac);
  hdr(i).sampint=in(40); sampint=in(40)/1e6;
  sampintdec=in(40)*decfac;
  hdr(i).year=in(60);
  hdr(i).day=in(61);
  hdr(i).hour=in(62);
  hdr(i).min=in(63);
  hdr(i).sec=in(64);
  hdr(i).time_basis=in(65);
  hdr(i).bigtime=date2secnds(in(60),in(61),in(62),in(63),in(64)+in(85)/1000);
  hdr(i).recid=in(68);
  hdr(i).siteid=in(69);
  hdr(i).msec=in(85);
  
  if hdr(i).units == 2,
      hdr(i).rec_lat = hdr(i).rec_lat/3600;
      hdr(i).rec_long = hdr(i).rec_long/3600;      
      hdr(i).source_lat = hdr(i).source_lat/3600;
      hdr(i).source_long = hdr(i).source_long/3600;      
  end
  
  %%% *** TO set RANGES correctly for SIO SEGY on 9/19/09  ***
  %%% start
  hdr(i).range = sqrt( ((hdr(i).rec_lat-hdr(i).source_lat)*111)^2 + ...
      ((hdr(i).rec_long - hdr(i).source_long)*111*cos(48*pi/180))^2 )*1000;
  %%% end 
  
  if i==1; hdr(length(trace)).trace=0; end
  %----Determine number of samples to skip and read if time limits set
  if any(tlim)
    if vred
       nskip=max( round((tlim(1) + abs(hdr(i).range/1000/vred) ...
          - hdr(i).delay/1000) / sampint),0);
       % modified by Qin Li  08/10/2000
       hdr(i).delay=max(tlim(1)*1000+round(abs(hdr(i).range/vred)),hdr(i).delay);
%      nread=max( min( floor((tlim(2) + abs(hdr(i).range/1000/vred) ...
%            - hdr(i).delay/1000) / sampint)+1,hdr(i).nsamp)-nskip,0);
    else
       nskip=max( round((tlim(1) - hdr(i).delay/1000) / sampint),0);
       % modified by Qin Li  08/10/2000
       hdr(i).delay=max(tlim(1)*1000,hdr(i).delay);
%      nread=max( min( floor((tlim(2) - hdr(i).delay/1000) ...
%                  / sampint)+1,hdr(i).nsamp)-nskip, 0);
    end
    nread=min(round((tlim(2)-tlim(1))/sampint),reel.nsamp-nskip);
%   hdr(i).delay=hdr(i).delay+round(1000*nskip*sampint);  
    hdr(i).nsamp=nread;
  else
    nskip=0;
    nread=hdr(i).nsamp;
  end
  %----Read in data 
  if nread
    if fseek(fid,nskip*bytes_per_sample,'cof');
      disp(['  Failed to skip to seismogram for trace ',itrace])
      disp(' '); drawnow
      fclose(fid);
      return
    end
% Non-standard SEGY reel format IEEE (big-endian) 4 byte floating
    if reel.format==5
      [seis1,count]=fread(fid,nread,'float32');
      if count~=nread
        disp(['  Failed to read data point ',int2str(count)+1,' for trace ', ...
               int2str(itrace)]);
        disp(' ')
      disp(' '); drawnow
        fclose(fid);
        return
      end
% SEGY reel format = 2 - 4 byte integer
    elseif reel.format==2
      [seis1,count]=fread(fid,nread,'int32');
      if count~=nread
        disp(['  Failed to read data point ',int2str(count)+1,' for trace ', ...
               int2str(itrace)]);
        disp(' '); drawnow
        fclose(fid);
        return
      end
% SEGY reel format = 3 - 2 byte integer
    elseif reel.format==3
      [seis1,count]=fread(fid,nread,'int16');
      if count~=nread
        disp(['  Failed to read data point ',int2str(count)+1,' for trace ', ...
               int2str(itrace)]);
        disp(' '); drawnow
        fclose(fid);
        return
      end
% SEGY reel format = 1 -  non-IEEE (IBM) 32 bit floating point format;
    elseif reel.format==1
      [c,count]=fread(fid,nread*4,'uchar');
      if count~=4*nread
        disp(['  Failed to read data point ',int2str(count)+1,' for trace ', ...
               int2str(itrace)]);
        disp(' '); drawnow
        fclose(fid);
        return
      end
      qf=((c(4:4:nread*4)/256+c(3:4:nread*4))/256+c(2:4:nread*4))/256;
      qc=c(1:4:nread*4);
      qs=(qc<128)*2-1;
      qc=qc-128.*(qs<1)-64;
      seis1=qs.*16 .^(qc) .*qf;   
    end
    if decfac==1
      seis(1:nread,i)=seis1(:);
    else
      seis1=resample(seis1,1,decfac);
      seis(1:ceil(nread/decfac),i)=seis1(:);
      hdr(i).sampint=sampintdec;
      hdr(i).nsamp=nsampdec;
    end
    fprintf('%i ',hdr(i).shot)
  end
end
fclose(fid);
fprintf('\n\n')
disp('  Finished reading traces')
disp(' '); drawnow

%----Eliminate traces with no data inside time limits
if any(tlim);
  blank=~[hdr.nsamp];
  if all(blank)
    disp('  READ_SEGY: None of the traces contain data in time limits')
    disp(' '); drawnow
    seis=[]; hdr=[]; reel=[];
    return
  elseif any(blank);
    disp('  READ_SEGY: Eliminating traces with no data in time limits')
    disp(' '); drawnow
    seis=seis(:,blank); hdr=hdr(blank);
  end
end

%Modify number of samples in reel header
if any(tlim); reel.nsamp=nread; end;

%----Sign ranges for unsigned data and eliminate traces outside signed ranges
if ~RANGESIGNED
  range=[hdr.range]/1e3;
  if any(range<0)
    disp('  READ_SEGY - Note: Range appears to be signed but will be treated as unsigned');
    range=abs(range);
    disp(' '); drawnow
  end
  sequential=diff([hdr.shot])==1;
  dr=diff(range).*sequential;
  rsign=sign([0 dr]+[dr 0]);
  if any(~rsign)
    disp('  READ_SEGY - Warning some ranges may be improperly signed')
    disp(' '); drawnow
    rsign(~rsign)=rsign(~rsign)+1;
  end
%----Remove range sign flips in first and last 5% or record section
%----Often occurs as a result of ships turning at start + end of line
  if DELTARANGEMAX
    drmax=DELTARANGEMAX/1000;
  else
    drmax=2;
  end
  range=range.*rsign;
  dr=diff(range);
  index=find(abs(dr)>drmax);
  fracpos=index/length(dr);
  low=index(find(fracpos<=0.05));
  if length(low)
    low=max(low); 
    rsign(1:low)=ones(size(rsign(1:low)))*rsign(low+1);
    range(1:low)=rsign(1:low).*range(1:low);
    disp('  Fixing range inconsistencies in the first 5% of shots');
    disp(' '); drawnow
  end
  high=index(find(fracpos>=0.95));
  if length(high)
    high=min(high)+1; 
    rsign(high:end)=ones(size(rsign(high:end)))*rsign(high-1);
    range(high:end)=rsign(high:end).*range(high:end);
    disp('  Fixing range inconsistencies in the last 5% of shots');
    disp(' '); drawnow
  end
%----Scrap range signing if still more than one change in sign
  dr=diff(range);
  index=find(abs(dr)>drmax);
  if length(index)>1
    disp('  Geometry too complicated for meaningful range signs - All ranges positive')
    rsign=abs(rsign);
    disp(' '); drawnow
  end
%----Assign range signs and eliminate traces outside range limits
  for i=1:length(hdr)
    hdr(i).range=abs(hdr(i).range)*rsign(i);
  end
  if any(rlim);
    range=range.*(rsign);
    l=range>=rlim(1) & range<=rlim(2);
    if ~any(l)
      disp('  READ_SEGY: No traces within signed range limits')
      disp(' '); drawnow
      seis=[]; hdr=[]; reel=[];
      return
    elseif ~all(l)
      disp('  READ_SEGY: Eliminating traces outside signed range limits')
      disp(' '); drawnow
      seis=seis(:,l);
      hdr=hdr(l);
    end
  end
end
