sgywrite.m

这是matlab在地球物理数据处理方面的源码

function fid = sgywrite(sgyfile,traces,t,offset)
% Program to automatically write  a SEGY file.
% Program returns matrix traces with amplitudes
% and time base t.
% Can easily add more variables if required
% Refer to SEG digital tape standards booklet.
% Written by D. Schmitt, January 1, 1999
% The input file must have the full filename including
% the standard .sgy designation at the end.
% To run [fileid] = sgywrite('filename',tracematrix,timebase,offset);
% Trace and Binary Headers are described. No Ascii Header is 
% written but just filled.


fclose('all');
if length(sgyfile) > 12
	error('Input filename is too long, please rewrite')
end
if sgyfile((length(sgyfile)-3):length(sgyfile)) ~= ['.sgy'] %This doesn't work
	error('Input filename must have .sgy appended')
end	

% fid = fopen(sgyfile,'w','b')  % for workstation mode
fid = fopen(sgyfile,'w','l');  % for standard PC Vista mode

[m,n] = size(traces);

filler = ['0'];
for i = 1:(3200)  % 3200 byte ascii header
	count = fwrite(fid,filler,'char');  % 3200
end
 

% Binary header information, 400 bytes

binhead = [ 1  	% jobid
         1  	% line number *
	 1  	% reel number *
	n	% Number of Data Traces *
 	0 	% Number of Auxiliary Traces *
	(t(2)-t(1))*1e6 	% Sample interval in microseconds *
 	(t(2)-t(1))*1e6 	% Original sample interval in microseconds
 	m 	% Number of samples per trace *
	m 	% Original number of saples per trace
	1 	% Data Sample Format 1=Float(4), 2=Fixed(4),3 = Fixed(2)*
	1	% CDP Fold (expected per CDP ensemble)
 	1	% Trace sorting code, 1 = as recorded
 	10	% Vertical Sum Code (Stack?), N = number of vertical sums. 
 	0	% Sweep Frequency at Start
	0	% Sweep Frequency at End
	0	% Sweep Length (seconds)
	0	% Sweep Type Code
	0	% Trace Number of Sweep Channel
	0	% Sweep Taper Length at start (in ms)
	0	% Sweep Taper Length at End (in ms)
	0	% Taper Type
	0	% Correlated Data Traces 1 = No 2 = Yes
	0	% Binary Gain Recovered 1 = Yes  2 = No
	0	% Amplitude Recovery Method 1=none, 2=spherical, 3=AGC, 4=other
	1	% Measurement system 1= meters, 2 = feet
	1	% Polarity 1 = increase P or upward case motion gives negative
	0];	% Vibratory polarity code - see SEG specs.
binhead(28:197) = zeros(170,1);  % Fill out remainer of binary header

count = fwrite(fid,binhead(1:3),'int32');  %  The first three are 4 byte
count = fwrite(fid,binhead(4:197),'int16');  % Remainder are in 2 byte form 

% Make traceheader information

trahead = zeros(147,n);  % Vista SEGY dictionary has 147 items

trahead(1,:) = 1:n;  	% Trace Sequence Number *
trahead(2,:) = 1:n;  	% Reel trace sequence number *
trahead(3,:) = ones(1,n);    	% Orignal Field Record Number *
trahead(4,:) = 1:n;  	% Trace Number within original field record *
trahead(5,:) = ones(1,n);	% Energy Source Number (Shot Point Number)
trahead(6,:) = ones(1,n);	% CDP Ensemble Number
trahead(7,:) = 1:n;		% Trace number within CDP ensemble
trahead(8,:) = ones(1,n); 	% Trace Identification * see tech info
trahead(9,:) = 10*ones(1,n);	% Number of vertically summed traces (shot stack)
trahead(10,:) = ones(1,n);	% Number of horizontally summed traces	
trahead(11,:) = 2*ones(1,n);	%Data use, 1 = production, 2 = test
trahead(12,:) = offset; 	% Trace-Shot Point Offset in meters
trahead(13,:) = zeros(1,n);	% Receiver Elevation (properly from sea level)
trahead(14,:) = zeros(1,n);	% Source Elevation at surface
trahead(15,:) = zeros(1,n);	% Source Depth (a positive number)
trahead(16,:) = zeros(1,n);	% Receiver Datum elevation
trahead(17,:) = zeros(1,n); 	% Source Datum elevation
trahead(18,:) = zeros(1,n);	% Water Depth at Source
trahead(19,:) = zeros(1,n);	% Water Depth at Group
trahead(20,:) = ones(1,n);	% Scaler to multiply elevations and depths above
trahead(21,:) = ones(1,n);	% Scaler to multiply co-ordinates below
				% Note, scalars are plus/minus 1, 10, 100, ...
				% Positive means multiply, negative means divide
