function [] = fn_writeseg2(Data,smrate,c,ns,H,path)
%     Version 1.1 
%   ` Write seg-2 formated data
%     Matlab version by Andrew Lamb, BSU, August 2011
%     References:
%     Pullan, S. E., Recommended standard for seismic (/radar) data files 
%     in the personal computer environment. Geophysics 55, 1260-1271 (1990)
%     SegyMAT code: http://segymat.sourceforge.net/ 


% Function Input Parameters
% Data = recorded data for each trace in the SegyMAT format
% smrate = sample rate
% c = number of channels/traces 
% ns = number of samples per trace
% H=SegyTraceHeaders structure created by SegyMAT's ReadSegy command
% path = path to write seg2 file according to definition by Pullan 1990
%
% Function Output
% seg2 file with '[H(1).EnergySourcePoint].sg2' as the filename
%

%% SEG2 file format schematic (example)
% Schematic is based on example using c=48, ns=150
% and strings with * as recommended by Pullan 1990.
% These strings are 'TRACE_SORT AS_ACQUIRED' and 'UNITS METERS' for the 
% File Descriptor Block (total 44 bytes) and DELAY, RECEIVER_LOCATION and 
% SAMPLE_INTERVAL for the Trace Descriptor Block (total 64 bytes)
% ---------------------------------------
% |*******FILE DESCRIPTOR BLOCK*********|
% |           MAIN HEADERS              |
% |         bytes 1 - 16 (16)           |
% |             RESERVED                |
% |             17 - 32 (16)            |
% |         TRACE POINTER SUB-BLOCK     |
% |             33 - 224 (192)          |
% |  assumes 4 bytes/trace & 48 traces  |
% |               STRINGS               |
% |             225 - 268 (44)          |
% |             [TOTAL=268]             |
% ---------------------------------------
% ---------------------------------------
% |******TRACE DESCRIPTOR BLOCK 1*******|
% |           MAIN HEADERS              |
% |         bytes 269 -284 (16)         |
% |             RESERVED                |
% |            285 - 300 (16)           |
% |               STRINGS               |
% |             301 - 364 (64)          |
% |            [SUBTOTAL=96]            |
% ---------------------------------------
% ---------------------------------------
% |************DATA BLOCK 1*************|
% |           365 - 964 (600)           |
% |            [SUBTOTAL=600]           |
% ---------------------------------------
% |       [BLOCK 1 TOTAL=696]           |
% |||||||||||||||||||||||||||||||||||||||
% TRACE DESCRIPTOR AND DATA BLOCKS 2 - 47
% |||||||||||||||||||||||||||||||||||||||
% ---------------------------------------
% |******TRACE DESCRIPTOR BLOCK 48******|
% |           MAIN HEADERS              |
% |      bytes 32981 - 32996 (16)       |
% |             RESERVED                |
% |          32997 - 33012 (16)         |
% |               STRINGS               |
% |          33013 - 33076 (64)         |
% |            [SUBTOTAL=96]            |
% ---------------------------------------
% ---------------------------------------
% |************DATA BLOCK 48************|
% |         33077 - 33676 (600)         |
% |            [SUBTOTAL=600]           |
% ---------------------------------------
% FILE SIZE = 33,676 bytes (as written to disk)

fid=fopen(strcat(path,num2str([H(1).EnergySourcePoint]),'.sg2'),'w');
i=1;

%% Prepare Strings for File Descriptor Block
fds1='TRACE_SORT AS_ACQUIRED';                    %File Descriptor String 1
fds2='UNITS METERS';                              %File Descriptor String 2

%Next statement calculates the combined length of strings including the two
%empty bytes required at the end of the string sub-block to mark the end
fdslen=length(fds1)+3+length(fds2)+3+2; 

%Next line checks how many bytes are needed to pad block to a division of 4
fdspad=4-((fdslen)/4-floor((fdslen)/4))/0.25;

% Next statement is needed because if tspad==4 then there is no need to pad
if fdspad==4 
    fdspad=0; 
end

%% Prepare Strings for Trace Descriptor Block
% This section estimates the lengths of each string in the Trace
% Descriptor Block and assigns an initial value
tds{1}=['DELAY ' sprintf('%8.4f',[H(i).DelayRecordingTime])];
tds{2}=['RECEIVER_LOCATION ' sprintf('%15.3f',[H(i).GroupX])];
tds{3}=['SAMPLE_INTERVAL ' sprintf('%0.4f',smrate)];
tds{4}=['CDP_NUMBER ' sprintf('%08d',[H(i).cdp])];
tds{5}=['CDP_TRACE ' sprintf('%08d',[H(i).cdpTrace])];
% I've multiplied the next line by zero to force the receiver inline offset
% to be zero for reading into Rayfract    
tds{6}=['RECEIVER_GEOMETRY ' sprintf('%15.3f',[H(i).GroupX])];
tds{7}=['SHOT_SEQUENCE_NUMBER ' sprintf('%08d',[H(i).EnergySourcePoint])]; %Rayfract: 'File number'
tds{8}=['SOURCE_LOCATION ' sprintf('%15.3f',[H(i).EnergySourcePoint])];
% I've multiplied the next line by zero to force the shot inline offset
% to be zero for reading into Rayfract
tds{9}=['SOURCE_GEOMETRY ' sprintf('%15.3f',[H(i).SourceX]*0)]; %Rayfract: 'Shot inline offset'
tds{10}=['SOURCE_STATION_NUMBER ' sprintf('%08d',[H(i).SourceX])]; %Rayfract: 'Shot pos. [Station no.]'
tds{11}=['CHANNEL_NUMBER ' sprintf('%08d',[H(i).cdpTrace])]; 
tds{12}=['RAW_RECORD ' sprintf('%08d',[H(i).FieldRecord])]; %Rayfract: 'Record number'
tds{13}=['SOURCE VIBROSEIS'];

%Next statement calculates the combined length of strings including the two
%empty bytes required at the end of the string sub-block to mark the end
tdslen=sum(cellfun(@length,tds))+3*size(tds,2)+2;

%Next line checks how many bytes are needed to pad block to a division of 4
tdspad=4-((tdslen)/4-floor((tdslen)/4))/0.25; 

% Next statement is needed because if tspad==4 then there is no need to pad
if tdspad==4 
    tdspad=0; 
end

%% Write FILE DESCRIPTOR BLOCK
% This section writes the File Descriptor Block using the definitions from
% the SEG2 standard definition (Pullan 1990).
% Write Main Headers
%fwrite(fid, A, precision, skip, machineformat)
fwrite(fid, hex2dec('3A55'), 'uint16', 0, 'n'); %File descriptor 0-1 55 1st => PC
fwrite(fid, 1, 'uint16', 0, 'n'); %Revision number. 2-3 int16
fwrite(fid, c*4, 'uint16', 0, 'n'); %Size of Trace Pointer SubBlock. 4-5
fwrite(fid, c, 'uint16', 0, 'n'); %No. of traces in this file 6-7
fwrite(fid, hex2dec('01'), 'uint8', 0, 'n'); %Size of string term. 8
fwrite(fid, 0, 'uint8', 0, 'n'); %Null. 9 int8
fwrite(fid, hex2dec('02'), 'uint8', 1, 'n'); %Size of string term. 11
fwrite(fid, hex2dec('0A'), 'uint8', 0, 'n'); %1st line of terminator (LF). 12
fwrite(fid, hex2dec('0D'), 'uint8', 0, 'n'); %1st line of terminator (CR). 13

%Write trace pointer sub-blocks. These give byte location of trace = i
for i =1:c
    if i==1 % When c==1 the pointer skips over the RESERVED Bytes
        fwrite(fid, 32+(4*(c))+fdslen+fdspad+((32+tdslen+tdspad+(ns*4))*...
           (i-1)), 'ulong', 18, 'n');
    else  %Writes byte location of remaining trace pointer sub-blocks
        fwrite(fid, 32+(4*(c))+fdslen+fdspad+((32+tdslen+tdspad+(ns*4))*...
            (i-1)), 'ulong', 0, 'n');
    end
end

% Write Strings for File Descriptor block
% Write TRACE_SORT string
fwrite(fid, length(fds1)+3, 'uint16' ,0,'n');
fwrite(fid, fds1,'char',0,'n');
fwrite(fid, 0, 'uint8', 0,'n');
% Write UNITS string
fwrite(fid, length(fds2)+3, 'uint16' ,0,'n');
fwrite(fid, fds2,'char',0,'n');
fwrite(fid, 0, 'uint8', 0,'n');

%% Write TRACE DESCRIPTOR and DATA BLOCK loops
% Write Main Headers
for i=1:c
    % Reassign trace descriptor strings with i value
    tds{1}=['DELAY ' sprintf('%8.4f',[H(i).DelayRecordingTime])];
    tds{2}=['RECEIVER_LOCATION ' sprintf('%15.3f',[H(i).GroupX])];
    tds{3}=['SAMPLE_INTERVAL ' sprintf('%0.4f',smrate)];
    tds{4}=['CDP_NUMBER ' sprintf('%08d',[H(i).cdp])];
    tds{5}=['CDP_TRACE ' sprintf('%08d',[H(i).cdpTrace])];
    % I've multiplied the next line by zero to force the receiver inline offset
    % to be zero for reading into Rayfract    
    tds{6}=['RECEIVER_GEOMETRY ' sprintf('%15.3f',[H(i).GroupX])];
    tds{7}=['SHOT_SEQUENCE_NUMBER ' sprintf('%08d',[H(i).EnergySourcePoint])]; %Rayfract: 'File number'
    tds{8}=['SOURCE_LOCATION ' sprintf('%15.3f',[H(i).EnergySourcePoint])];
    % I've multiplied the next line by zero to force the shot inline offset
    % to be zero for reading into Rayfract
    tds{9}=['SOURCE_GEOMETRY ' sprintf('%15.3f',[H(i).SourceX]*0)]; %Rayfract: 'Shot inline offset'
    tds{10}=['SOURCE_STATION_NUMBER ' sprintf('%08d',[H(i).SourceX])]; %Rayfract: 'Shot pos. [Station no.]'
    tds{11}=['CHANNEL_NUMBER ' sprintf('%08d',[H(i).cdpTrace])]; 
    tds{12}=['RAW_RECORD ' sprintf('%08d',[H(i).FieldRecord])]; %Rayfract: 'Record number'
    tds{13}=['SOURCE VIBROSEIS'];

    if i==1
        skip=2+fdspad; %Number of bytes to skip to ensure file descriptor
                       %block is a division of 4
    else
        skip=0; %No skipping is required one first byte of Trace descriptor
                %block has been written
    end
    
    fwrite(fid, hex2dec('4422'), 'uint16', skip, 'n'); %Trace descriptor 0-1
    fwrite(fid, 32+tdslen+tdspad, 'uint16', 0, 'n'); %X=Size of this block 2-3
    fwrite(fid, ns*4, 'ulong', 0, 'n'); %Y=Size of data block 4-7
    fwrite(fid, ns, 'ulong', 0, 'n'); %NS=Number of samples in data block 8-11
    fwrite(fid, hex2dec('04'), 'uint8', 0, 'n'); %Data format
    
    %Write Strings for Trace Descriptor block
    for k=1:size(tds,2)
        if k==1
            skip=19; %This skips over the 20 reserve bytes before writing strings
        else
            skip=0;
        end
        fwrite(fid, length(tds{k})+3, 'uint16' ,skip,'n');% Skip 20 reserved bytes
        fwrite(fid, tds{k}, 'char',0,'n');
        fwrite(fid, 0, 'uint8', 0,'n');
    end
%     % Write DELAY string
%     fwrite(fid, length(tds1)+3, 'uint16' ,19,'n');% Skip 20 reserved bytes
%     fwrite(fid, tds1, 'char',0,'n');
%     fwrite(fid, 0, 'uint8', 0,'n');
%     % Write RECEIVER_LOCATION string
%     fwrite(fid, length(tds2)+3, 'uint16' ,0,'n');
%     fwrite(fid, tds2, 'char',0,'n');
%     fwrite(fid, 0, 'uint8', 0,'n');
%     % Write SAMPLE_INTERVAL string
%     fwrite(fid, length(tds3)+3, 'uint16' ,0,'n');
%     fwrite(fid, tds3, 'char',0,'n');
%     fwrite(fid, 0, 'uint8', 0,'n');
%     % Write CDP_NUMBER string
%     fwrite(fid, length(tds4)+3, 'uint16' ,0,'n');
%     fwrite(fid, tds4, 'char',0,'n');
%     fwrite(fid, 0, 'uint8', 0,'n');
    
   % DATA BLOCK loop
    for j=1:ns
        if j==1
            fwrite(fid, Data(j,i), 'float32', 2+tdspad, 'n');
        else
            fwrite(fid, Data(j,i), 'float32', 0, 'n');
        end
    end
end
fclose(fid);