open_segy_file4reading.m

实现地震勘探中

   end
   param.precision='uint32=>uint32';

elseif binary_header(10) == 5
   if strcmpi(param.format,'header')  || strcmpi(param.format,'ieee')
%      param.fmt='ieee-be'; 
      param.format='ieee-be';
   else
      display('IBM format requested')
      mywarning(['Data apparently stored as 32-bit IEEE big-endian floating point numbers; ' ...
         'data sample format code = ',num2str(binary_header(10))]);      
   end
   param.precision='single';

else
   if param.debug
      disp(['Data in unsupported format; ' ...
            'data sample format code = ',num2str(binary_header(10))]);
      param.precision='single';

   else
      myerror(['Data in unsupported format; ' ...
            'data sample format code = ',num2str(binary_header(10))]);
      ABORTED=true;
      fclose(fid)
      return
   end
end

%	Extract units of measurement for distances
if binary_header(25) == 1
   param.units='m';
elseif binary_header(25) == 2
   param.units='ft';
else
   param.units='unknown';
end

ABORTED=false;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [seismic,idx4times]=create_seismic_structure_no5(binary_header,param)
% Create the fields of the seismic structure and populate them as far as possible
% INPUT
% binary_header   binary header
% param           structure with input parameters for "read_segy_file"
% OUTPUT
% seismic         seismic structure
% idx4time        index of time samples to keep

global S4M

step=binary_header(6)/1000;
no_samples=binary_header(8);
last=step*(no_samples-1);


%	Compute indices of time range 
if ~isempty(param.times)    % Compute index vector for time samples
   ita=ceil(param.times(1)/step);
   param.times(1)=ita*step;
   
   param.times(2)=min(last,param.times(2));
   ite=fix(param.times(2)/step);
   param.times(2)=ite*step;

   idx4times=ita+1:ite+1;

else
  param.times=[0,last];
  idx4times=(1:no_samples)';   % Keep all time samples

end

%	Create seismic structure
seismic.type='seismic';
seismic.tag='unspecified';
[dummy,seismic.name]=fileparts(S4M.filename);  %#ok First output variable is not required
seismic.from=fullfile(S4M.pathname,S4M.filename);
seismic.line_number=binary_header(2);
seismic.reel_number=binary_header(3);
seismic.traces_per_record=binary_header(4);
seismic.aux_per_record=binary_header(5);
if binary_header(11) ~= 0
  seismic.cdp_fold=binary_header(11);
end

