plotslice.m

水声模型 很不错的东西

function [NewDir, Status, RxDepth] = PlotSlice(DefaultDir, FBase, FigNum, NRows, NCols, SubplNum, ...
   HoldCurrent, Style, RangeLimits, PltPhase, RxDepth)
% PlotSlice()
% plots a single TL slice from the shade file
% Modified by AJD 13/6/02 to read binary shade file 
%DefaultDir

Status = 1;
if nargin > 1
   filename = [DefaultDir FBase '.shd'];
   NewDir = DefaultDir;
   GotFileName = 1;
   GotRangeLimits = 1;
else
   GotRangeLimits = 0;
   if nargin < 1
       ChangeDir = 0;
   else
       StartDir = pwd;
       DirStat = 1;
       if ~exist(DefaultDir, 'dir')
           DirStat = LibMakeDirectory(DefaultDir);
       end
       if DirStat
           cd(DefaultDir);
           ChangeDir = 1;
       else
           ChangeDir = 0;
       end
    end

	[FName, Path] = uigetfile('*.shd', 'Select file to plot');

	if ChangeDir
   	    cd(StartDir);
	end

	if FName == 0
      NewDir = DefaultDir;
      GotFileName = 0;
	else
      GotFileName = 1;
   
   	filename = [Path FName];
      NewDir = Path;
   end
end

if GotFileName;
	isd = 1;  %Source number to plot

   [ pltitl, freq, nsd, nrd, nrr, sd, rd, rr, tlt, FileStatus ] = ReadShadeBin( filename );
   
   if FileStatus
      
	   Range = rr;

		if nrd > 1
            if (nargin < 11) | ((nargin >= 11) & isempty(RxDepth))
       	        Ans = inputdlg({'Enter desired depth (m) '}, ...
          	    ['Receiver depths are from ' num2str(rd(1)) ' to ' num2str(rd(end)) ' m'], 1, {'20'});
	            if isempty(Ans)
   	                return;
                end
      
   		        RxDepth = str2num(Ans{1});
            end
	   	    [MinDiff, ird] = min(abs(RxDepth - rd));
		else
   		    ird = 1;
            RxDepth = rd(1);
		end
   
   	if nargin == 1
		   Ans = inputdlg({...
   		      'Figure number', ...
      		   'Hold current plot (y/n)', ...
         	   'Line style', ...
	         'Plot phase as well as TL (y/n)'}, '', 1, {int2str(gcf), 'n', 'b-', 'n'});
		   if isempty(Ans)
      		return;
	   	end
   
   		FigNum = str2num(Ans{1});
   		HoldCurrent = strcmpi(Ans{2}, 'y');
	   	Style = Ans{3};
   	   PltPhase = strcmpi(Ans{4}, 'y');
      	NCols = 1;
	      NRows = 1;
   	   SubplNum = 1;
   	end
   
   	figure(FigNum);
	   if ~HoldCurrent      
            if exist('GUI_SetupGraphMenu.m', 'file')
                GUI_SetupGraphMenu('plot');
            end
	   end
   	if PltPhase
      	subplot(3,1,1);
	   else
   	   subplot(NRows, NCols, SubplNum);
	   end
   
   	if HoldCurrent      
      	hold on;
	   else
   	   hold off;
	   end
   
   
		plot( Range, -20*log10(abs(tlt( ird, : ))), Style)
	   if GotRangeLimits
   	   axis([RangeLimits(1) RangeLimits(2) -inf inf]);
	   end
		view( 0, -90 );
	   set(gca, 'Box', 'on');
   	xlabel( 'Range (m)' ); ylabel( 'Transmission Loss (dB)' );
	   title( [deblank( pltitl ) '. f=' num2str(freq) 'Hz, Zs=' num2str(sd(isd)) 'm, Zr=' num2str(rd(ird)) 'm.'] , ...
   	   'interpreter', 'none');
   
   	if PltPhase
      	K = 2*pi*freq/1500;
	      PrComplex = tlt(:, ird).';
   	   PrComplexRel = PrComplex .* exp(i*K*Range.');  %Removes factor of exp(ikr)
      	Phase = unwrap(angle( PrComplex));
	      PhaseRate = gradient(Phase) ./ gradient(Range.');
      
   	   PhaseR = unwrap(angle( PrComplexRel));
      	PhaseRRate = gradient(PhaseR) ./ gradient(Range.');
	      %PhCompare = unwrap(angle(exp(-i*K*rkm*1000)));
      
   	   subplot(3,1,2);
      	if HoldCurrent      
      		hold on;
	   	else
   	   	hold off;
	   	end

			plot( Range, Phase, Style);
      
   	   if ~HoldCurrent
	   	   hold on;
         	plot( Range, PhaseR,  'k:');
	         hold off;
   	   end
      
	   	if GotRangeLimits
   	   	axis([RangeLimits(1) RangeLimits(2) -inf inf]);
		   end
   	   set(gca, 'Box', 'on')
      	xlabel( 'Range (m)' ); 
	      ylabel( 'Phase (rad)' );
      
      
   	   subplot(3,1,3);
      	if HoldCurrent      
      		hold on;
	   	else
   	   	hold off;
	   	end

			plot( Range, PhaseRate, Style);
      
   	   if ~HoldCurrent
	   	   hold on;
         	plot( Range, PhaseRRate,  'k:');
	         hold off;
   	   end
      
	   	if GotRangeLimits
   	   	axis([RangeLimits(1) RangeLimits(2) -inf inf]);
		   end
   	   set(gca, 'Box', 'on');
      	xlabel( 'Range (m)' ); 
	      ylabel( 'Rate of phase change (rad/m)' );
      
   	end
   end
   
end