tp@0: function [outputfile] = EDB1readcad(CADfile,outputfile,planecornerstype,checkgeom)
tp@0: % EDB1readcad - Reads a file of type .CAD (made by e.g. CATT-Acoustic) and saves all the geometry data in a mat-file. 
tp@0: % CAD-files v6,v7,v8 are read.
tp@0: %
tp@0: % Input parameters:
tp@0: %  CADfile     	(optional) The input file, with or without the .CAD extension.
tp@0: %					If this file is not specified, a file opening window will appear.
tp@0: %  outputfile      (optional) The output file. If not specified, it will be given the
tp@0: %                   same name as the CAD file with an '_cadgeo.mat' ending instead of the '.CAD'.
tp@0: %  planecornerstype (optional) Could have the value 'zero' or 'circ'.Default: 'circ'.
tp@0: %               	Affects the matrix planecorners, see below.
tp@0: %  checkgeom		(optional) If this parameter is given the value 'check', then a few checks
tp@0: %                   of the geometry consistency will be done: a check for duplicate corners
tp@0: %                   for redundant corners and for corners that are connected to only one plane.
tp@0: %					Only warnings are given. As default, no check is done.
tp@0: %  SHOWTEXT (global)	If this global parameter has the value 3 or higher, informative text will
tp@0: %                   be printed on the screen.
tp@0: %
tp@0: % Output parameters:
tp@0: %	outputfile		The name of the outputfile, generated either automatically or specified as 
tp@0: %					an input parameter.
tp@0: %
tp@0: % Output parameters in the outputfile:
tp@0: % 	corners         Matrix [ncorners,3] with the corner coordinates
tp@0: %   cornernumbers   Vector [ncorners,1] with the corner numbers used in the CAD-file. That is
tp@0: %                   in the outputfile, the corners will effectively be renumbered from 1 to ncorners
tp@0: %                   in the order they appeared in the CAD-file. 
tp@0: %   planecorners    Matrix [planes,nmaxcornersperplane] with the corner numbers that make up each plane.
tp@0: %                   Since planes can have different numbers of corners, the number of columns is the 
tp@0: %                   maximum number of corners that any plane has. NB! The corner numbers are in the
tp@0: %                   renumbered system, not in the CAD-file system.
tp@0: %   planenames      Matrix [nplanes,nmaxcharacters1] (sparse), with the planenames in the CAD file.
tp@0: %   planeabstypes   Matrix [nplanes,nmaxcharacters2] (sparse), with the absorber names in the CAD file. 
tp@0: % 	planenumbers    Vector [nplanes,1] with the plane numbers used in the CAD file.
tp@0: %   cornercrossref  Vector [nmaxnumberinCADfile,1] (sparse), which gives the matlab-file corner
tp@0: %                   numbers for each CAD-file corner number. For instance, if a CAD-file corner
tp@0: %                   number 1200 is translated to corner number 72 in the matlab-file, then
tp@0: %                   cornercrossref(1000) = 72.
tp@0: %   Snumbers        Vector [nsources,1] of the source numbers in the CAD-file. NB! Only CAD-files
tp@0: %                   v6 use source numbers; later versions use source names instead of numbers. 
tp@0: %                   For CAD-file v7 and v8, Snumbers is empty.
tp@0: %   Snames          Matrix [nsources,2] of the source names in the CAD-file. NB! Only CAD-files v7
tp@0: %                   or later use source names (such as 'A0' etc). Older versions use source numbers, 
tp@0: %                   in which case Snames is empty.
tp@0: %   Sdirectivitynames   Matrix [nsources,nmaxcharacters3] (sparse) with the source directivity names 
tp@0: %                   in the CAD file. 
tp@0: %   Sdelays         Vector [nsources,1] of the source delays in the CAD-file, in milliseconds.
tp@0: % 	Scoords         Matrix [nsources,3] of the source coordinates.
tp@0: %   Sdirections     Matrix [nsources,3] of the aim point coordinates for the sources.
tp@0: %   Srotations      Vector [nsources,1] of the source rotations, in degrees.
tp@0: %   Slevels         Matrix [nsources,6/8] of the source levels for six (CAD-files v6) or eight
tp@0: %                   octave bands (CAD-files v7 or later).
tp@0: %   Rnumbers        Vector [nreceivers,1] of the receiver numbers in the CAD-file.
tp@0: %   Rcoords         Matrix [nreceivers,3] of the receiver coordinates.
tp@0: %   CATTversionnumber   6,7 or 8
tp@0: %	planeeqs        Matrix [nplanes,4] of the plane equations as derived  from the plane definitions. 
tp@0: %                   Each row has the values [A B C D] of the plane equation on the form Ax + By + Cz = D
tp@0: %   planenvecs      Matrix [nplanes,3] of the normalized normal vectors of the planes.
tp@0: %   ncornersperplanevec     Vector [nplanes,1] which gives the number of corners for each plane.
tp@0: %   planesatcorners     Matrix [ncorners,4] which gives the plane numbers that each corner is connected to.
tp@0: %                   NB! Here it is assumed that each corner is connected to maximum four planes.
tp@0: %   nplanespercorners   Vector [ncorners,1] which gives the number of planes that each corner is connected to.
tp@0: % 	minvals         Matrix [nplanes,3] which for each plane gives the smallest coordinate in the x-, y- and z-direction.
tp@0: % 	maxvals         Matrix [nplanes,3] which for each plane gives the largest coordinate in the x-, y- and z-direction.
tp@0: %   planehasindents Vector [nplanes,1] which for each plane gives the
tp@0: %                   number of indeting corners
tp@0: %   indentingcorners  Matrix [nplanes,max(ncornersperplanevec)] which for each plane gives the number to the first corner
tp@0: %                   in a corner triplet which identifies an indenting corner. The number of the corner is the
tp@0: %                   order given in planecorners for that plane.
tp@0: %
tp@0: % Uses the functions EDB1extrnums, EDB1cross
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) 20090718
tp@0: %
tp@0: % [outputfile] = EDB1readcad(CADfile,outputfile,planecornerstype,checkgeom);
tp@0: 
tp@0: global SHOWTEXT
tp@0: 
tp@0: if nargin == 0
tp@0: 	CADfile = '';
tp@0: 	outputfile = '';
tp@0: 	planecornerstype = '';
tp@0: 	checkgeom = '';
tp@0: elseif nargin == 1
tp@0: 	outputfile = '';
tp@0: 	planecornerstype = '';
tp@0: 	checkgeom = '';
tp@0: elseif nargin == 2
tp@0: 	planecornerstype = '';
tp@0: 	checkgeom = '';
tp@0: elseif nargin == 3
tp@0: 	checkgeom = '';
tp@0: end
tp@0: 
tp@0: geomacc = 1e-10;
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % If no CAD-file was specified, present a file opening window
tp@0: 
tp@0: if isempty(CADfile)
tp@0: 	[CADfile,filepath] = uigetfile('*.cad','Please select the CADfile');
tp@0:     [filepath,temp1,temp2] = fileparts(filepath);
tp@0: 	if ~isstr(CADfile)
tp@0: 		return
tp@0: 	end
tp@0:     [temp1,filestem,temp2] = fileparts(CADfile);
tp@0: 
tp@0:     CADfile = [[filepath,filesep],filestem,'.cad'];
tp@0: else
tp@0: 	[filepath,filestem,fileext] = fileparts(CADfile);
tp@0:     CADfile = [[filepath,filesep],filestem,fileext];
tp@0: end
tp@0: 
tp@0: if exist(CADfile) ~= 2
tp@0: 	error(['ERROR: CAD-file: ',CADfile,' can not be found'])
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % If no output file was specified, construct an automatic file name
tp@0: 
tp@0: if isempty(outputfile)
tp@0: 	outputfile = [[filepath,filesep],filestem,'_cadgeo.mat'];
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Read in the entire file into the string B
tp@0: % Find the line starts and ends
tp@0: 
tp@0: fid = fopen([CADfile],'r');
tp@0: if fid == 0
tp@0: 	error(['ERROR: The CADfile ',CADfile,' could not be opened.'])
tp@0: end
tp@0: B = fread(fid,inf,'char').';
tp@0: fclose(fid);
tp@0: 
tp@0: % Cut out sequential spaces to save some space.
tp@0: iv = find(B(1:length(B)-1)==32 & B(2:length(B))==32);
tp@0: B(iv) = [];
tp@0: 
tp@0: % If the text file was generated on a Mac, each line has ASCII 13 at the end.
tp@0: % On a PC, each line ends with ASCII 13, and the next line starts with ASCII 10.
tp@0: % In Unix, each line ends with ASCII 10.
tp@0: %
tp@0: % Change this to a single ASCII 13
tp@0: 
tp@0: % New 17 Jul 09: use the builtin function regexprep instead:
tp@0: 
tp@0: % 1. Replace all ASCII 10 by ASCII 13
tp@0: 
tp@0: B = regexprep(setstr(B),'\n','\r');
tp@0: 
tp@0: % 2. Cut out sequential CRs (ASCII13) to save some space.
tp@0: iv = find(B(1:length(B)-1)==13 & B(2:length(B))==13);
tp@0: B(iv) = [];
tp@0: 
tp@0: % 3. Cut out sequences of CR, space, CR because they "survived" the
tp@0: %    search for consecutive blanks.
tp@0: iv = find(B(1:length(B)-1)==32 & B(2:length(B))==13);
tp@0: B(iv) = [];
tp@0: iv = find(B(1:length(B)-1)==13 & B(2:length(B))==13);
tp@0: B(iv) = [];
tp@0: 
tp@0: % 4. Convert all text to lower case
tp@0: B = lower(setstr(B));
tp@0: 
tp@0: % Some special characters might give negative ASCII values
tp@0: % (Comment by PS 17Jul09: Will this ever happen??)
tp@0: 
tp@0: iv = find(B<0);
tp@0: if ~isempty(iv)
tp@0: 	B(iv) = B(iv)+256;   
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Look for the text 'CATT-Acoustic v' anywhere, and read the numerical
tp@0: % value right after this text.
tp@0: 
tp@0: iv = regexp(B,'catt-acoustic v');
tp@0: 
tp@0: if isempty(iv)
tp@0: 	CATTversionnumber = [];
tp@0: else
tp@0:  	CATTversionnumber = str2num(setstr(B(iv(1)+15)));
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Find the four sections of the file:
tp@0: %    %CORNERS      %PLANES      %SOURCES   %RECEIVERS
tp@0: % (or small letters)
tp@0: % We assume that they come in this order
tp@0: 
tp@0: stringposCORNERS = regexp(B,'%corners');
tp@0: if isempty(stringposCORNERS)
tp@0:  	   error('ERROR: The .cad file must contain the word %CORNERS');    
tp@0: end
tp@0: 
tp@0: stringposPLANES = regexp(B,'%planes');
tp@0: if isempty(stringposPLANES)
tp@0:  	   error('ERROR: The .cad file must contain the word %PLANES');    
tp@0: end
tp@0: 
tp@0: stringposSOURCES = regexp(B,'%sources');
tp@0: 
tp@0: stringposRECEIVERS = regexp(B,'%receivers');
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Look for the corners definitions between the lines containing
tp@0: % %SOURCES and another keyword  (%PLANES or %SOURCES or %RECEIVERS)
tp@0: % The lines should look like:
tp@0: %	1   -0.2000   0.34   2.3
tp@0: 
tp@0: if stringposPLANES < stringposCORNERS
tp@0:     error('ERROR: Sorry, but in the .cad file, the %CORNERS section must come before the %PLANES section')
tp@0: end
tp@0: 
tp@0: C = textscan(B(stringposCORNERS:stringposPLANES-1),'%d%f%f%f','Headerlines',1);
tp@0: cornernumbers = (C{1});
tp@0: corners = [C{2} C{3} C{4}];
tp@0: ncorners = length(cornernumbers);
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Look for the planes definitions between the lines containing
tp@0: % %PLANES and %SOURCES, or the end
tp@0: % The lines should look like:
tp@0: %  1 / planename /RIGID
tp@0: %  1  4  3  2
tp@0: %  2 / /RIGID
tp@0: 
tp@0: if ~isempty(stringposSOURCES)
tp@0:     if stringposSOURCES < stringposPLANES
tp@0:         error('ERROR: Sorry, but in the .cad file, the %PLANES section must come before the %SOURCES section')
tp@0:     end    
tp@0:     Str = B(stringposPLANES:stringposSOURCES-1);
tp@0: else
tp@0:     Str = B(stringposPLANES:end);    
tp@0: end
tp@0: 
tp@0: % To make it easier with textscan, we put the two types of line onto
tp@0: % a single line, with a new '/' inserted. Then all lines have the same
tp@0: % format. 
tp@0: 
tp@0: % Find all CR and replace every second with a '/'
tp@0: ivendlines = find(Str==13);
tp@0: Str(ivendlines(2:2:end)) = '/';
tp@0: 
tp@0: C = textscan(Str,'%d%s%s%s','Headerlines',1,'Delimiter','/','BufSize',32768);
tp@0: 
tp@0: planenumbers = C{1};
tp@0: Str2 = C{2};
tp@0: Str3 = C{3};
tp@0: Str = C{4};
tp@0: 
tp@0: nplanes = size(Str,1);
tp@0: ncornersperplanevec = zeros(nplanes,1);
tp@0: planecorners = zeros(nplanes,4);
tp@0: planenames = zeros(nplanes,30);
tp@0: planeabstypes = zeros(nplanes,6);
tp@0: 
tp@0: for ii = 1:nplanes
tp@0:     listofplanes = EDB1extrnums(Str{ii});
tp@0:     ncornersperplanevec(ii) = length(listofplanes);
tp@0:     if ncornersperplanevec(ii) > size(planecorners,2)
tp@0:        planecorners = [planecorners zeros(nplanes, ncornersperplanevec(ii)-size(planecorners,2))];
tp@0:     end
tp@0:     planecorners(ii,1:ncornersperplanevec(ii)) = listofplanes;
tp@0: 
tp@0:     txtstr2 = Str2{ii};
tp@0:     if length(txtstr2) > size(planenames,2)
tp@0:        planenames = [planenames zeros(nplanes,length(txtstr2)-size(planenames,2))]; 
tp@0:     end
tp@0:     planenames(ii,1:length(txtstr2)) = txtstr2;
tp@0: 
tp@0:     txtstr3 = Str3{ii};
tp@0:     if length(txtstr3) > size(planeabstypes,2)
tp@0:        planeabstypes = [planeabstypes zeros(nplanes,length(txtstr3)-size(planeabstypes,2))]; 
tp@0:     end
tp@0:     planeabstypes(ii,1:length(txtstr3)) = txtstr3;
tp@0: end
tp@0: 
tp@0: if max(max(planecorners)) > max(cornernumbers)
tp@0:     error('ERROR: One plane definition in the CAD-file used a higher corner number than was defined in the CORNERS section')
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Look for the sources definitions between the lines containing
tp@0: % %SOURCES and %RECEIVERS, or the end
tp@0: % The lines should look like: (v6)
tp@0: % 0 OMNI:SD0
tp@0: %   0.0000000   0.0000000   1.0000000
tp@0: %   0.0000000   0.0000000   0.0000000
tp@0: %  85.0  88.0  91.0  94.0  97.0  100.0
tp@0: %
tp@0: % or: (v7)
tp@0: % A0 OMNI:SD0
tp@0: %   0.0000000   0.0000000   1.0000000
tp@0: %   0.0000000   0.0000000   0.0000000  15.00000
tp@0: %  85.0  88.0  91.0  94.0  97.0  100.0 : 103.0 106.0
tp@0: 
tp@0: if ~isempty(stringposSOURCES)
tp@0: 
tp@0:     if ~isempty(stringposRECEIVERS)
tp@0:         if stringposRECEIVERS < stringposSOURCES
tp@0:             error('ERROR: Sorry, but in the .cad file, the %SOURCES section must come before the %RECEIVERS section')
tp@0:         end    
tp@0:         Str = B(stringposSOURCES:stringposRECEIVERS-1);
tp@0:     else
tp@0:         Str = B(stringposSOURCES:end);    
tp@0:     end
tp@0: 
tp@0:     % To make it easier with textscan, we put the four types of line onto
tp@0:     % a single line, with new '/' inserted. Then all lines have the same
tp@0:     % format. 
tp@0:     % Find all CR and replace number 2,3,4, 6,7,8, etc with a '/'
tp@0:     
tp@0:     ivendlines = find(Str==13);
tp@0:     ivkeep = [1:4:length(ivendlines)];
tp@0:     ivremove = [1:length(ivendlines)];
tp@0:     ivremove(ivkeep) = [];
tp@0:     Str(ivendlines(ivremove)) = '/';
tp@0: 
tp@0:     C = textscan(Str,'%s%s%s%s','Headerlines',1,'Delimiter','/');
tp@0: 
tp@0:     nsources = size( C{1},1 );
tp@0:     Scoords = zeros(nsources,3);
tp@0:     Sdirections = zeros(nsources,3);
tp@0:     Srotations = zeros(nsources,1);
tp@0:     Slevels = zeros(nsources,8);
tp@0:     Snumbers = zeros(nsources,1);
tp@0:     Snames = zeros(nsources,4);
tp@0:     Sdirectivitynames = ( zeros(nsources,30) );
tp@0:     Sdelays = zeros(nsources,1);
tp@0: 
tp@0:     Str1 = C{1};
tp@0:     Str2 = C{2};
tp@0:     Str3 = C{3};
tp@0:     Str4 = C{4};
tp@0:     for ii = 1:nsources    
tp@0:         Scoords(ii,:) = EDB1extrnums(Str2{ii});
tp@0:         listofvalues3 = EDB1extrnums(Str3{ii});
tp@0:         listofvalues4 = EDB1extrnums(Str4{ii});
tp@0:         if ii==1
tp@0:            if length(listofvalues3)==3
tp@0:                CATTversionnumber = 6;
tp@0:            elseif length(listofvalues3)==4
tp@0:                CATTversionnumber = 7;
tp@0:            else
tp@0:                error('ERROR: The source data should have three or four values for the direction')
tp@0:            end
tp@0:         end
tp@0: 
tp@0:         Sdirections(ii,1:3) = listofvalues3(1:3);
tp@0: 
tp@0:         % The first field contains three items which are a mix of text and numerical values.
tp@0:         % First we remove possible consecutive blanks.
tp@0:         textstr = Str1{ii};
tp@0:         iv = find(textstr(1:length(textstr)-1)==32 & textstr(2:length(textstr))==32);
tp@0:         textstr(iv) = [];
tp@0:         iv = find(textstr==32);
tp@0: 
tp@0: 
tp@0:         if CATTversionnumber == 7
tp@0:            Srotations(ii) = listofvalues3(4); 
tp@0:            Snames(ii,1:length(textstr(1:iv(1)-1))) = textstr(1:iv(1)-1);
tp@0:             Sdirectivitynames(ii,1:iv(2)-iv(1)-1) = textstr(iv(1)+1:iv(2)-1);
tp@0:             Sdelays(ii) = str2num( textstr(iv(2)+1:end) );
tp@0:         else
tp@0:             Snumbers(ii) = str2num( textstr(1:iv(1)-1) );
tp@0:             Sdirectivitynames(ii,1:length(textstr(iv(1)+1:end))) = textstr(iv(1)+1:end);
tp@0:         end
tp@0:         Slevels(ii,1:length(listofvalues4)) = listofvalues4;      
tp@0:     end
tp@0:     Snames = setstr(Snames);
tp@0:     Sdirectivitynames = setstr(Sdirectivitynames);
tp@0: 
tp@0: else
tp@0:     
tp@0:     Scoords = [];
tp@0:     Sdirections = [];
tp@0:     Srotations = [];
tp@0:     Slevels = [];
tp@0:     Snumbers = [];
tp@0:     Snames = [];
tp@0:     Sdirectivitynames = [];
tp@0:     Sdelays = [];    
tp@0:     
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Look for the receivers definitions after the line containing
tp@0: % %RECEIVERS
tp@0: % The lines should look like:
tp@0: %	1   -0.2000   0.34   2.3
tp@0: 
tp@0: if ~isempty(stringposRECEIVERS)
tp@0:     C = textscan(B(stringposRECEIVERS:end),'%d%f%f%f','Headerlines',1);
tp@0:     Rnumbers = (C{1});
tp@0:     Rcoords= [C{2} C{3} C{4}];
tp@0: else
tp@0:     Rnumbers = [];
tp@0:     Rcoords = [];
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Change the numbering in the CAD file to a contiguous one
tp@0: 
tp@0: % Corner numbers
tp@0: % First we make an inelegant crossreference vector giving the
tp@0: % resulting corner number in the position given by the CAD-file
tp@0: % corner number.
tp@0: 
tp@0: maxnumb = max(cornernumbers);
tp@0: cornercrossref = sparse(zeros(1,maxnumb));
tp@0: 
tp@0: % Temporary fix: if one corner has the number zero, then we make a 
tp@0: % separate crossref variable for that one.
tp@0: if cornernumbers(1) == 0
tp@0:     cornercrossrefzero = 1;
tp@0:     cornercrossref(cornernumbers(2:end)) = [2:ncorners];
tp@0: else
tp@0:     cornercrossrefzero = [];
tp@0:     cornercrossref(cornernumbers) = [1:ncorners];
tp@0: end
tp@0: 
tp@0: [nplanes,ncornersperplane] = size(planecorners);
tp@0: 
tp@0: planecorners_CATTnumbers = planecorners;
tp@0: 
tp@0: iv = find(planecorners_CATTnumbers~=0);
tp@0: planecorners(iv) = full(cornercrossref(planecorners(iv)));
tp@0: iv = find(planecorners_CATTnumbers==0);
tp@0: planecorners(iv) = cornercrossrefzero;
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Go through all planes. If there is a plane definition including
tp@0: % zeros, and planecornerstype == 'circ', expand it repeating the
tp@0: % same corner order again.
tp@0: 
tp@0: if isempty(planecornerstype)
tp@0: 	planecornerstype = 'circ';
tp@0: else
tp@0: 	planecornerstype = setstr(lower(planecornerstype(1)));
tp@0: 	if planecornerstype(1) == 'z'
tp@0: 		planecornerstype = 'zero';
tp@0: 	else
tp@0: 		planecornerstype = 'circ';
tp@0: 	end
tp@0: end
tp@0: 
tp@0: if planecornerstype == 'circ'
tp@0: 	for ii = 1:nplanes
tp@0: 		iv = find( planecorners(ii,:) ~= 0);
tp@0: 		ncornersatplane = length(iv);
tp@0: 		if ncornersatplane ~= ncornersperplane
tp@0: 			pattern = planecorners(ii,iv);
tp@0: 			nrepeatings = ceil(ncornersperplane/ncornersatplane);
tp@0: 			for jj = 1:nrepeatings-1
tp@0: 				pattern = [pattern planecorners(ii,iv)];
tp@0: 			end
tp@0: 			planecorners(ii,:) = pattern(1:ncornersperplane);
tp@0: 		end
tp@0: 	end
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Find the normal vectors to the planes using the cross products
tp@0: %
tp@0: % The normal vector is basically the cross product between two vectors
tp@0: % connecting three consecutive points. If there are indents though
tp@0: % they will cause reversed normal vectors, so one must go through all
tp@0: % combinations and choose the majority normal vector.
tp@0: %
tp@0: % 26mar09  Use the fact described above for detecting indention corners
tp@0: 
tp@0: planenvecs = zeros(nplanes,3);
tp@0: planehasindents = zeros(nplanes,1);
tp@0: indentingcorners = sparse(zeros(nplanes,max(ncornersperplanevec)));
tp@0: 
tp@0: for ii = 1:nplanes
tp@0:     
tp@0: 	iv = find(planecorners(ii,:)~=0);
tp@0: 	cornerlist = planecorners(ii,iv);
tp@0: 	iv = find(cornerlist == cornerlist(1));
tp@0: 	if length(iv) > 1
tp@0: 		cornerlist = cornerlist(1:iv(2)-1);
tp@0: 	end
tp@0: 	ncorners = length( cornerlist );
tp@0: 	cornerlist = [cornerlist cornerlist(1) cornerlist(2)];
tp@0: 
tp@0: 	nvectorlist = zeros(ncorners,3);
tp@0: 	nveclen = zeros(ncorners,1);	
tp@0: 	
tp@0: 	for jj = 1:ncorners
tp@0: 		co1numb = cornerlist(jj);
tp@0: 		co2numb = cornerlist(jj+1);
tp@0: 		co3numb = cornerlist(jj+2);
tp@0: 		vec1 = corners(co1numb,:) - corners(co2numb,:);
tp@0: 		vec2 = corners(co3numb,:) - corners(co2numb,:);
tp@0: 		nvec = EDB1cross(vec1.',vec2.').';
tp@0: 		nveclen(jj) = norm(nvec);
tp@0: 		if nveclen(jj) > 0
tp@0: 			nvectorlist(jj,:) = nvec./nveclen(jj);
tp@0: 		end
tp@0: 	end
tp@0: 	
tp@0: 	iv = find(nveclen < max(nveclen)*0.001);
tp@0: 	nvectorlist(iv,:) = [];
tp@0: 	nvecref = nvectorlist(1,:);
tp@0: 	
tp@0:     [n1,slask] = size(nvectorlist);
tp@0: 	nvecsigns = round(sum(     (nvectorlist.').*(nvecref(ones(n1,1),:).')      ));
tp@0:     
tp@0:     if sum(nvecsigns) == 0
tp@0:         disp(' ')        
tp@0:        error(['ERROR: Plane ',int2str(planenumbers(ii)),' (plane numbering as in the CAD file) seems to be twisted.'])        
tp@0:     end
tp@0:     
tp@0:     
tp@0:     if abs(sum(nvecsigns)) ~= n1
tp@0:        disp(['Plane ',int2str(ii),' has indents!']) 
tp@0:        nindents = (n1 - abs(sum(nvecsigns)))/2;
tp@0:        planehasindents(ii) = nindents;
tp@0:        ivindent = find(nvecsigns == -1);
tp@0:        if length(ivindent) == nindents
tp@0:             indentingcorners(ii,ivindent) = 1;
tp@0:        else
tp@0:            if length(ivindent) == ncorners - nindents
tp@0:                ivindent = find(nvecsigns == 1);
tp@0:                 indentingcorners(ii,ivindent) = 1;           
tp@0:            else
tp@0:               error(['ERROR: An unexpected problem. Please report to the developer']) 
tp@0:            end
tp@0:        end
tp@0:     end
tp@0:     
tp@0:     nvecdiff = [nvectorlist(2:n1,1).*nvecsigns(2:n1).' nvectorlist(2:n1,2).*nvecsigns(2:n1).' nvectorlist(2:n1,3).*nvecsigns(2:n1).'] - nvecref(ones(n1-1,1),:);
tp@0:         
tp@0:     if n1 > 2
tp@0:         nvecdiff = sum(nvecdiff.'.^2).';    
tp@0:     else
tp@0:         nvecdiff = norm(nvecdiff);
tp@0:     end
tp@0:     if any(nvecdiff>1e-4),        
tp@0:         nvecdiff
tp@0:         error(['ERROR: Normal vectors for plane ',int2str(planenumbers(ii)),' (in the CAD file, = ',int2str(ii),' in the EDB1 file), get different normal vectors for consecutive corner triplets. Check the geometry in the CAD-file'])
tp@0:     elseif any(nvecdiff>1e-8)
tp@0:         nvecdiff
tp@0:         disp(['WARNING: Normal vectors for plane ',int2str(planenumbers(ii)),' (in the CAD file, = ',int2str(ii),' in the EDB1 file), get somewhat different normal vectors for consecutive corner triplets. Check the geometry in the CAD-file'])
tp@0:     end
tp@0:     
tp@0: 	if ncorners > 5 & abs(sum(nvecsigns)) <= 1
tp@0: 		disp(['WARNING for plane number ',int2str(planenumbers(ii)),' in the CAD-file'])
tp@0: 		disp(['   with the name ',strtrim(setstr(full(planenames(ii,:))))])
tp@0: 		disp('   The normal vector can not be determined for this plane because there are')
tp@0: 		disp('   the same number of inwards and outwards corners')
tp@0: 		disp('   This is a list of the possible normal vectors:')
tp@0: 		[nv1,slask] = size(nvectorlist);
tp@0: 		for kk = 1:nv1
tp@0: 			vecstr = ['   ',int2str(kk),'. ',num2str(-nvectorlist(kk,1)),' ',num2str(-nvectorlist(kk,2)),' ',num2str(-nvectorlist(kk,3))];
tp@0: 			disp(vecstr)
tp@0: 		end
tp@0: 		disp(' ')
tp@0: 	
tp@0:       preferredsign = input('   Give the number of a correct normal vector for this plane please ');
tp@0:       switchsign = nvecsigns.'./nvecsigns(preferredsign);
tp@0: 		nvectorlist = nvectorlist.*switchsign(:,ones(1,3));	
tp@0: 	else
tp@0: 		mostcommonsign = sign(sum(nvecsigns));
tp@0: 
tp@0: 		switchsign = nvecsigns.'./mostcommonsign;
tp@0: 		nvectorlist = nvectorlist.*switchsign(:,ones(1,3));
tp@0:     end
tp@0: 
tp@0: 	planenvecs(ii,:) = mean(nvectorlist);
tp@0: end
tp@0: 
tp@0: planenvecs = -planenvecs;
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Plane equations, as Ax + By + Cz = D for each plane
tp@0: 
tp@0: planeeqs = zeros(nplanes,4);
tp@0: planeeqs(:,1:3) = planenvecs;
tp@0: planeeqs(:,4) = sum( (planenvecs.').*(corners(planecorners(:,1),:).')  ).';
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Useful data: planesatcorners, minvals and maxvals
tp@0: 
tp@0: [ncorners,slask] = size(corners);
tp@0: planesatcornerhits = zeros(ncorners,nplanes);
tp@0: 
tp@0: for ii = 1:nplanes
tp@0: 	cornerlist = planecorners(ii,1:ncornersperplanevec(ii));
tp@0: 	planesatcornerhits(cornerlist,ii) = planesatcornerhits(cornerlist,ii) + 1;
tp@0: end
tp@0: 
tp@0: maxplanespercorner = 0;
tp@0: for ii = 1:ncorners
tp@0: 	nplanes = length(find(planesatcornerhits(ii,:) ~= 0));
tp@0: 	if nplanes > maxplanespercorner
tp@0: 		maxplanespercorner = nplanes;
tp@0: 	end	
tp@0: end
tp@0: 
tp@0: planesatcorners = zeros(ncorners,maxplanespercorner);
tp@0: nplanespercorners = zeros(ncorners,1);
tp@0: for ii = 1:ncorners
tp@0: 	iv = find(planesatcornerhits(ii,:)~=0);
tp@0: 	planesatcorners(ii,1:length(iv)) = iv;
tp@0: 	nplanespercorners(ii) = length(iv);
tp@0: end
tp@0: 
tp@0: % find cubic boxes inside which the planes are placed
tp@0: 
tp@0: [nplanes,slask] = size(planeeqs);
tp@0: 
tp@0: minvals = zeros(nplanes,3);
tp@0: maxvals = zeros(nplanes,3);
tp@0: 
tp@0: for ii = 1:nplanes
tp@0: 	cornerlist = planecorners(ii,:);
tp@0: 	cornerlist = cornerlist( find(cornerlist~= 0) );
tp@0: 	cornercoord = corners(cornerlist,:);
tp@0: 	minvals(ii,:) = min(cornercoord);
tp@0: 	maxvals(ii,:) = max(cornercoord);
tp@0: end
tp@0: 
tp@0: minvals = minvals - geomacc;
tp@0: maxvals = maxvals + geomacc;
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Check the geometry a bit
tp@0: 
tp@0: if isempty(checkgeom)
tp@0: 	checkgeom = 0;
tp@0: else
tp@0: 	if setstr(checkgeom(1)) == 'c'
tp@0: 		checkgeom = 1;
tp@0: 	else
tp@0: 		checkgeom = 0;
tp@0: 	end
tp@0: end
tp@0: 
tp@0: if checkgeom
tp@0: 
tp@0: 	badproblem = 0;
tp@0: 	disp(' ')
tp@0: 	disp('   Checking if there are corners with identical coordinates')
tp@0: 	for ii = 1:ncorners-1
tp@0: 		othercorners = [ii+1:ncorners];
tp@0: 		iv = find((corners(othercorners,1)==corners(ii,1)) & (corners(othercorners,2)==corners(ii,2)) & (corners(othercorners,3)==corners(ii,3)));
tp@0: 		if ~isempty(iv)
tp@0: 			disp(['      WARNING: corner ',int2str(ii),' and ',int2str(othercorners(iv(1))),' have the same coordinates'])
tp@0: 			disp(['      This should be fixed in the CAD-file or in the CATT-model'])
tp@0: 			badproblem = 1;
tp@0: 		end
tp@0: 	end
tp@0: 
tp@0: 	disp('   Checking if there are unused corners')
tp@0: 	iv = find(planesatcorners(:,1)==0);
tp@0: 	if ~isempty(iv)
tp@0: 		disp('      WARNING: The following corners in the CAD-file are not used:')
tp@0: 		printvec = int2str(cornernumbers(iv(1)));
tp@0: 		for iiprint = 2:length(iv)
tp@0: 			printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0: 		end
tp@0: 		disp(['            ',printvec])
tp@0: 		disp(['      This is not a big problem'])
tp@0: 		planesatcorners(iv,:) = [];
tp@0: 	end
tp@0: 
tp@0: 	disp('   Checking if any corners seem redundant')
tp@0: 	iv = find(planesatcorners(:,2)==0);
tp@0: 	if ~isempty(iv)
tp@0: 		disp('      WARNING: The following corners in the CAD-file belong to only one plane:')
tp@0: 		printvec = int2str(cornernumbers(iv(1)));
tp@0: 		for iiprint = 2:length(iv)
tp@0: 			printvec = [printvec,' ',int2str(cornernumbers(iv(iiprint)))];
tp@0: 		end
tp@0: 		disp(['            ',printvec])
tp@0: 		disp(['      This should be fixed in the CAD-file or in the CATT-model'])
tp@0: 		badproblem = 1;
tp@0: 	end
tp@0: 
tp@0: 	if badproblem == 1
tp@0: 		disp(' ');
tp@0: 		error(['ERROR: Problems with the geometry defined in the CAD-file: ',CADfile])
tp@0: 	end
tp@0: 
tp@0: end
tp@0: 
tp@0: %---------------------------------------------------------------
tp@0: % Save the relevant variables in the output file
tp@0: 
tp@0: % NB! Sparse matrices can not get a non-double format
tp@0: 
tp@0: if ncorners < 256
tp@0:     planecorners = uint8(planecorners);
tp@0: elseif ncorners < 65536
tp@0:     planecorners = uint16(planecorners);    
tp@0: end   
tp@0: 
tp@0: if max(ncornersperplanevec) <= 255
tp@0:     ncornersperplanevec = uint8(ncornersperplanevec);
tp@0:     planehasindents = uint8(planehasindents);
tp@0: else
tp@0:     ncornersperplanevec = uint16(ncornersperplanevec);
tp@0:     planehasindents = uint16(planehasindents);
tp@0: end
tp@0: 
tp@0: Varlist = [' corners cornernumbers planecorners planenames planeabstypes '];
tp@0: Varlist = [Varlist,' planenumbers cornercrossref cornercrossrefzero Snumbers Snames Sdirectivitynames Sdelays'];
tp@0: Varlist = [Varlist,' Scoords Sdirections Srotations Slevels Rnumbers Rcoords CATTversionnumber'];
tp@0: Varlist = [Varlist,' planeeqs planenvecs ncornersperplanevec planesatcorners nplanespercorners'];
tp@0: Varlist = [Varlist,' minvals maxvals planehasindents indentingcorners'];
tp@0: 
tp@0: planenames = sparse(planenames+1-1);
tp@0: planeabstypes = sparse(planeabstypes+1-1);
tp@0: Sdirectivitynames = sparse(Sdirectivitynames+1-1);
tp@0: 
tp@0: eval(['save ',outputfile,Varlist])