tp@0: function snapshot = EDB1makesnapshot(firstirfile,xvec,yvec,timestep,filtir,maxvalue,viewpoint,...
tp@0:     outlinecorners,noutlinepoints,sourcepoint,typeofir,windowsize,ampexp)
tp@0: % EDB1makesnapshot - Makes a snapshot of a sound field based on IRs.
tp@0: %
tp@0: % Input parameters:
tp@0: %   firstirfile     The first of the IR files
tp@0: %   xvec, yvec      The ranges of x- and y-values of the receiver positions
tp@0: %   timestep        The time step to plot
tp@0: %   filtir          A window to filter with
tp@0: %   maxvalue        [1,2] The minimum and maximum amplitude in the plot
tp@0: %   viewpoint       The point to view from [x y z]
tp@0: %   outlinecorners  A matrix of corners which will be plotted as an outline
tp@0: %                   in the form [x1 y1;x2 y2;x3 y3;...]
tp@0: %                   If the matrix has more columns than two, each pair of
tp@0: %                   column will be plotted as a separate outline, and they
tp@0: %                   will not be tied together.
tp@0: %   noutlinepoints  A list with the number of points in each column pair.
tp@0: %   sourcepoint     The [x y z] coordinates of the source.
tp@0: %   typeofir        't' (default), 'g' (geom), 'f' (direct sound), 'd'
tp@0: %                   (diffracted), 'c' (direct sound + geom)
tp@0: %   windowsize      's', 'm' or 'l' or 'x'
tp@0: %   ampexp          (optional) If a value is given here, the signal will be
tp@0: %                   plotted as abs(pressure)^ampexp, so that if ampexp =
tp@0: %                   0.5, sqrt(abs(pressure)) will be plotted. If no value
tp@0: %                   is given, or the value zero is given, the linear pressure will be plotted.
tp@0: %
tp@0: % Output parameters:
tp@0: %   snapshot           A snapshot, in the format which can be played by the
tp@0: %                   command surf(xvec,yvec,snapshot).
tp@0: %
tp@0: % Uses functions EDB1strpend
tp@0: %
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: %   This file is part of the Edge Diffraction Toolbox by Peter Svensson.                       
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify       
tp@0: %   it under the terms of the GNU General Public License as published by the Free Software     
tp@0: %   Foundation, either version 3 of the License, or (at your option) any later version.        
tp@0: %                                                                                              
tp@0: %   The Edge Diffraction Toolbox is distributed in the hope that it will be useful,       
tp@0: %   but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS  
tp@0: %   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.             
tp@0: %                                                                                              
tp@0: %   You should have received a copy of the GNU General Public License along with the           
tp@0: %   Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.                 
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20051013
tp@0: %
tp@0: % M = EDB1snapshot(firstirfile,xvec,yvec,starttimestep,ntimesteps,filtir,maxvalue,viewpoint
tp@0: %    outlinecorners,noutlinepoints,sourcepoint,typeofir,windowsize,ampexp);
tp@0: 
tp@0: disp('*************************************')
tp@0: disp('*')
tp@0: disp('*   Creating a snapshot....')
tp@0: disp(' ')
tp@0: 
tp@0: 
tp@0: %--------------------------------------------------------------
tp@0: % Extract the filename stem
tp@0: 
tp@0: if nargin == 0
tp@0: 	[firstirfile,irfilepath] = uigetfile('*ir.mat','Please select the first irfile');
tp@0:     [irfilepath,temp1,temp2] = fileparts(irfilepath);
tp@0: 	if ~isstr(firstirfile) | isempty(firstirfile)
tp@0: 		return
tp@0: 	end
tp@0: else
tp@0:     [irfilepath,firstirfile,fileext] = fileparts(firstirfile);
tp@0:     irfilepath = [irfilepath,filesep];
tp@0:     firstirfile = [firstirfile,fileext];
tp@0: end
tp@0: 
tp@0: if nargin < 11
tp@0:     typeofir = 't';    
tp@0:     windowsize = 's';
tp@0:     ampexp = 0;
tp@0: else
tp@0:     typeofir = lower(typeofir(1));    
tp@0:     if nargin < 12
tp@0:         windowsize = 's';        
tp@0:         ampexp = 0;
tp@0:     else
tp@0:         windowsize = lower(windowsize(1));                  
tp@0:         if nargin < 13
tp@0:             ampexp = 0;
tp@0:         else
tp@0:             if ampexp < 0
tp@0:                 error('ERROR: ampexp must be larger than, or equal to, zero')    
tp@0:             end
tp@0:         end
tp@0:     end
tp@0: end
tp@0: 
tp@0: filestem = EDB1strpend(firstirfile,'_ir');
tp@0: iv = find(filestem=='_');
tp@0: iv = iv(length(iv));
tp@0: firstnumber = str2num(filestem(iv+1:length(filestem)));
tp@0: filestem = filestem(1:iv);
tp@0: 
tp@0: nx = length(xvec);
tp@0: ny = length(yvec);
tp@0: nfiles = nx*ny;
tp@0: 
tp@0: loadcommandbase = ['load ',irfilepath,filestem];
tp@0:     
tp@0: %--------------------------------------------------------------
tp@0: % Read the files, filter and store
tp@0: 
tp@0: Bigir = zeros(nfiles,1);
tp@0: 
tp@0: disp('   Loading the files, starting with:')
tp@0: for ii = 1:nfiles
tp@0:     eval([loadcommandbase,sprintf('%.0f',ii),'_ir.mat'])
tp@0:     if typeofir == 'f'
tp@0:         irtot = irdirect;    
tp@0:     elseif typeofir == 'g'
tp@0:         irtot = irgeom;
tp@0:     elseif typeofir == 'c'
tp@0:         irtot = irgeom + irdirect;
tp@0:     elseif typeofir == 'd',     
tp@0:         irtot = irdiff;
tp@0:     end
tp@0:         
tp@0:     if any(irtot)
tp@0:         
tp@0:         if length(irtot) > timestep
tp@0:             irtot = irtot(1:timestep);    
tp@0:         elseif length(irtot) < timestep
tp@0:             irtot = [irtot;zeros(timestep-length(irtot),1)];            
tp@0:         end
tp@0: 
tp@0:         irsamplevalue = sum(irtot([timestep:-1:timestep-length(filtir)+1]).*filtir);
tp@0:     else
tp@0:         irsamplevalue = 0;
tp@0:     end
tp@0:     
tp@0:     if ampexp == 0
tp@0:         Bigir(ii) = irsamplevalue;
tp@0:     else
tp@0:         Bigir(ii) = (abs(irsamplevalue)).^ampexp;
tp@0:     end
tp@0:     
tp@0: end
tp@0: 
tp@0: disp('   ... files are loaded.')
tp@0: 
tp@0: 
tp@0: if length(maxvalue) < 2
tp@0:     maxvalue = [maxvalue(1) max(Bigir)];    
tp@0: end
tp@0: 
tp@0: %--------------------------------------------------------------
tp@0: % Prepare the outline plotting
tp@0: 
tp@0: opengl neverselect
tp@0: 
tp@0: snapshot = reshape(Bigir,ny,nx);
tp@0: 
tp@0: [noutlinerows,ncolumns] = size(outlinecorners);
tp@0: noutlines = floor(ncolumns/2);
tp@0: if noutlines == 1
tp@0:     outlinecorners = [outlinecorners;outlinecorners(1,:)];
tp@0:     noutlinepoints = noutlinepoints+1;
tp@0: end
tp@0: 
tp@0: axisvalues = [min(xvec) max(xvec) min(yvec) max(yvec) maxvalue(1) maxvalue(2)];
tp@0: figure(1)
tp@0: clf
tp@0: 
tp@0: xyaspectratio = abs(max(xvec)-min(xvec))/abs(max(yvec)-min(yvec));
tp@0: 
tp@0: windowheight = 675;
tp@0: windowwidth = windowheight*xyaspectratio;
tp@0: 
tp@0: if windowwidth > 800
tp@0:     scaledownfactor = windowwidth/800;
tp@0:     windowwidth = windowwidth/scaledownfactor;
tp@0:     windowheight = windowheight/scaledownfactor;
tp@0: end
tp@0: 
tp@0: windowpos = [380 80];
tp@0: if windowsize == 's'
tp@0:     windowwidth = windowwidth/2;
tp@0:     windowheight = windowheight/2;
tp@0: elseif windowsize == 'm'
tp@0:     windowwidth = windowwidth/1.41;
tp@0:     windowheight = windowheight/1.41;
tp@0: elseif windowsize == 'x'
tp@0:     windowwidth = windowwidth*1.41;
tp@0:     windowheight = windowheight*1.41;
tp@0:     windowpos(1) = 100;
tp@0: end
tp@0: 
tp@0: set(1,'Position',[windowpos(1:2)   windowwidth   windowheight])
tp@0: 
tp@0: snapshot = reshape(Bigir,ny,nx);
tp@0: surf(xvec,yvec,snapshot);
tp@0: 
tp@0: shading interp
tp@0: colormap('jet')    
tp@0: caxis([maxvalue(1) maxvalue(2)])
tp@0: axis(axisvalues);
tp@0: axis off
tp@0: view(viewpoint)
tp@0: 
tp@0: if noutlines > 0
tp@0: 	for kk = 1:noutlines
tp@0: 		for ll = 1:noutlinepoints(kk)-1
tp@0: 			H = line([outlinecorners(ll,(kk-1)*2+1) outlinecorners(ll+1,(kk-1)*2+1)],[outlinecorners(ll,(kk-1)*2+2) outlinecorners(ll+1,(kk-1)*2+2)],[0 0]);
tp@0: 			set(H,'LineWidth',4);
tp@0: 			set(H,'Color',[1 1 0]);
tp@0: 			H = line([outlinecorners(ll,(kk-1)*2+1) outlinecorners(ll+1,(kk-1)*2+1)],[outlinecorners(ll,(kk-1)*2+2) outlinecorners(ll+1,(kk-1)*2+2)],2*[maxvalue(2) maxvalue(2)]);
tp@0: 			set(H,'LineWidth',4);
tp@0: 			set(H,'Color',[1 1 0]);
tp@0: 			H = line([outlinecorners(ll,(kk-1)*2+1) outlinecorners(ll,(kk-1)*2+1)],[outlinecorners(ll,(kk-1)*2+2) outlinecorners(ll,(kk-1)*2+2)],[0 2*maxvalue(2)]);        
tp@0: 			set(H,'LineWidth',4);
tp@0: 			set(H,'Color',[1 1 0]);
tp@0: 		end
tp@0: 		H = line([outlinecorners(ll+1,(kk-1)*2+1) outlinecorners(ll+1,(kk-1)*2+1)],[outlinecorners(ll+1,(kk-1)*2+2) outlinecorners(ll+1,(kk-1)*2+2)],[0 2*maxvalue(2)]);        
tp@0: 			set(H,'LineWidth',4);
tp@0: 			set(H,'Color',[1 1 0]);
tp@0: 	end
tp@0: end
tp@0: if ~isempty(sourcepoint)
tp@0: 	H = line([sourcepoint(1)*0.95 sourcepoint(1)*1.05],[sourcepoint(2)*0.95 sourcepoint(2)*1.05],[maxvalue(2) maxvalue(2)]);
tp@0: 			set(H,'LineWidth',4);
tp@0: 			set(H,'Color',[1 0 0]);        
tp@0: 	H = line([sourcepoint(1)*0.95 sourcepoint(1)*1.05],[sourcepoint(2)*1.05 sourcepoint(2)*0.95],[maxvalue(2) maxvalue(2)]);
tp@0: 			set(H,'LineWidth',4);
tp@0: 			set(H,'Color',[1 0 0]);        
tp@0: end