Mercurial > hg > human-echolocation
view private/EDB1makeirs.m @ 18:2d5f50205527 jabuilder_int tip
Escape the trailing backslash as well
| author | Chris Cannam |
|---|---|
| date | Tue, 30 Sep 2014 16:23:00 +0100 |
| parents | 90220f7249fc |
| children |
line wrap: on
line source
function edirfile = EDB1makeirs(edpathsfile,specorder,Rstart,... EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,... edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,reftoshortlistE,re1sho,re2sho,... thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,includedirectsound,... saveindividualdiffirs) % EDB1makeirs - Constructs impulse responses from a list of paths in a file. % % Input parameters: % edpathsfile The name of the file containing the plane and edge data % specorder The maximum order of reflections that is calculated. % Rstart The reference distance that the time zero of the impulse % response refers to. % EDcalcmethod 'n' or 'v' or 'k', meaning: % 'n': the new method by Svensson et al % 'v': Vanderkooys method % 'k': the Kirchhoff diffraction approximation % edgestartcoords A matrix [nedges,3] of the startpoint coordinates % of % the edges. % edgeendcoords A matrix [nedges,3] of the endpoint coordinates of % the edges. % edgenvecs A matrix [nedges,3] of the normal vectors of % the reference plane of each edge. % edgelengthvec A list [nedge,1] of the lenghts of each edge. % planeeqs A matrix, [nplanes,4], of all the plane equations. % approxplanemidpoints A matrix, [nplanes,3], of all plane midpoint % coordinates. It may be an approximate location of % the midpoint. % reflfactors A list [nplanes,1] of the reflection factors of % each plane. % closwedangvec A list [nedge,1] of the close-wedge angle of each % edge. % planesatedge A matrix [nedge,2] of the two planes of each edge % with the reference plane first. % elemsize A list [1,difforder] with the relative edge element % sizes that will be used for each order. The first % value, for diffraction order 1, is redundant and % not used. For higher orders, a value of, e.g., 0.5 % means that edges will be subdivided into elements % such that the wavelength at half the sampling % frequency is 0.5 times the element size. A lower % value gives faster but less accurate results. % reftoshortlistE A matrix, [nedges,nedges] with a pointer to the % short lists that contain edge-to-edge data. % re1sho A short list, [nshort,1] of the cylindrical radii % of each edge startpoint relative to all other % edges. % re2sho A short list, [nshort,1] of the cylindrical radii % of each edge endpoint relative to all other % edges. % thetae1sho A short list, [nshort,1] of the theta angle % of each edge startpoint relative to all other % edges. % thetae2sho A short list, [nshort,1] of the theta angle % of each edge endpoint relative to all other % edges. % ze1sho A short list, [nshort,1] of the z-value % of each edge startpoint relative to all other % edges. % ze2sho A short list, [nshort,1] of the z-value % of each edge endpoint relative to all other % edges. % edgeseespartialedge A matrix, [nedges,nedges], with edge-to-edge % visibility values. The values are 0 (invisible) % up to (2^(nedgesubs)-1)^2 (full visibility). % A pos. value indicates that the edge-to-edge path % runs along a plane; a neg. avlue indicates that it % does not run along a plane. % edgeplaneperptoplane1 A matrix, [nplanes,nedges], with 0 or 1. A % value 1 indicates that one of the edge's two % defining planes is perpendicular to another plane. % desiredirfile The file name that will be given to the % output file. % guiderowstouse (optional) A list of values, indicating which rows in the % mainlistguide and mainlistguidepattern that should % be used. This way it is possible to select only % diffraction, for instance. If this list is not % specified, all components will be used. % includedirectsound (optional) 0 or 1, indicating whether the direct % sound should be included or not. Default: 1 % saveindividualdiffirs (optional) Vector of length 2 with values 0 or 1 % Pos 1 indicating whether % 0 - all orders of diffraction IR:s will be summed % in a single vector 'irdiff', or % 1 - each order of diffraction IR:s will be placed % in individual columns in a matrix 'irdiff'. % Pos 2 indicating whether % 0 - only the total diffraction ir will be saved in the file. % 1 - all individual diffraction irs will be saved in a large % matrix alldiffirs. % Default: [0 0]. % FSAMP, CAIR, SHOWTEXT Global parameters. % % Output parameters: % edirfile The filename used for the output file that contains % the impulse responses. % % Data in the output file: % irdirect The IR containing the direct sound, size [nirlength,1]. % irgeom The IR containing the specular reflections, size [nirlength,1]. % irdiff The IR containing the diffracted components, size [nirlength,1]. % irtot The IR containing the sum of all components, size [nirlength,1]. % alldiffirs (optional) Matrix with all individual first-order diffraction IRs % on individual lines. % alldiffirsextradata (optional) Vector with combination number (form reflpaths) that % matches alldiffirs. % FSAMP Rstart CAIR Same as the input parameters. % % Uses functions EDB1irfromslotlist EDB1calcdist EDB1coordtrans1 EDB1coordtrans2 % EDB1wedge1st_int EDB1wedge2nd % % ---------------------------------------------------------------------------------------------- % This file is part of the Edge Diffraction Toolbox by Peter Svensson. % % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by the Free Software % Foundation, either version 3 of the License, or (at your option) any later version. % % The Edge Diffraction Toolbox is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. % % You should have received a copy of the GNU General Public License along with the % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>. % ---------------------------------------------------------------------------------------------- % Peter Svensson (svensson@iet.ntnu.no) 20061208 % % edirfile = EDB1makeirs(edpathsfile,specorder,... % Rstart,EDcalcmethod,edgestartcoords,edgeendcoords,edgenvecs,... % edgelengthvec,planeeqs,approxplanemidpoints,reflfactors,closwedangvec,planesatedge,elemsize,reftoshortlistE,re1sho,re2sho,... % thetae1sho,thetae2sho,ze1sho,ze2sho,edgeseespartialedge,edgeplaneperptoplane1,desiredirfile,guiderowstouse,includedirectsound,... % saveindividualdiffirs); global SHOWTEXT FSAMP CAIR BIGEDGESTEPMATRIX global IRDIFFVEC eval(['load ',edpathsfile]) %-------------------------------------------------------------------------- % We look for the optional vector multfactors in the edpathsfile. % If it was there, then those values will be used to switch off, or boost % any reflection path. % Those values will only be used for diffraction paths. if exist('multfactors') == 0 multfactors = ones(size(reflpaths,1),1); else if size(multfactors,1) ~= size(reflpaths,1) error(['The edpaths file contained a vector multfactors which does not have the same length as reflpaths.']) end end %-------------------------------------------------------------------------- edirfile = desiredirfile; Sdirection = [1 0 0]; maxnedgesorplanes = max(max(reflpaths)); if ~isempty(maxnedgesorplanes) multfac = 10^ceil(log10(double(maxnedgesorplanes)+1)); else multfac = 0; end if nargin < 27 saveindividualdiffirs = [0 0]; end if nargin < 26 includedirectsound = 1; end if nargin < 25 guiderowstouse = []; end if isempty(guiderowstouse) usesubset = 0; if SHOWTEXT > 2 disp('Using all components') end else usesubset = 1; if SHOWTEXT > 2 disp('Using only some components') end end nyvec = pi./(2*pi - closwedangvec); onesvec1 = ones(1,3); souspecboost = 1; if ~isempty(Sinsideplanenumber) if reflfactors(Sinsideplanenumber(1)) ~= 0 souspecboost = 2; end end recspecboost = 1; if ~isempty(Rinsideplanenumber) if reflfactors(Rinsideplanenumber(1)) ~= 0 recspecboost = 2; end end if isempty(re1sho) userwantsdiff2 = 0; else userwantsdiff2 = 1; end %------------------------------------------------------- if exist('mainlistguide') ~= 1 mainlistguide = []; else mainlistguide = double(mainlistguide); end [ncomponents,ncols] = size(reflpaths); [nrowsinlist,slask] = size(mainlistguide); if ncols > 1 difforderinlist = sum(mainlistguidepattern.'=='d').'; else difforderinlist = (mainlistguidepattern=='d'); end lastNdiffrow = zeros(1,specorder); for ii = 1:specorder iv = find(difforderinlist <= ii ); if ~isempty(iv) lastNdiffrow(ii) = iv(end); end end nrefltogetherwdiff = sum(mainlistguidepattern.'=='d').'.*( sum(mainlistguidepattern.'=='s').'); if SHOWTEXT >= 3 disp([' Constructing IR. ',int2str(ncomponents),' components.']) end irdirect = [0 0].'; irgeom = [0 0].'; irdiff = [0 0].'; irtot = [0 0].'; Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR EDcalcmethod']; if ncomponents == 0 eval(['save ',edirfile,Varlist]) return end %############################################################################## %############################################################################## %############################################################################## %############################################################################## % % Diffraction once, possibly with pre- and post-specular reflections % %############################################################################## directsoundonboundary = 0; specreflonboundary = zeros(size(planeeqs,1),1); if firstdiffrow ~= 0 % First we remove the 'd' for all the rows in the mainlistguidepattern that % we do not want to use. if usesubset == 1 singdiffrows = find(difforderinlist==1); rowstoclear = setdiff(singdiffrows,guiderowstouse); if ~isempty(rowstoclear) mainlistguidepattern(rowstoclear,:) = mainlistguidepattern(rowstoclear,:)*0; end end % Here we can use the single-diffraction calculation for all % combinations. The source should either be the original source or % an image source. The receiver should either be the original % receiver or an image receiver. % First we create lists that specify whether there are specular % reflections before the edge diffraction, and after the edge % diffraction. Also, a list which gives the column where the edge % diffraction is present can be used to efficiently extract the % edge number. % NB! The lists prespeclist, postspeclist and singdiffcol all have % the (short) length of the mainlistguide. % The data in these lists will then be used to create the long % lists prespecpresent and postspecpresent which simply are 0 or 1. % A matrix called edgemask will have the same number of columns as % the reflpaths matrix. For each row, there will be a 1 in the % column where the edge diffraction is so we can efficiently find % the edge number. multmatrix = [1:specorder]; if size(mainlistguidepattern,2) > specorder multmatrix = [1:size(mainlistguidepattern,2)]; end multmatrix = multmatrix( ones(nrowsinlist,1),: ); if ncols > 1 singdiffcol = sum( (multmatrix.*(mainlistguidepattern=='d')).' ).'; singdiffcol = singdiffcol.*(sum(mainlistguidepattern.'=='d').'<=1); nrefl = sum( (mainlistguidepattern=='s' | mainlistguidepattern=='d').' ).'; else singdiffcol = mainlistguidepattern=='d'; nrefl = singdiffcol | mainlistguidepattern=='s'; end prespeclist = singdiffcol-1; prespeclist = prespeclist.*(prespeclist>0); postspeclist = (nrefl-singdiffcol).*(singdiffcol~=0); prespecpresent = zeros(ncomponents,1); postspecpresent = zeros(ncomponents,1); edgemask = zeros(ncomponents,specorder); diffrowsthatareOK = [firstdiffrow:lastNdiffrow(1)]; if usesubset == 1 diffrowsthatareOK = intersect(diffrowsthatareOK,guiderowstouse); end for ii = diffrowsthatareOK, %firstdiffrow:lastNdiffrow(1) iv = [(mainlistguide(ii,2)):(mainlistguide(ii,3))]; onesvec2 = ones(length(iv),1); edgemask(iv,singdiffcol(ii)) = onesvec2; prespecpresent(iv) = (prespeclist(ii)>0)*onesvec2; postspecpresent(iv) = (postspeclist(ii)>0)*onesvec2; end prespecpresent = prespecpresent(:,onesvec1); postspecpresent = postspecpresent(:,onesvec1); if ncols > 1 edgenumb = sum( (edgemask.*double(reflpaths(:,1:specorder))).').'; else edgenumb = edgemask.*double(reflpaths); end nedgeirs = size(edgenumb,1); nnonzeroedgeirs = sum(edgenumb(:,1)>0); ircounter = 0; edgelist = unique(edgenumb); for ii = 1: length(edgelist) edge = edgelist(ii); if edge~= 0 iv = find(edgenumb==edge); ncombs = length(iv); if ncombs > 0 & any(multfactors(iv)) if SHOWTEXT >= 4 disp([' Edge ',int2str(edge),': ',int2str(ncombs),' combinations']) end onesvec3 = ones(ncombs,1); IS = full(specextradata(iv,1:3)); IR = full(specextradata(iv,4:6)); edgestart = full(edgeextradata(iv,1)); edgeend = full(edgeextradata(iv,2)); % Calculate the cyl coordinates of all IS % % Calculate the cyl coord of all IR edgecoords = [edgestartcoords(edge,:);edgeendcoords(edge,:)]; [rs,thetas,zs,rr,thetar,zr] = EDB1coordtrans2(IS,IR,edgecoords,edgenvecs(edge,:)); bc = reflfactors(planesatedge(edge,:)).'; if prod(double(bc==[1 1])) ~= 1 disp(' Non-rigid wedge') end for jj = 1:ncombs if multfactors(iv(jj)) > 0 if EDcalcmethod(1) == 'n' [irnew,slask,singularterm] = EDB1wedge1st_int(FSAMP,closwedangvec(edge),rs(jj),thetas(jj),zs(jj),rr(jj),thetar(jj),zr(jj),... edgelengthvec(edge)*[edgestart(jj) edgeend(jj)],EDcalcmethod,Rstart,bc); if any(singularterm) if singularterm(2) | singularterm(3) directsoundonboundary = 1; elseif singularterm(1) if specorder == 1 specreflonboundary(planesatedge(edge,2)) = 1; else disp('WARNING! A specular refl. of order > 1 is half obscured but this is not handled yet!'); end elseif singularterm(4) if specorder == 1 specreflonboundary(planesatedge(edge,1)) = 1; else disp('WARNING! A specular refl. of order > 1 is half obscured but this is not handled yet!'); end end end else % ... if EDcalcmethod(1) == 'n' [irnew,slask,singularterm] = EDB1wedge1stcombo(FSAMP,closwedangvec(edge),rs(jj),thetas(jj),zs(jj),rr(jj),thetar(jj),zr(jj),... edgelengthvec(edge)*[edgestart(jj) edgeend(jj)],EDcalcmethod,Rstart,bc); end % ... if EDcalcmethod(1) == 'n' % Decide whether the IR should be boosted or not % % The factors souspecboost and recspecboost have % values other than 1 if the source or receiver is directly % at a plane, for the boosting of the direct sound, and % specular reflections. % % This boost factor should be used also for edge % diffraction if: % 1. There are some specular reflections between the % source and the edge % or % 2. There are no specular reflections between the % source and the edge, but the source/receiver is at % a plane which does not belong to the edge. if souspecboost ~= 1 | recspecboost ~= 1 boostfactor = 1; if prespecpresent(iv(jj)) == 1 boostfactor = souspecboost; else if ~isempty(Sinsideplanenumber) if Sinsideplanenumber(1)~=planesatedge(edge,1) & Sinsideplanenumber(1)~=planesatedge(edge,2), boostfactor = souspecboost; end end end if postspecpresent(iv(jj)) == 1 boostfactor = boostfactor*recspecboost; else if ~isempty(Rinsideplanenumber) if Rinsideplanenumber(1)~=planesatedge(edge,1) & Rinsideplanenumber(1)~=planesatedge(edge,2), boostfactor = boostfactor*recspecboost; end end end if boostfactor ~= 1 irnew = irnew*boostfactor; end end % .... if souspecboost ~= 1 | recspecboost ~= 1 ndiff = length(irdiff); nnew = length(irnew); ircounter = ircounter + 1; if nnew > ndiff irdiff = [irdiff;zeros(nnew-ndiff,1)]; end irdiff(1:nnew) = irdiff(1:nnew) + irnew*double(multfactors(iv(jj))); if saveindividualdiffirs(2) == 1 if exist('alldiffirs','var') == 0 alldiffirs = sparse(zeros(nnonzeroedgeirs,nnew)); alldiffirs(1,:) = irnew*double(multfactors(iv(jj))).'; alldiffirsextradata = zeros(nnonzeroedgeirs,1); alldiffirsextradata(1) = iv(jj); else if nnew > ndiff alldiffirs = [alldiffirs zeros(nnonzeroedgeirs,nnew-ndiff)]; end alldiffirs(ircounter,1:nnew) = irnew*double(multfactors(iv(jj))).'; alldiffirsextradata(ircounter) = iv(jj); end end if SHOWTEXT >= 7 figure(1) plot(irnew) figure(2) plot(irdiff) save ~/Documents/Temp/irdiff2.mat pause end end % ... if multfactors(iv) > 0 end % ..... for jj = 1:ncombs end %.... if ncombs > 0 & any(multfactors(iv)) end % .... if edge~= 0 end % .... for ii = 1: length(edgelist) end % ...if firstdiffrow ~= 0 nspecreflonboundary = sum(specreflonboundary); %############################################################################## %############################################################################## %############################################################################## %############################################################################## % % The direct sound % %############################################################################## if usesubset == 1 & directsoundrow == 1 if any(guiderowstouse == 1) == 0 directsoundrow = 0; end end if directsoundrow == 1 & includedirectsound == 1 dist = norm(R-S); slotnumberfrac = (dist-Rstart)/CAIR*FSAMP+1; amp = souspecboost*recspecboost./dist; if directsoundonboundary amp = amp/2; if SHOWTEXT >= 4 disp('HALVING DIRECT SOUND') end end slotnumber = floor(slotnumberfrac); if any(slotnumber<1) error('ERROR: Rstart is set to too low a value!') end slotnumberfrac = slotnumberfrac - slotnumber; irdirect = zeros( slotnumber+2,1 ); irdirect(slotnumber) = amp.*(1-slotnumberfrac) ; irdirect(slotnumber+1) = amp.*slotnumberfrac; if ncomponents == 1 nnew = length(irdirect); if nnew > 2 irdiff = [irdiff;zeros(nnew-2,1)]; irtot = [irtot;zeros(nnew-2,1)]; irgeom = [irgeom;zeros(nnew-2,1)]; end irtot(1:nnew) = irtot(1:nnew) + irdirect; irtot = sparse(irtot); irgeom = sparse(irgeom); irdirect = sparse(irdirect); irdiff = sparse(irdiff); eval(['save ',edirfile,Varlist]) return end else if directsoundonboundary == 1 dist = norm(R-S); slotnumberfrac = (dist-Rstart)/CAIR*FSAMP+1; amp = souspecboost*recspecboost./dist; if directsoundonboundary amp = amp/2; if SHOWTEXT >= 4 disp('HALVING DIRECT SOUND') end end slotnumber = floor(slotnumberfrac); if any(slotnumber<1) error('ERROR: Rstart is set to too low a value!') end slotnumberfrac = slotnumberfrac - slotnumber; irdirect = zeros( slotnumber+2,1 ); irdirect(slotnumber) = amp.*(1-slotnumberfrac) ; irdirect(slotnumber+1) = amp.*slotnumberfrac; end if SHOWTEXT >= 4 disp([' No direct sound']) end end %############################################################################## %############################################################################## %############################################################################## %############################################################################## % % All-specular s, ss, sss, etc % %############################################################################## if allspecrows(1) ~= 0 if usesubset == 1 specrowsthatareOK = intersect(allspecrows,guiderowstouse); else specrowsthatareOK = [allspecrows(1):allspecrows(2)]; end if ~isempty(specrowsthatareOK) if usesubset == 1 temp = mainlistguide(specrowsthatareOK,2:3); ivspec = [double(mainlistguide(specrowsthatareOK(1),2)):double(mainlistguide(specrowsthatareOK(1),3))]; for ii = 2:size(temp,1); ivspec = [ivspec [double(mainlistguide(specrowsthatareOK(ii),2)):double(mainlistguide(specrowsthatareOK(ii),3))]]; end else ivspec = [mainlistguide(allspecrows(1),2):mainlistguide(allspecrows(2),3)]; end if nspecreflonboundary > 0 listofspecreflonboundary = find(specreflonboundary); specrefltoadd = setdiff(listofspecreflonboundary,ivspec) if ~isempty(specrefltoadd) error(['ERROR: We need to add some code here, to reinsert pruned spec. refl.']) end end if SHOWTEXT >= 4 disp([' ',int2str(length(ivspec)),' specular reflections']) end dist = EDB1calcdist(full(specextradata(ivspec,1:3)),R); if nspecreflonboundary > 0 specamp = ones(size(souspecboost)); amplitudeshouldbehalf = ismember(reflpaths(ivspec,1),listofspecreflonboundary); specamp = specamp - amplitudeshouldbehalf*0.5; irnew = specamp*souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist); else irnew = souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist); end ngeom = length(irgeom); nnew = length(irnew); if nnew > ngeom irgeom = [irgeom;zeros(nnew-ngeom,1)]; end irgeom(1:nnew) = irgeom(1:nnew) + irnew; else if nspecreflonboundary > 0 listofspecreflonboundary = find(specreflonboundary); error(['ERROR: We need to add some code here, to reinsert pruned spec. refl., pos. 2']) end end else if nspecreflonboundary > 0 listofspecreflonboundary = find(specreflonboundary); [xis] = EDB1findis(S,listofspecreflonboundary,planeeqs,1,[1 1 1]); disp(['WARNING: We have some new code here, to reinsert pruned spec. refl., pos. 3']) dist = EDB1calcdist(xis,R); irnew = 0.5*souspecboost*recspecboost*EDB1irfromslotlist((dist-Rstart)/CAIR*FSAMP+1,1./dist); ngeom = length(irgeom); nnew = length(irnew); if nnew > ngeom irgeom = [irgeom;zeros(nnew-ngeom,1)]; end irgeom(1:nnew) = irgeom(1:nnew) + irnew; end end %############################################################################## %############################################################################## %############################################################################## %############################################################################## % % Multiple diffraction, possibly with pre- and post-specular reflections % %############################################################################## if userwantsdiff2 == 1 JJ = setstr(32*ones(specorder,1)); for jj=1:specorder jjstr = int2str(jj); JJ(jj,1:length(jjstr)) = jjstr; end [n1,n2] = size(JJ); for Ndifforder = 2:specorder if SHOWTEXT >= 3 disp([' Diffraction order ',int2str(Ndifforder)]) end if any(difforderinlist==Ndifforder) & elemsize(Ndifforder) > 0 % Calculate some general parameters that are shared for all % N-diffraction calculations divmin = CAIR/(FSAMP*elemsize(Ndifforder)); ndivvec = ceil(abs( edgelengthvec.' )/divmin); dzvec = (edgelengthvec.')./ndivvec; ncylrows = 4*(Ndifforder-1); % Here we can use the double-diffraction calculation for all % combinations. The source should either be the original source or % an image source. The receiver should either be the original % receiver or an image receiver. noffset = lastNdiffrow(Ndifforder-1); ivNdiff = [noffset+1:lastNdiffrow(Ndifforder)]; ndiffcombs = length(ivNdiff); zerosvec1 = zeros(ndiffcombs,1); ndoubrows = length(ivNdiff); previousrow = lastNdiffrow(Ndifforder-1); if previousrow > 0 noffsetlonglist = mainlistguide(previousrow,3); else noffsetlonglist = 0; end nremaining = ncomponents - (mainlistguide(lastNdiffrow(Ndifforder),3)); nlonglist = ncomponents-noffsetlonglist-nremaining; ivlonglist = [mainlistguide(ivNdiff(1),2):mainlistguide(ivNdiff(end),3)].'; diffcols = zerosvec1(:,ones(1,Ndifforder)); for ii = ivNdiff(1):ivNdiff(end) diffcols(ii-ivNdiff(1)+1,:) = find(mainlistguidepattern(ii,:)=='d'); end nreflorder = sum( (mainlistguidepattern(ivNdiff(1):ivNdiff(end),:)=='d'|mainlistguidepattern(ivNdiff(1):ivNdiff(end),:)=='s').').'; nprespecs = diffcols(:,1)-1; nmidspecs = diffcols(:,2)-diffcols(:,1)-1; npostspecs = nreflorder-diffcols(:,Ndifforder); if ndiffcombs > 1 ivreftolonglist = ivlonglist(:,ones(1,ndiffcombs)); comppattern = mainlistguide(ivNdiff,2).'; comppattern = comppattern(ones(nlonglist,1),:); rowinpatternlist = sum( (ivreftolonglist>=comppattern).' ).'; else rowinpatternlist = ones(size(ivlonglist)); end % Construct a long matrix with the edge numbers extracted from the % reflpaths matrix. longnmidspec = zeros(nlonglist,1); iv = find(nmidspecs>0); for ii = 1:length(iv) ivreftolonglist = [mainlistguide(ivNdiff(iv(ii)),2):mainlistguide(ivNdiff(iv(ii)),3)] - noffsetlonglist; longnmidspec(ivreftolonglist) = nmidspecs(iv(ii)); end %------------------------------------------------------------------ % edgepattern will be a matrix which, for each reflection combination, contains % the 2,3,4,... edge numbers that are involved in each path % regardless of there are specular reflections before, in the % middle, or after. edgepattern = zeros(nlonglist,Ndifforder); countvec = [1:nlonglist].'; reflpathscut = reflpaths(ivlonglist,:); for ii = 1:Ndifforder ivrefcol = countvec + (diffcols(rowinpatternlist,ii)-1)*nlonglist; edgepattern(:,ii) = reflpathscut(ivrefcol); end %------------------------------------------------------------------ % midspecpattern will be a matrix which, for each reflection combination, contains % the 1,2,3,4,... specular reflection numbers (i.e., plane numbers) that are involved % in each path between diffraction 1 and diffraction 2. % First, the reflection component immediately after diffraction 1 % is selected for every reflection. Then there is a masking which % zeroes the midspecpattern value for all the combinations that % don't have any midspec. % NB! For combinations like ssssdd, the first step will point to a % column which is outside the matrix. This is fixed by maximizing % the column number. midspecpattern = zeros(nlonglist,max([specorder-Ndifforder 1])); maxpointvalue = nlonglist*max([specorder-Ndifforder 1]); for ii = 1:specorder-Ndifforder ivrefcol = countvec + (diffcols(rowinpatternlist,1)+(ii-1))*nlonglist; ivrefcol = mod(ivrefcol-1,maxpointvalue)+1; midspecpattern(:,ii) = double(reflpathscut(ivrefcol)).*(longnmidspec>=ii); end diff1col = diffcols(rowinpatternlist,1); % Many combinations include the same edge-pair/N-let, so we extract the % unique edge pairs/N-lets, and go through them in the ii-loop [edgeshortlist,ivec,jvec] = unique(edgepattern,'rows'); lastndivcomb = zeros(1,Ndifforder); for ii = 1: length(ivec) if SHOWTEXT >= 3 if round(ii/ceil(length(ivec)/10))*ceil(length(ivec)/10) == ii disp([' Combination no. ',int2str(ii),' of ',int2str(length(ivec))]) end end iv = find(jvec==ii); ncombs = length(iv); onesvec3 = ones(ncombs,1); if SHOWTEXT >= 4 numvec = int2str(edgeshortlist(ii,1)); for ll = 2:Ndifforder numvec = [numvec,' ', int2str(edgeshortlist(ii,ll))]; end disp([' Edges ',numvec]) end if Ndifforder >= 2 & any(multfactors(ivlonglist(iv))) newndivvec = ndivvec(edgeshortlist(ii,:)); if any(lastndivcomb~=newndivvec) if Ndifforder == 2 ivmatrix = EDB1creindexmatrix(newndivvec); else ivmatrix = EDB1creindexmatrix(newndivvec(2:end)); end [nedgeelcombs,slask] = size(ivmatrix); if Ndifforder == 2 BIGEDGESTEPMATRIX = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),:); else BIGEDGESTEPMATRIX = (double(ivmatrix)-0.5)./newndivvec(uint8(ones(nedgeelcombs,1)),2:end); end clear ivmatrix lastndivcomb = newndivvec; end end pathalongplane = (edgeseespartialedge(edgeshortlist(ii,2),edgeshortlist(ii,1))>0); for jj = 3:Ndifforder pathalongplane = [pathalongplane,(edgeseespartialedge(edgeshortlist(ii,jj),edgeshortlist(ii,jj-1))>0)]; end if any(multfactors(ivlonglist(iv))) IS = full(specextradata(ivlonglist(iv),1:3)); IR = full(specextradata(ivlonglist(iv),4:6)); firstedgestart = full(edgeextradata(ivlonglist(iv),1)); firstedgeend = full(edgeextradata(ivlonglist(iv),2)); lastedgestart = full(edgeextradata(ivlonglist(iv),3)); lastedgeend = full(edgeextradata(ivlonglist(iv),4)); % Calculate the cyl coordinates of all IS/S and IR/R cylS = zeros(ncombs,3); edgecoords = [edgestartcoords(edgeshortlist(ii,1),:);edgeendcoords(edgeshortlist(ii,1),:)]; [cylS(:,1),cylS(:,2),cylS(:,3)] = EDB1coordtrans1(IS,edgecoords,edgenvecs(edgeshortlist(ii,1),:)); cylR = zeros(ncombs,3); edgecoords = [edgestartcoords(edgeshortlist(ii,Ndifforder),:);edgeendcoords(edgeshortlist(ii,Ndifforder),:)]; [cylR(:,1),cylR(:,2),cylR(:,3)] = EDB1coordtrans1(IR,edgecoords,edgenvecs(edgeshortlist(ii,Ndifforder),:)); bc = ones(Ndifforder,2); % Check real reflfactors!!! bc = reshape(bc.',2*Ndifforder,1); % Pick out the edge-to-edge coordinates for this specific % edge pair/N-let. if ~isempty(reftoshortlistE), if Ndifforder >= 2 index1 = reftoshortlistE(edgeshortlist(ii,2),edgeshortlist(ii,1)); cylE2_r1 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)]; index2 = reftoshortlistE(edgeshortlist(ii,1),edgeshortlist(ii,2)); cylE1_r2 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)]; end if Ndifforder >= 3 index1 = reftoshortlistE(edgeshortlist(ii,3),edgeshortlist(ii,2)); cylE3_r2 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)]; index2 = reftoshortlistE(edgeshortlist(ii,2),edgeshortlist(ii,3)); cylE2_r3 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)]; end if Ndifforder >= 4 index1 = reftoshortlistE(edgeshortlist(ii,4),edgeshortlist(ii,3)); cylE4_r3 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)]; index2 = reftoshortlistE(edgeshortlist(ii,3),edgeshortlist(ii,4)); cylE3_r4 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)]; end if Ndifforder >= 5 index1 = reftoshortlistE(edgeshortlist(ii,5),edgeshortlist(ii,4)); cylE5_r4 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)]; index2 = reftoshortlistE(edgeshortlist(ii,4),edgeshortlist(ii,5)); cylE4_r5 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)]; end if Ndifforder >= 6 index1 = reftoshortlistE(edgeshortlist(ii,6),edgeshortlist(ii,5)); cylE6_r5 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)]; index2 = reftoshortlistE(edgeshortlist(ii,5),edgeshortlist(ii,6)); cylE5_r6 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)]; end if Ndifforder >= 7 index1 = reftoshortlistE(edgeshortlist(ii,7),edgeshortlist(ii,6)); cylE7_r6 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)]; index2 = reftoshortlistE(edgeshortlist(ii,6),edgeshortlist(ii,7)); cylE6_r7 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)]; end if Ndifforder >= 8 index1 = reftoshortlistE(edgeshortlist(ii,8),edgeshortlist(ii,7)); cylE8_r7 = [re1sho(index1,:) thetae1sho(index1,:) ze1sho(index1,:);re2sho(index1,:) thetae2sho(index1,:) ze2sho(index1,:)]; index2 = reftoshortlistE(edgeshortlist(ii,7),edgeshortlist(ii,8)); cylE7_r8 = [re1sho(index2,:) thetae1sho(index2,:) ze1sho(index2,:);re2sho(index2,:) thetae2sho(index2,:) ze2sho(index2,:)]; end else error('STRANGE TO END UP HERE??') if Ndifforder == 2 cylE1_r2 = [0 0 edgelengthvec(edgeshortlist(ii,jj));0 0 edgelengthvec(edgeshortlist(ii,jj))]; cylE2_r1 = [0 0 edgelengthvec(edgeshortlist(ii,jj));0 0 edgelengthvec(edgeshortlist(ii,jj))]; else error(['ERROR: Geometries with a single edge can not handle difforder >= 3']) end end for jj = 1:ncombs if multfactors(ivlonglist(iv(jj))) if Ndifforder == 2 cylE1_r2frac = cylE1_r2; e1length = edgelengthvec(edgeshortlist(ii,1)); cylE2_r1frac = cylE2_r1; e2length = edgelengthvec(edgeshortlist(ii,2)); if firstedgestart(jj) ~= 0 cylE1_r2frac(1,3) = cylE1_r2frac(1,3) + e1length*firstedgestart(jj); end if firstedgeend(jj) ~= 1 cylE1_r2frac(2,3) = cylE1_r2frac(1,3) + e1length*(1-firstedgeend(jj)); end if lastedgestart(jj) ~= 0 cylE2_r1frac(1,3) = cylE2_r1frac(1,3) + e2length*lastedgestart(jj); end if lastedgeend(jj) ~= 1 cylE2_r1frac(2,3) = cylE2_r1frac(1,3)+ e2length*(1-lastedgeend(jj)); end if midspecpattern(iv(jj))~=0 % For the cases with specular reflections % in-between a double diffraction we must % mirror the two edges to calculate the % relative-to-edge cylindrical coordinates. nspec = sum(midspecpattern(iv(jj),:)>0); % Mirror edge 2 in all the specular reflection % planes, in reversed order edgestartmirr = edgestartcoords(edgeshortlist(ii,2),:); edgeendmirr = edgeendcoords(edgeshortlist(ii,2),:); edgevector = edgeendmirr - edgestartmirr; if lastedgestart(jj) ~= 0 edgestartmirr = edgestartmirr + edgevector*lastedgestart(jj); else edgestartmirr = edgestartmirr + edgevector*1e-6; end if lastedgeend(jj) ~= 1 edgeendmirr = edgestartmirr + edgevector*(1-lastedgeend(jj)); else edgeendmirr = edgestartmirr + edgevector*(1-1e-6); end edgerefcoords = [edgestartcoords(edgeshortlist(ii,1),:);edgeendcoords(edgeshortlist(ii,1),:)]; edgerefnvec = edgenvecs(edgeshortlist(ii,1),:); % If we have a specular reflection in a plane % which is perpendicular to the edge plane, we % should nudge the mirrored edge out a bit so % that there is no 0/(2*pi) mistake if nspec == 1 & ( edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,1)) | edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,2)) ) vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgestartmirr; edgestartmirr = edgestartmirr + vectowardsmidpoint*1e-10; vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgeendmirr; edgeendmirr = edgeendmirr + vectowardsmidpoint*1e-10; end xis = [edgestartmirr;edgeendmirr]; for kk = nspec:-1:1 xis = EDB1findis(xis,[midspecpattern(iv(jj),kk);midspecpattern(iv(jj),kk)],planeeqs,2,onesvec1); end [rstart,thetastart,zstart,rend,thetaend,zend] = EDB1coordtrans2(xis(1,:),xis(2,:),edgerefcoords,edgerefnvec); cylE2mirr_r1 = [rstart thetastart zstart;rend thetaend zend]; % Mirror edge 1 in all the specular reflection % planes, in forward order edgestartmirr = edgestartcoords(edgeshortlist(ii,1),:); edgeendmirr = edgeendcoords(edgeshortlist(ii,1),:); edgevector = edgeendmirr - edgestartmirr; if firstedgestart(jj) ~= 0 edgestartmirr = edgestartmirr + edgevector*firstedgestart(jj); else edgestartmirr = edgestartmirr + edgevector*1e-6; end if firstedgeend(jj) ~= 1 edgeendmirr = edgestartmirr + edgevector*(1-firstedgeend(jj)); else edgeendmirr = edgestartmirr + edgevector*(1-1e-6); end edgerefcoords = [edgestartcoords(edgeshortlist(ii,2),:);edgeendcoords(edgeshortlist(ii,2),:)]; edgerefnvec = edgenvecs(edgeshortlist(ii,2),:); % Normally, when there is a specular reflection % in-between, the diffraction path will not be % along a plane, unless: % 1. The same edge is involved twice, with a % reflection in a perpendicular plane % in-between. % 2. See below pathalongplanewmidspec = 0; if edgeshortlist(ii,1) == edgeshortlist(ii,2) if thetastart == 0 | thetastart == (2*pi-closwedangvec(edgeshortlist(ii,2))) pathalongplanewmidspec = 1; end end % If we have a specular reflection in a plane % which is perpendicular to the edge plane, we % should nudge the mirrored edge out a bit so % that there is no 0/(2*pi) mistake % % We also have case 2 here (see above) for when % we could have a pathalongplanewmidspec % 2. When two different edges have a perpendicular reflection % in-between if nspec == 1 & ( edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,1)) | edgeplaneperptoplane1(midspecpattern(iv(jj),1),edgeshortlist(ii,2)) ) vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgestartmirr; edgestartmirr = edgestartmirr + vectowardsmidpoint*1e-10; vectowardsmidpoint = approxplanemidpoints(midspecpattern(iv(jj),1),:) - edgeendmirr; edgeendmirr = edgeendmirr + vectowardsmidpoint*1e-10; pathalongplanewmidspec = 1; end xis = [edgestartmirr;edgeendmirr]; for kk = 1:nspec xis = EDB1findis(xis,[midspecpattern(iv(jj),kk);midspecpattern(iv(jj),kk)],planeeqs,2,onesvec1); end [rstart,thetastart,zstart,rend,thetaend,zend] = EDB1coordtrans2(xis(1,:),xis(2,:),edgerefcoords,edgerefnvec); cylE1mirr_r2 = [rstart thetastart zstart;rend thetaend zend]; [irnew,slask] = EDB1wedge2nd(cylS(jj,:),cylR(jj,:),cylE2mirr_r1,cylE1mirr_r2,... nyvec(edgeshortlist(ii,:)),[edgelengthvec(edgeshortlist(ii,1))*[firstedgestart(jj) firstedgeend(jj)];edgelengthvec(edgeshortlist(ii,2))*[lastedgestart(jj) lastedgeend(jj)]],dzvec(edgeshortlist(ii,:)),... EDcalcmethod,pathalongplanewmidspec,Rstart,bc,CAIR,FSAMP); else % .... if midspecpattern(iv(jj))~=0 [irnew,slask] = EDB1wedge2nd(cylS(jj,:),cylR(jj,:),cylE2_r1frac,cylE1_r2frac,... nyvec(edgeshortlist(ii,:)),[edgelengthvec(edgeshortlist(ii,1))*[firstedgestart(jj) firstedgeend(jj)];edgelengthvec(edgeshortlist(ii,2))*[lastedgestart(jj) lastedgeend(jj)]],dzvec(edgeshortlist(ii,:)),... EDcalcmethod,pathalongplane,Rstart,bc,CAIR,FSAMP); end % .... if midspecpattern(iv(jj))~=0 if SHOWTEXT >= 6 figure(1) plot(irnew) figure(2) plot(irdiff) pause end IRDIFFVEC = [IRDIFFVEC;sum(irnew)]; elseif Ndifforder == 3, % if Ndifforder == 2 for kk = 1:newndivvec(1) if SHOWTEXT >= 5 disp([' ',int2str(kk),' of ',int2str(newndivvec(1))]) end BIGEDGE1stvalue = (kk-0.5)./newndivvec(1); wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3]; [irnewpartition,slask] = EDB1wedgeN(cylS(jj,:),cylR(jj,:),wedgeparams,ncylrows,... nyvec(edgeshortlist(ii,:)),edgelengthvec(edgeshortlist(ii,:)).',... dzvec(edgeshortlist(ii,:)),EDcalcmethod,pathalongplane,nedgeelcombs,Rstart,bc,CAIR,FSAMP,BIGEDGE1stvalue); irnewpartition = real(irnewpartition); if kk == 1 irnew = irnewpartition; else lengthaddition = length(irnewpartition); lengthaccum = length(irnew); if lengthaddition > lengthaccum irnew = [irnew;zeros(lengthaddition-lengthaccum,1)]; end irnew(1:lengthaddition) = irnew(1:lengthaddition) + irnewpartition; end end if SHOWTEXT >= 6 sum(irnew) plot(irnew) pause end end if Ndifforder == 4 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4]; elseif Ndifforder == 5 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5]; elseif Ndifforder == 6 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6]; elseif Ndifforder == 7 wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6;cylE7_r6;cylE6_r7]; elseif Ndifforder == 8, wedgeparams = [cylE2_r1;cylE1_r2;cylE3_r2;cylE2_r3;cylE4_r3;cylE3_r4;cylE5_r4;cylE4_r5;cylE6_r5;cylE5_r6;cylE7_r6;cylE6_r7;cylE8_r7;cylE7_r8]; end if Ndifforder >= 4, for kk = 1:newndivvec(1) if SHOWTEXT >= 5 disp([' ',int2str(kk),' of ',int2str(newndivvec(1))]) end BIGEDGE1stvalue = (kk-0.5)./newndivvec(1); [irnewpartition,slask] = EDB1wedgeN(cylS(jj,:),cylR(jj,:),wedgeparams,ncylrows,... nyvec(edgeshortlist(ii,:)),edgelengthvec(edgeshortlist(ii,:)).',... dzvec(edgeshortlist(ii,:)),EDcalcmethod,pathalongplane,nedgeelcombs,Rstart,bc,CAIR,FSAMP,BIGEDGE1stvalue); irnewpartition = real(irnewpartition); if kk == 1 irnew = irnewpartition; else lengthaddition = length(irnewpartition); lengthaccum = length(irnew); if lengthaddition > lengthaccum irnew = [irnew;zeros(lengthaddition-lengthaccum,1)]; end irnew(1:lengthaddition) = irnew(1:lengthaddition) + irnewpartition; end end end if SHOWTEXT >= 5 varname = ['allirs',int2str(Ndifforder)]; if exist(varname) == 0 eval([varname,' = irnew;']) else lengthaddition = length(irnew); eval(['lengthaccum = size(',varname,',1);']) if lengthaddition > lengthaccum eval([varname,' = [',varname,';zeros(lengthaddition-lengthaccum,size(',varname,',2))];']) end eval([varname,' = [',varname,' zeros(size(',varname,',1),1)];']) eval([varname,'(1:length(irnew),end) = irnew;']) end end % Decide whether the IR should be boosted or not boostfactor = 1; if souspecboost ~= 1 | recspecboost ~= 1 if prespecpresent(iv(jj)) == 1 boostfactor = souspecboost; else if ~isempty(Sinsideplanenumber) if Sinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,1),1) & Sinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,1),2), boostfactor = souspecboost; end end end if postspecpresent(iv(jj)) == 1 boostfactor = boostfactor*recspecboost; else if ~isempty(Rinsideplanenumber) if Rinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,Ndifforder),1) & Rinsideplanenumber(1)~=planesatedge(edgeshortlist(ii,Ndifforder),2), boostfactor = boostfactor*recspecboost; end end end end % ... if souspecboost ~= 1 | recspecboost ~= 1 % For thin plates, we must have a boost factor! % This is because there will be multiple equivalent % combinations passing on the rear side of the thin plate if all( nyvec(edgeshortlist(ii,:)) == 0.5 ) boostfactor = boostfactor*2^(Ndifforder-1); end if boostfactor ~= 1 irnew = irnew*boostfactor; end irnew = irnew*double(multfactors(ivlonglist(iv(jj)))); [ndiff,ncolsdiff] = size(irdiff); nnew = size(irnew,1); if saveindividualdiffirs(1) == 0 if nnew > ndiff irdiff = [irdiff;zeros(nnew-ndiff,1)]; end irdiff(1:nnew) = irdiff(1:nnew) + irnew; else if Ndifforder > ncolsdiff irdiff = [irdiff zeros(ndiff,Ndifforder-ncolsdiff)]; ncolsdiff = ncolsdiff+1; end if nnew > ndiff irdiff = [irdiff;zeros(nnew-ndiff,Ndifforder)]; end irdiff(1:nnew,Ndifforder) = irdiff(1:nnew,Ndifforder) + irnew; end end % ... if multfactors(ivlonglist(iv(jj))) end % ....for jj = 1:ncombs else % ... if any(multfactors(ivlonglist(iv))) if SHOWTEXT >= 4 disp([' Combination not computed because of repetition']) end end % ... if any(multfactors(ivlonglist(iv))) end % .... for ii = 1: length(ivec) end % ... if any(difforderinlist==Ndifforder) if size(irdiff,2) < Ndifforder & saveindividualdiffirs(1) == 1 irdiff = [irdiff zeros(size(irdiff,1),Ndifforder-size(irdiff,2))]; end end % ... for Ndifforder = 2:specorder end % .... if userwantsdiff2 == 1 ntot = length(irtot); ndirect = length(irdirect); ngeom = size(irgeom,1); ndiff = size(irdiff,1); if ndirect > ntot irtot = [irtot;zeros(ndirect-ntot,1)]; ntot = length(irtot); end irtot(1:ndirect) = irtot(1:ndirect) + irdirect; if ngeom > ntot irtot = [irtot;zeros(ngeom-ntot,1)]; ntot = length(irtot); end irtot(1:ngeom) = irtot(1:ngeom) + irgeom; if ndiff > ntot irtot = [irtot;zeros(ndiff-ntot,1)]; ntot = length(irtot); end if saveindividualdiffirs(1) == 0 | userwantsdiff2 == 0 irtot(1:ndiff) = irtot(1:ndiff) + irdiff; else irtot(1:ndiff) = irtot(1:ndiff) + sum(irdiff.').'; end %------------------------------------------------------- % Make the IRs the same length nmax = max([ length(irtot) size(irgeom,1) size(irdiff,1) length(irdirect)]); if length(irtot) < nmax irtot = [irtot;zeros(nmax-length(irtot),1)]; end if length(irdirect) < nmax irdirect = [irdirect;zeros(nmax-length(irdirect),1)]; end if length(irgeom) < nmax irgeom = [irgeom;zeros(nmax-size(irgeom,1),size(irgeom,2))]; end if length(irdiff) < nmax irdiff = [irdiff;zeros(nmax-size(irdiff,1),size(irdiff,2))]; end %------------------------------------------------------- % Save the IRs irtot = sparse(irtot); irgeom = sparse(irgeom); irdirect = sparse(irdirect); irdiff = sparse(irdiff); Varlist = [' irdirect irgeom irdiff irtot FSAMP Rstart CAIR EDcalcmethod']; if length(saveindividualdiffirs) > 1 if saveindividualdiffirs(2) == 1 Varlist = [Varlist,' alldiffirs alldiffirsextradata']; end end if SHOWTEXT >= 5 if specorder >= 2 if exist('allirs2') == 0 allirs2 = []; end Varlist = [Varlist,' allirs2']; end if specorder >= 3 if exist('allirs3') == 0 allirs3 = []; end Varlist = [Varlist,' allirs3']; end if specorder >= 4 if exist('allirs4') == 0 allirs4 = []; end Varlist = [Varlist,' allirs4']; end if specorder >= 5 if exist('allirs5') == 0 allirs5 = []; end Varlist = [Varlist,' allirs5']; end if specorder >= 6 if exist('allirs6') == 0 allirs6 = []; end Varlist = [Varlist,' allirs6']; end if specorder >= 7 if exist('allirs7') == 0 allirs7 = []; end Varlist = [Varlist,' allirs7']; end if specorder >= 8 if exist('allirs8') == 0 allirs8 = []; end Varlist = [Varlist,' allirs8']; end end eval(['save ',edirfile,Varlist])