seismic.first=param.times(1);
seismic.last=param.times(2);
seismic.step=step;
seismic.units='ms';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [header_info,indices,true4four,constraint_info,param]=select_trace_headers2read_no6(param)
% Create a cell array with ibformation about the headers that are to be read
% their bit location and length
% INPUT
% param   input parameters of read_segy_file (plus ome parameters 
%         acquired along the way
% OUTPUT
% header_info
%	The 7 columns of header_info are:
%       mnemonic, units, description, first byte, No. of bytes,
%       bool ("true" if user-requested), and index.
%       If the 240-byte header is read into a 60  element array then the index 
%       defines in which element a 4-byte header is: for a 4-byte header word
%       starting at byte 13 index=4 (index=(start_byte-1)/4+1)
%       The index for a two-byte header is similarly defined.
% indices    vector of header indices
% true4four  boolean vector true if index is for a four-byte header.

global S4M

constraint_info=[];

if isempty(param.headers)
   nh1=0;
elseif iscell(param.headers{1})
   nh1=length(param.headers);  % User-requested headers
else
   param.headers={param.headers};
   nh1=1;
end

%         Select standard headers and add user-defined headers
%	Descriptions for standard binary trace headers
nh0=22;		% Default headers
nh=nh0+nh1;     % Total number of headers
header_info=cell(nh,7);

header_info(1:nh0,6)=deal({false});

header_info(1:nh0,3)=[...
      {'Trace sequence number within line'}; ...        % 4-byte headers 
      {'Original Field record number'}; ...
      {'Trace sequence number within original field record'}; ...
      {'Energy source point number'}; ...
      {'CDP number'}; ...
      {'Trace sequence number within CDP ensemble'}; ...
      {'Offset'}; ...
      {'Source depth below surface'}; ...
      {'Water depth at source'}; ...
      {'Water depth at receiver group'}; ...
      {'X coordinate of source'}; ...
      {'Y coordinate of source'}; ...
      {'Sourface elevation at source'}; ...
      {'X coordinate of receiver'}; ...
      {'Y coordinate of receiver'}; ...
      {'Receiver elevation'}; ...
      {'X-coordinate of CDP'}; ...
      {'Y-coordinate of CDP'}; ...
      {'In-line number'}; ...
      {'Cross-line number'}; ...
      {'Trace type (1=live,2=dead,3=dummy,...)'}; ...   % 2-byte headers
      {'Lag time between shot and recording start'}; ...
];

header_info(1:nh0,[1,4,5]) =  { ...  % 4-byte headers
      'ds_seqno',      1,4; ... 
      'ffid',          9,4; ...
      'o_trace_no',   13,4; ...  
      'source',       17,4; ...   
      'cdp',          21,4; ...      
      'seq_cdp',      25,4; ... 
      'offset',       37,4; ...     
      'depth',        49,4; ...  
      'sou_h2od',     61,4; ...  
      'rec_h2od',     65,4; ... 
      'sou_x',        73,4; ...     
      'sou_y',        77,4; ...      
      'sou_elev',     45,4; ...      
      'rec_x',        81,4; ...       
      'rec_y',        85,4; ...     
      'rec_elev',     41,4; ...
      'cdp_x',       181,4; ...
      'cdp_y',       185,4; ...
      'iline_no',    189,4; ...
      'xline_no',    193,4; ...
      'trc_type',     29,2; ...  % 2-byte headers
      'lag',         109,2};

header_info(1:nh0,2)=deal({'n/a'});
header_info(7:16,2)={param.units};
header_info(22,2)={'ms'};

if nh1 > 0

%   nh=length(param.headers);
%   nh=nh1;

%      Defaults for some commonly requested headers
   defaults.iline_no={'n/a','In-line number'};
   defaults.xline_no={'n/a','Cross-line number'};
   defaults.cdp_x={param.units,'X-coordinate of CDP'}; 
   defaults.cdp_y={param.units,'Y-coordinate of CDP'}; 
   defaults.rec_x={param.units,'X-coordinate of receiver'}; 
   defaults.rec_y={param.units,'Y-coordinate of receiver'}; 
   defaults.sou_x={param.units,'X-coordinate of source'}; 
   defaults.sou_y={param.units,'Y-coordinate of source'}; 

%	Handle header supplided by users

   header_info(nh0+1:nh,6)=deal({true});

   for ii=nh0+1:nh
      nhii=length(param.headers{ii-nh0});
      if nhii == 5
         header_info(ii,[1,4,5,2,3])=param.headers{ii-nh0};
         if isempty(header_info{ii,2}); 
            header_info{ii,2}='n/a';
         end
         if isempty(header_info{ii,3}); 
            header_info{ii,3}=header_info{ii,1};
         end

      elseif nhii == 3
         header_info(ii,[1,4,5])=param.headers{ii-nh0};
         if any(ismember(fieldnames(defaults),lower(header_info{ii,1})))
            header_info(ii,2:3)=defaults.(lower(header_info{ii,1}));
         else
            header_info{ii,2}='n/a';
       	    header_info{ii,3}=header_info{ii,1};
         end

      else
         error('The description of an additional header to be read must be a cell vector with 3 or 5 entries.') 
      end
   end
   ierr=false;
   for jj=nh0+1:nh
      if mod(header_info{jj,4},header_info{jj,5}) ~= 1
         disp(['First byte for header ',header_info{jj,1}, ' is wrong.'])
         ierr=true;
      end
   end

   if ierr
      error('Abnormal termination.')
   end
end

%	Keep unique headers
[dummy,index]=myunique(header_info(:,1)); %#ok First output argument is not required
header_info=header_info(index,:);
nh=length(index);

indices=zeros(nh,1);
true4four=true(nh,1);
for ii=1:nh
   indices(ii)=(header_info{ii,4}-1)/header_info{ii,5}+1;
   header_info{ii,7}=indices(ii); 
   true4four(ii)=header_info{ii,5} == 4;
end

%	If traces are selected via constraints on headers then 
%	create a cell array with header names and associate locations
%       in the "headers" matrix
if ~isempty(param.traces)  &&  ischar(param.traces)
%	Find headers in constraint string
   expression=param.traces;
   words=symvar(expression);

   if ~S4M.case_sensitive
     %		Change expression
     for ii=1:length(words)
         expression=strrep(expression,words{ii},lower(words{ii}));
      end
      param.traces=expression;
      words=lower(words);
      headers=lower(header_info(:,1));
    else
      headers=header_info(:,1);
   end
   words=unique(words);
   index=find(ismember(headers,words));
   bool=ismember('trace_no',words);
   if ~bool  && isempty(index)
      disp([' No headers found in trace constraint "',expression,'".'])
      myerror(' Reading of SEG-Y file terminated abnormally.')
   else
      constraint_info=cell(length(index)+bool,2);
      for ii=1:length(index)
         constraint_info(ii,:)=[headers(index(ii)),{num2str(index(ii))}];
      end
      if bool
         constraint_info(end,:)={'trace_no',0};
      end
   end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ntr=get_no_of_traces_no7(filename,nsamples)
% Use number of bytes to compute number of traces

ll=dir(filename);
nbytes=ll.bytes;
ntr=0.25*(nbytes-3600)/(nsamples+60);

if ~isnearinteger(ntr)
   mywarning(['Number of bytes in file "',filename,'", (',num2str(nbytes), ...
      '), is not compatible with constant-length traces'])
   ntr=fix(ntr);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ascii=ebcdic2ascii(ebcdic)
% Function converts EBCDIC string to ASCII
% see http://www.room42.com/store/computer_center/code_tables.shtml
%
% Written by: E. R.: Feb. 20, 2000
% Last updated:
%
%           ascii=ebcdic2ascii(ebcdic)
% INPUT
% ebcdic    EBCDIC string
% OUTPUT
% ascii	   ASCII string


pointer= ...
[ 0    16    32    46    32    38    45    46    46    46    46    46   123   125    92    48
  1    17    33    46    46    46    47    46    97   106   126    46    65    74    46    49
  2    18    34    50    46    46    46    46    98   107   115    46    66    75    83    50
  3    19    35    51    46    46    46    46    99   108   116    46    67    76    84    51
  4    20    36    52    46    46    46    46   100   109   117    46    68    77    85    52
  5    21    37    53    46    46    46    46   101   110   118    46    69    78    86    53
  6    22    38    54    46    46    46    46   102   111   119    46    70    79    87    54
  7    23    39    55    46    46    46    46   103   112   120    46    71    80    88    55
  8    24    40    56    46    46    46    46   104   113   121    46    72    81    89    56
  9    25    41    57    46    46    46    46   105   114   122    46    73    82    90    57
 10    26    42    58    46    33   124    58    46    46    46    46    46    46    46    46
 11    27    43    59    46    36    44    35    46    46    46    46    46    46    46    46
 12    28    44    60    60    42    37    64    46    46    46    46    46    46    46    46
 13    29    45    61    40    41    95    39    46    46    91    93    46    46    46    46
 14    30    46    46    43    59    62    61    46    46    46    46    46    46    46    46
 15    31    47    63   124    94    63    34    46    46    46    46    46    46    46    46];

pointer=reshape(pointer,1,256);

ascii=pointer(ebcdic+1);