readgreen.m

水声模型 很不错的东西

function [SourceDepth, ReceiverDepth, Freq, Kr, GKr, Title] = ReadGreen(FileName)
%[SourceDepth, ReceiverDepth, Freq, Kr, GKr, Title] = ReadGreen(FileName)
%
%Function to read in a Greens function file generated by Scooter
%This is in binary format, so this may not work on all platforms
%
%FileName - full filename, including path, of greens function file
%
%SourceDepth - row vector of source depths (m)
%ReceiverDepth - row vector of receiver depths (m)
%Freq - frequency (Hz)
%Kr row vector of horizontal wavenumbers
%GKr - Greens function.  This is of dimension GKr(NKr, NReceiverDepths, NSourceDepths).
%Title - title read from Greens function file
%
%Alec Duncan,
%Centre for Marine science and Technology
%Curtin University
%Western Australia
%May 2002

Debug = 0;

Real = 'float32';
Char = 'uchar';
Int = 'int32';

InFile = fopen(FileName, 'r');
if InFile < 0
    uiwait(warndlg(['Couldn''t open file: ' FileName]));
else
    RecordLength = fread(InFile, 1, Int);
    Title = deblank(char(fread(InFile, [1, 80], Char)));

    lSeekRecord(InFile, 2, RecordLength);
    PltType = deblank(char(fread(InFile, [1, 10], Char)));
    
    Dat = fread(InFile, [1,3], Real);  %Xs, Ys, Theta
    XsYsTheta = Dat;
    
    lSeekRecord(InFile, 3, RecordLength);
    Freq = fread(InFile, 1, Real);
    Dat = fread(InFile, [1,4], Int);
    NSourceDepths = Dat(1);
    NReceiverDepths = Dat(2);
    NKPts = Dat(3);
    NKCore = Dat(4);
    
    Dat = fread(InFile, [1,2], Real);
    DeltaK = Dat(1);
    Atten = Dat(2);
    
    NRecSr = fread(InFile, 1, Int);  %Number of records for each Greens fn

    lSeekRecord(InFile, 4, RecordLength);
    SourceDepth = fread(InFile, [1, NSourceDepths], Real);
    
    lSeekRecord(InFile, 5, RecordLength);
    ReceiverDepth = fread(InFile, [1, NReceiverDepths], Real);

    if Debug
        DebugStr = {' ', ['RecordLength = ', num2str(RecordLength)], ...
                ['Title = ' Title], ...
                ['PltType = ' PltType], ...
                ['Xs, Ys, Theta = ' num2str(XsYsTheta)], ...
                ['Freq = ' num2str(Freq)], ...
                ['NSourceDepths = ' num2str(NSourceDepths)], ...
                ['NReceiverDepths = ', num2str(NReceiverDepths)], ...
                ['NKPts = ', num2str(NKPts)], ...
                ['NKCore = ', num2str(NKCore)], ...
                ['DeltaK = ', num2str(DeltaK)], ...
                ['Atten = ', num2str(Atten)], ...
                ['NRecSr = ', num2str(NRecSr)], ...
                ['SourceDepth = ', num2str(SourceDepth)], ...
                ['ReceiverDepth = ', num2str(ReceiverDepth)]};
        disp(char(DebugStr));
    end
    Kr = zeros(1, NKPts);
    

    for IRec = 1:NRecSr
        IStart = (IRec - 1) * NKCore + 1;
        IEnd = min([IStart + NKCore - 1, NKPts]);
        NRead = IEnd - IStart + 1;
        
        lSeekRecord(InFile, IRec+5, RecordLength);
        Kr(1, IStart:IEnd)  = fread(InFile, [1, NRead], Real);
    end
    

    GKr = zeros(NKPts, NReceiverDepths, NSourceDepths);
    
    %Read in the Greens function using a direct copy of the record indexing code from Fields
    for ISource = 1:NSourceDepths
        for IReceiver = 1:NReceiverDepths
            ISR = (ISource - 1) * NReceiverDepths + IReceiver;  
            IRecIn = 5 + ISR * NRecSr;
            for IRec = 1:NRecSr
                
                IStart = (IRec - 1) * NKCore + 1;
                IEnd = min([IStart + NKCore - 1, NKPts]);
                NRead =(IEnd - IStart + 1) * 2;   %Complex data so have to read in twice as many points
                lSeekRecord(InFile, IRecIn + IRec, RecordLength);
                [InDat, Count] = fread(InFile, [NRead, 1], Real);  
                if Count == NRead
                    ThisGKr = InDat(1:2:end) + i*InDat(2:2:end);
                    GKr(IStart:IEnd, IReceiver, ISource) = ThisGKr(1:NRead/2);
                else
                    uiwait(warndlg([int2str(Count), ' points read, ', int2str(NRead), ' points expected, ISource = ', ...
                            int2str(ISource), ', IReceiver = ' int2str(IReceiver)]));
                end
            end
        end
    end
    
    fclose(InFile);
end

function lSeekRecord(InFile, RecNumber, RecordLength)
%Seek to the start of the specified record
NEOR = 0;  %Number of bytes for end of record mark
BytesPerRecord = RecordLength * 4 + NEOR;

ByteOffset = BytesPerRecord * (RecNumber - 1);   %Offset from start of file

fseek(InFile, ByteOffset, 'bof');