trahead(22,:) = zeros(1,n);	% X source co-ordinate
trahead(23,:) = zeros(1,n);	% Y Source co-ordinate
trahead(24,:) = offset;		% X receiver co-ordinate
trahead(25,:) = zeros(1,n);	% Y receiver co-ordinate
trahead(26,:) = ones(1,n);	% Co-ordinate units, 1(meter or feet),2(seconds of arc)
trahead(27,:) = 300*ones(1,n);	% Weatering velocity
trahead(28,:) = 2000*ones(1,n);	% Subweathering velocity
trahead(29,:) = zeros(1,n);	% Uphole time at source
trahead(30,:) = zeros(1,n);	% Uphole time at receiver
trahead(31,:) = zeros(1,n);	% Source static correction
trahead(32,:) = zeros(1,n);	% Reciever static correction
trahead(33,:) = zeros(1,n);	% Total static applied (zero if none)
trahead(34,:) = zeros(1,n);	% Lag time A in ms
trahead(35,:) = zeros(1,n);	% Lag time B in ms
trahead(36,:) = zeros(1,n);	% Delay Recording time in ms
trahead(37,:) = zeros(1,n); 	% Mute time Start
trahead(38,:) = zeros(1,n);	% Mute time end
trahead(39,:) = m*ones(1,n);	% Number of samples in trace *
trahead(40,:) = (t(2)-t(1))*1e6*ones(1,n); % Sample interval in microseconds for this trace
trahead(41,:) = ones(1,n);	% Gain type of field instruments, see tech standards
trahead(42,:) = ones(1,n);	% Instrument gain constant
trahead(43,:) = ones(1,n);	% Instrument early gain
trahead(44,:) = ones(1,n);	% Correlated 1 = no, 2 = yes
trahead(45,:) = zeros(1,n);	% Sweep freq at start
trahead(46,:) = zeros(1,n);	% Sweep freq at end
trahead(47,:) = zeros(1,n);	% Sweep length
trahead(48,:) = zeros(1,n);	% Sweep type, see tech standards
trahead(49,:) = zeros(1,n);	% Sweep taper at start in ms
trahead(50,:) = zeros(1,n);	% Sweep taper at end in ms
trahead(51,:) = ones(1,n);	% Sweep taper type
trahead(52,:) = zeros(1,n);	% Alias filter freqency if used
trahead(53,:) = zeros(1,n);	% Alias filter slope if used
trahead(54,:) = zeros(1,n);	% Notch filter frequency if used
trahead(55,:) = zeros(1,n);	% Notch filter slope
trahead(56,:) = zeros(1,n);	% Low cut if used
trahead(57,:) = zeros(1,n);	% High Cut if used
trahead(58,:) = zeros(1,n);	% Low cut slope
trahead(59,:) = zeros(1,n);	% High cut slope 
trahead(60,:) = 1998*ones(1,n); % Data recorded year
trahead(61,:) = 925*ones(1,n);	% Day of year 
trahead(62,:) = ones(1,n);	% Hour of day
trahead(63,:) = ones(1,n);	% Minute of hour
trahead(64,:) = ones(1,n);	% Second of minute
trahead(65,:) = ones(1,n);	% Time Basis code, 1 = local, 2 = GMT, 3 = other
trahead(66,:) = zeros(1,n);	% Trace weighting factor - see tech standards
trahead(67,:) = offset;		% Group number on roll switch position 1
trahead(68,:) = zeros(1,n);	% Group number of trace number one within orig field 
trahead(69,:) = ones(1,n);	% Group number of last trace within original field
trahead(70,:) = zeros(1,n);	% Gap Size (total number of groups dropped)
trahead(71,:) = zeros(1,n);	% Overtravel associated with taper 1=behind(down), 2=ahead(up)
%  Note - Normal SEG-Y format stops here, bytes 181-240 are unassigned.



% Loop through to write traces and trace headers.
for i = 1:60
	fil(i) = filler;
end % i - this makes a filler for the end of the binary trace header.

for i = 1:n

	count = fwrite(fid,trahead(1:7,i),'int32');  % Bytes 1-28
	count = fwrite(fid,trahead(8:11,i),'int16'); % Bytes 29-36
	count = fwrite(fid,trahead(12:19,i),'int32'); % Bytes 37-68
	count = fwrite(fid,trahead(20:21,i),'int16'); % Bytes 69-72
	count = fwrite(fid,trahead(22:25,i),'int32'); % Bytes 73-88
	count = fwrite(fid,trahead(26:71,i),'int16'); % Bytes 89-180
	count = fwrite(fid,fil,'char'); % Fill bytes 181-240 

	fwrite(fid,traces(:,i),'float'); %  Write trace data
end
fclose(fid);