diff options
Diffstat (limited to 'CSXCAD/matlab')
66 files changed, 4447 insertions, 0 deletions
diff --git a/CSXCAD/matlab/AddBox.m b/CSXCAD/matlab/AddBox.m new file mode 100644 index 0000000..d07274c --- /dev/null +++ b/CSXCAD/matlab/AddBox.m @@ -0,0 +1,42 @@ +function CSX = AddBox(CSX, propName, prio, start, stop, varargin) +% function CSX = AddBox(CSX, propName, prio, start, stop, varargin) +% +% Add a box to CSX and assign to a property with name <propName>. +% +% start: box start coordinates +% stop : box stop coordinates +% prio : primitive priority +% +% optional: +% Transformation: perform a transformation on a primitive by adding +% e.g.: 'Transform', {'Scale','1,1,2','Rotate_X',pi/4,'Translate','0,0,100'} +% Note: This will only affect the 3D material/metal discretisation +% +% example: +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddBox(CSX,'metal',10,[0 0 0],[100 100 200]); %assign box +% +% with transformation: +% CSX = AddBox(CSX,'metal',10,[0 0 0],[100 100 200], ... +% 'Transform', {Rotate_Z, pi/4}); +% +% See also AddCylinder, AddCylindricalShell, AddSphere, AddSphericalShell, +% AddCurve, AddWire, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +box.ATTRIBUTE.Priority = prio; + +box.P1.ATTRIBUTE.X=start(1); +box.P1.ATTRIBUTE.Y=start(2); +box.P1.ATTRIBUTE.Z=start(3); + +box.P2.ATTRIBUTE.X=stop(1); +box.P2.ATTRIBUTE.Y=stop(2); +box.P2.ATTRIBUTE.Z=stop(3); + +box = AddPrimitiveArgs(box,varargin{:}); + +CSX = Add2Property(CSX,propName, box, 'Box');
\ No newline at end of file diff --git a/CSXCAD/matlab/AddConductingSheet.m b/CSXCAD/matlab/AddConductingSheet.m new file mode 100644 index 0000000..f3e79b4 --- /dev/null +++ b/CSXCAD/matlab/AddConductingSheet.m @@ -0,0 +1,47 @@ +function CSX = AddConductingSheet(CSX, name, conductivity, thickness) +%function CSX = AddConductingSheet(CSX, name, conductivity, thickness) +% +% Add a conducting sheet property to CSX with the given name. +% Remember to add at least one or more 2D!! geometrical primitives to this +% property. +% +% Hint: +% Set the thickness to 0 to fall back to a perfect metal (AddMetal) +% +% example: +% % create the conducting material peroperty, e.g. 40um thick copper +% CSX = AddConductingSheet(CSX,'copper',56e6,40e-6); +% % assign box the 2D box --> 40um thick sheet +% CSX = AddBox(CSX,'copper',10,[0 -50 200],[1000 50 200]); +% +% See also AddMaterial, AddMetal, AddExcitation, AddBox +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig 2012 + + +% fall back to ideal pec for t==0 or c==0 +if ((thickness==0) || (conductivity==0)) + CSX = AddMetal(CSX,name); + return; +end + +if ((conductivity<0) || (thickness<0)) + error('CSXCAD:AddConductingSheet','a negative conductivity or thickness is invalid'); +end + +if (conductivity<1e6) + warning('CSXCAD:AddConductingSheet','a conductivity below 1MA/Vm is not recommended'); +end + +if (thickness>500e-6) + warning('CSXCAD:AddConductingSheet','a thickness greater than 500um is not recommended'); +end + + +if (thickness<1e-6) + warning('CSXCAD:AddConductingSheet','a thickness lower than 1um is not recommended'); +end + +CSX = AddProperty(CSX, 'ConductingSheet', name,'Conductivity',conductivity,'Thickness',thickness); diff --git a/CSXCAD/matlab/AddCurve.m b/CSXCAD/matlab/AddCurve.m new file mode 100644 index 0000000..2e4e5b2 --- /dev/null +++ b/CSXCAD/matlab/AddCurve.m @@ -0,0 +1,44 @@ +function CSX = AddCurve(CSX, propName, prio, points, varargin) +% function CSX = AddCurve(CSX, propName, prio, points, varargin) +% +% Add a curve to CSX and assign to a property with name <propName>. +% +% Warning: This is a 1D object, not all properties may be compatible with a +% 1D object, e.g. a material property. +% +% points: curve coordinates array +% prio : primitive priority +% +% example: +% %first point +% points(1,1) = 0; +% points(2,1) = 5; +% points(3,1) = 10; +% %second point +% points(1,2) = 0; +% points(2,2) = 10; +% points(3,2) = 10; +% %third point ... +% % create a thin metal wire... +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddCurve(CSX,'metal',10, points); +% +% See also AddBox, AddCylindricalShell, AddCylinder, AddSphere, +% AddSphericalShell, AddWire, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +curve.ATTRIBUTE.Priority = prio; + +curve.Vertex={}; +for s=1:size(points,2) + curve.Vertex{end+1}.ATTRIBUTE.X = points(1,s); + curve.Vertex{end}.ATTRIBUTE.Y = points(2,s); + curve.Vertex{end}.ATTRIBUTE.Z = points(3,s); +end + +curve = AddPrimitiveArgs(curve,varargin{:}); + +CSX = Add2Property(CSX,propName, curve, 'Curve');
\ No newline at end of file diff --git a/CSXCAD/matlab/AddCylinder.m b/CSXCAD/matlab/AddCylinder.m new file mode 100644 index 0000000..1998325 --- /dev/null +++ b/CSXCAD/matlab/AddCylinder.m @@ -0,0 +1,35 @@ +function CSX = AddCylinder(CSX, propName, prio, start, stop, rad, varargin) +% function CSX = AddCylinder(CSX, propName, prio, start, stop, rad, varargin) +% +% Add a cylinder to CSX and assign to a property with name <propName>. +% +% start: cylinder axis start coordinates +% stop : cylinder axis box stop coordinates +% rad : cylinder radius +% prio : primitive priority +% +% example: +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddCylinder(CSX,'metal',10,[0 0 0],[0 0 200],50); +% +% See also AddBox, AddCylindricalShell, AddSphere, AddSphericalShell, +% AddCurve, AddWire, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +cylinder.ATTRIBUTE.Priority = prio; +cylinder.ATTRIBUTE.Radius = rad; + +cylinder.P1.ATTRIBUTE.X=start(1); +cylinder.P1.ATTRIBUTE.Y=start(2); +cylinder.P1.ATTRIBUTE.Z=start(3); + +cylinder.P2.ATTRIBUTE.X=stop(1); +cylinder.P2.ATTRIBUTE.Y=stop(2); +cylinder.P2.ATTRIBUTE.Z=stop(3); + +cylinder = AddPrimitiveArgs(cylinder,varargin{:}); + +CSX = Add2Property(CSX,propName, cylinder, 'Cylinder');
\ No newline at end of file diff --git a/CSXCAD/matlab/AddCylindricalShell.m b/CSXCAD/matlab/AddCylindricalShell.m new file mode 100644 index 0000000..53ff5f4 --- /dev/null +++ b/CSXCAD/matlab/AddCylindricalShell.m @@ -0,0 +1,41 @@ +function CSX = AddCylindricalShell(CSX, propName, prio, start, stop, rad, shell_width, varargin) +%function CSX = AddCylindricalShell(CSX, propName, prio, start, stop, rad, shell_width, varargin) +% +% Add a cylinder shell to CSX and assign to a property with name <propName>. +% +% start: cylinder axis start coordinates +% stop : cylinder axis box stop coordinates +% rad : cylinder radius +% shell_width: cylinder shell width +% prio : primitive priority +% +% Note: +% the inner radius of this shell is rad-shell_width/2 +% the outer radius of this shell is rad+shell_width/2 +% +% example: +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddCylindricalShell(CSX,'metal',10,[0 0 0],[0 0 200],50,10); +% +% See also AddBox, AddCylinder, AddSphere, AddSphericalShell, +% AddCurve, AddWire, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +cylinder.ATTRIBUTE.Priority = prio; +cylinder.ATTRIBUTE.Radius = rad; +cylinder.ATTRIBUTE.ShellWidth = shell_width; + +cylinder.P1.ATTRIBUTE.X=start(1); +cylinder.P1.ATTRIBUTE.Y=start(2); +cylinder.P1.ATTRIBUTE.Z=start(3); + +cylinder.P2.ATTRIBUTE.X=stop(1); +cylinder.P2.ATTRIBUTE.Y=stop(2); +cylinder.P2.ATTRIBUTE.Z=stop(3); + +cylinder = AddPrimitiveArgs(cylinder,varargin{:}); + +CSX = Add2Property(CSX,propName, cylinder, 'CylindricalShell'); diff --git a/CSXCAD/matlab/AddDebyeMaterial.m b/CSXCAD/matlab/AddDebyeMaterial.m new file mode 100644 index 0000000..941ffe9 --- /dev/null +++ b/CSXCAD/matlab/AddDebyeMaterial.m @@ -0,0 +1,29 @@ +function CSX = AddDebyeMaterial(CSX, name, varargin) +% function CSX = AddDebyeMaterial(CSX, name, varargin) +% +% Add a Debye type dispersive material model. +% +% The Debye type frequency dependent material: +% eps_r(w) = eps_r + sum_p ( eps_r_delta,p / (1+jw*t_relex,p) ) - j*kappa/w +% +% with +% eps_r_delta,p: the delta electric relative permitivity +% t_relex,p: the electric relaxation time +% +% Use SetMaterialProperty to define the material constants: +% 'EpsilonDelta_p': p-th eps_r_delta,p +% 'EpsilonRelaxTime_p': p-th t_relex,p +% +% example: +% CSX = AddDebyeMaterial(CSX,'debye'); +% CSX = SetMaterialProperty(CSX,'debye','Epsilon',5,'EpsilonDelta_1',0.1,'EpsilonRelaxTime_1',1e-9); +% [..] +% CSX = AddBox(CSX,'debye', 10 ,start,stop); +% +% See also AddBox, AddMaterial, SetMaterialProperty, AddLorentzMaterial +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (2013) + +CSX = AddProperty(CSX, 'DebyeMaterial', name, varargin{:}); diff --git a/CSXCAD/matlab/AddDiscMaterial.m b/CSXCAD/matlab/AddDiscMaterial.m new file mode 100644 index 0000000..86f5213 --- /dev/null +++ b/CSXCAD/matlab/AddDiscMaterial.m @@ -0,0 +1,47 @@ +function CSX = AddDiscMaterial(CSX, name, varargin) +% function CSX = AddDiscMaterial(CSX, name, varargin) +% +% Add a discretized material model property to CSX with the given name. +% Discretized model has to be stored in an hdf5 file. +% Use Transform option to perfom some transformation actions. +% +% variable arguments (key/value) +% 'File' (mandatory): define the filename of the discrete material +% 'Scale': scale the discrete material +% e.g. to your drawing units: 'Scale', 1/unit +% 'Transform': Apply a transformation, see AddBox for more infos +% 'UseDBBackground': set to 0, to use the properties background material +% instead of the database material with index 0 (default) +% +% examples: +% %add human body model +% CSX = AddDiscMaterial(CSX, 'Ella', 'Scale', 1/unit, ... +% 'Transform', {'Rotate_Z',pi/2,'Translate','300,300,500'}, ... +% 'File', 'model_file_name.h5' ); +% start = [mesh.x(1) mesh.y(1) mesh.z(1)]; +% stop = [mesh.x(end) mesh.y(end) mesh.z(end)]; +% CSX = AddBox(CSX,'Ella', 0 ,start,stop); +% +% See also AddBox, AddMetal, AddExcitation, AddProbe, AddDump +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +transform = []; + +for n=1:2:(nargin-2) + if (strcmp(varargin{n},'Transform')) + transform = varargin([n:n+1]); + varargin([n:n+1]) = []; + break + end +end + +CSX = AddProperty(CSX, 'DiscMaterial', name, varargin{:}); + +if ~isempty(transform) + for n=1:2:numel(transform{2}) + CSX.Properties.DiscMaterial{end}.Transformation.(transform{2}{n}).ATTRIBUTE.Argument=transform{2}{n+1}; + end +end diff --git a/CSXCAD/matlab/AddDump.m b/CSXCAD/matlab/AddDump.m new file mode 100644 index 0000000..f454992 --- /dev/null +++ b/CSXCAD/matlab/AddDump.m @@ -0,0 +1,87 @@ +function CSX = AddDump(CSX, name, varargin) +% function CSX = AddDump(CSX, name, varargin) +% +% Add a dump property to CSX with the given name. +% +% possible arguments for useage with openEMS: +% DumpType: 0 for E-field time-domain dump (default) +% 1 for H-field time-domain dump +% 2 for electric current time-domain dump +% 3 for total current density (rot(H)) time-domain dump +% +% 10 for E-field frequency-domain dump +% 11 for H-field frequency-domain dump +% 12 for electric current frequency-domain dump +% 13 for total current density (rot(H)) frequency-domain dump +% +% 20 local SAR frequency-domain dump +% 21 1g averaging SAR frequency-domain dump +% 22 10g averaging SAR frequency-domain dump +% +% 29 raw data needed for SAR calculations (electric field FD, +% cell volume, conductivity and density) +% +% Frequency: specify a frequency vector (required for dump types >=10) +% +% DumpMode: 0 no-interpolation +% 1 node-interpolation (default, see warning below) +% 2 cell-interpolation (see warning below) +% +% FileType: 0 vtk-file dump (default) +% 1 hdf5-file dump +% +% SubSampling: field domain sub-sampling, e.g. '2,2,4' +% OptResolution: field domain dump resolution, e.g. '10' or '10,20,5' +% +% MultiGridLevel: Request to dump from a multigrid level (default is 0) +% Note: This only takes effect if the method supports and +% uses multiple grids! +% +% StartTime/StopTime: Define a start and/or stop time (in seconds) +% for this dump to be active. +% +% Warning: +% FDTD Interpolation abnormalities: +% - no-interpolation: fields are located in the mesh by the +% Yee-scheme, the mesh only specifies E- or H-Yee-nodes +% --> use node- or cell-interpolation or be aware of the offset +% - E-field dump & node-interpolation: normal electric fields on +% boundaries will have false amplitude due to forward/backward +% interpolation in case of (strong) changes in material +% permittivity or on metal surfaces +% --> use no- or cell-interpolation +% - H-field dump & cell-interpolation: normal magnetic fields on +% boundaries will have false amplitude due to forward/backward +% interpolation in case of (strong) changes in material permeability +% --> use no- or node-interpolation +% +% e.g. AddDump(CSX,'Et'); +% CSX = AddBox(CSX,'Et',10,[0 0 0],[100 100 200]); %assign box +% +% or AddDump(CSX,'Ef',DumpType, 10, 'Frequency',[1e9 2e9]); +% CSX = AddBox(CSX,'Ef',10,[0 0 0],[100 100 200]); %assign box +% +% or AddDump(CSX,'Ht','SubSampling','2,2,4','DumpType',1); +% CSX = AddBox(CSX,'Ht',10,[0 0 0],[100 100 200]); %assign box +% +% See also AddMaterial, AddExcitation, AddProbe, AddMetal, AddBox +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +[CSX pos] = AddProperty(CSX, 'DumpBox', name); + +% make Node-Interpolation the default +CSX.Properties.DumpBox{pos}.ATTRIBUTE.DumpMode = 1; + +for n=1:numel(varargin)/2 + if ~ischar(varargin{2*n-1}) + error(['CSXCAD::AddDump: not an attribute: ' varargin{2*n-1}]); + end + if strcmp(varargin{2*n-1},'Frequency') + CSX.Properties.DumpBox{pos}.FD_Samples=varargin{2*n}; + else + CSX.Properties.DumpBox{pos}.ATTRIBUTE.(varargin{2*n-1})=varargin{2*n}; + end +end diff --git a/CSXCAD/matlab/AddExcitation.m b/CSXCAD/matlab/AddExcitation.m new file mode 100644 index 0000000..45a1a2c --- /dev/null +++ b/CSXCAD/matlab/AddExcitation.m @@ -0,0 +1,32 @@ +function CSX = AddExcitation(CSX, name, type, excite, varargin) +% function CSX = AddExcitation(CSX, name, type, excite, varargin) +% +% Creates an E-field or H-field excitation. +% +% CSX: CSX-struct created by InitCSX +% name: property name for the excitation +% type: 0=E-field soft excitation 1=E-field hard excitation +% 2=H-field soft excitation 3=H-field hard excitation +% 10=plane wave excitation +% +% excite: e.g. [2 0 0] for excitation of 2 V/m in x-direction +% +% additional options for openEMS: +% 'Delay' : setup an excitation time delay in seconds +% 'PropDir': direction of plane wave propagation (plane wave excite only) +% +% example: +% CSX = AddExcitation( CSX, 'infDipole', 1, [1 0 0] ); +% start = [-dipole_length/2, 0, 0]; +% stop = [+dipole_length/2, 0, 0]; +% CSX = AddBox( CSX, 'infDipole', 1, start, stop ); +% +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig +% +% See also SetExcitationWeight, AddMetal, AddExcitation, AddProbe, +% AddDump, AddBox + +CSX = AddProperty(CSX, 'Excitation', name, 'Type', type, 'Excite', excite, varargin{:}); diff --git a/CSXCAD/matlab/AddLinPoly.m b/CSXCAD/matlab/AddLinPoly.m new file mode 100644 index 0000000..b93a4ac --- /dev/null +++ b/CSXCAD/matlab/AddLinPoly.m @@ -0,0 +1,41 @@ +function CSX = AddLinPoly( CSX, materialname, prio, normDir, elevation, points, Length, varargin) +% CSX = AddLinPoly( CSX, materialname, prio, normDir, elevation, points, Length, varargin) +% +% CSX: CSX-object created by InitCSX() +% materialname: created by AddMetal() or AddMaterial() +% prio: priority +% normDir: normal direction of the polygon, +% e.g. 'x', 'y' or 'z', or numeric 0..2 +% elevation: position of the polygon plane +% points: two-dimensional coordinates +% length: linear extrution in normal direction, starting at elevation +% +% Warning: Polygon has to be defined using Cartesian Coords +% for use with cylindrical mesh, set 'CoordSystem',0 +% +% example: +% p(1,1) = 0; % x-coord point 1 +% p(2,1) = 0; % y-coord point 1 +% p(1,2) = 10; % x-coord point 2 +% p(2,2) = 20; % y-coord point 2 +% % normal direction: z (2) +% CSX = AddLinPoly( CSX, 'PEC', 1, 2, 254, p , 10, 'CoordSystem',0) +% +% 2011, Thorsten Liebig <thorsten.liebig@gmx.de> +% +% See also InitCSX AddMetal AddMaterial AddPolygon AddRotPoly + +polygon.ATTRIBUTE.Priority = prio; +polygon.ATTRIBUTE.Elevation = elevation; +polygon.ATTRIBUTE.Length = Length; +polygon.ATTRIBUTE.NormDir = DirChar2Int(normDir); + +polygon.Vertex = {}; +for s=1:size(points,2) + polygon.Vertex{end+1}.ATTRIBUTE.X1 = points(1,s); + polygon.Vertex{end}.ATTRIBUTE.X2 = points(2,s); +end + +polygon = AddPrimitiveArgs(polygon,varargin{:}); + +CSX = Add2Property( CSX, materialname, polygon, 'LinPoly'); diff --git a/CSXCAD/matlab/AddLorentzMaterial.m b/CSXCAD/matlab/AddLorentzMaterial.m new file mode 100644 index 0000000..b0e9629 --- /dev/null +++ b/CSXCAD/matlab/AddLorentzMaterial.m @@ -0,0 +1,60 @@ +function CSX = AddLorentzMaterial(CSX, name, varargin) +% function CSX = AddLorentzMaterial(CSX, name, varargin) +% +% Add a Drude/Lorentz type dispersive material model. +% Note: openEMS currently only supports a drude material type. +% +% The drude type frequency dependent material: +% eps_r(f) = eps_r* ( 1 - f_eps_plasma^2/(f*(f-j/t_eps_r)) ) +% mue_r(f) = mue_r* ( 1 - f_mue_plasma^2/(f*(f-j/t_mue_r)) ) +% with +% f_eps_plasma: the respective electric angular plasma frequency +% f_mue_plasma: the respective magnetic angular plasma frequency +% t_eps_r: the respective electric relaxation time +% t_mue_r: the respective magnetic relaxation time +% +% Use SetMaterialProperty to define the material constants: +% 'EpsilonPlasmaFrequency': electric plasma frequency (f_eps_plasma) +% 'MuePlasmaFrequency': magnetic plasma frequency (f_mue_plasma) +% 'EpsilonRelaxTime': electric plasma relaxation time (losses) +% 'MueRelaxTime': magnetic plasma relaxation time (losses) +% +% Note: all properties must be positive values +% +% Higher order Drude type: +% 'EpsilonPlasmaFrequency_<n>': n-th order electric plasma frequency (f_eps_plasma) +% 'MuePlasmaFrequency_<n>': n-th order magnetic plasma frequency (f_mue_plasma) +% 'EpsilonRelaxTime_<n>': n-th order electric plasma relaxation time (losses) +% 'MueRelaxTime_<n>': n-th order magnetic plasma relaxation time (losses) +% +% The Lorentz type frequency dependent material: +% eps_r(f) = eps_r* ( 1 - f_eps_plasma^2/(f^2-f_eps_Lor_Pole^2-jf^2*/t_eps_r)) ) +% mue_r(f) = mue_r* ( 1 - f_mue_plasma^2/(f^2-f_mue_Lor_Pole^2-jf^2*/t_mue_r)) ) +% with the additional parameter (see above) +% f_eps_Lor_Pole: the respective electric angular lorentz pole frequency +% f_mue_Lor_Pole: the respective magnetic angular lorentz pole frequency +% +% Use SetMaterialProperty to define the material constants: +% 'EpsilonLorPoleFrequency': electric lorentz pole frequency (f_eps_Lor_Pole) +% 'MueLorPoleFrequency': magnetic lorentz pole frequency (f_mue_Lor_Pole) +% +% Note: all properties must be positive values +% +% Higher order Drude type: +% 'EpsilonLorPoleFrequency_<n>': n-th order electric lorentz pole frequency (f_eps_plasma) +% 'MueLorPoleFrequency_<n>': n-th order magnetic lorentz pole frequency (f_mue_plasma) +% +% example: +% CSX = AddLorentzMaterial(CSX,'drude'); +% CSX = SetMaterialProperty(CSX,'drude','Epsilon',5,'EpsilonPlasmaFrequency',5e9,'EpsilonRelaxTime',1e-9); +% CSX = SetMaterialProperty(CSX,'drude','Mue',5,'MuePlasmaFrequency',5e9,'MueRelaxTime',1e-9); +% [..] +% CSX = AddBox(CSX,'drude', 10 ,start,stop); +% +% See also AddBox, AddMaterial, SetMaterialProperty +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +CSX = AddProperty(CSX, 'LorentzMaterial', name, varargin{:}); diff --git a/CSXCAD/matlab/AddLumpedElement.m b/CSXCAD/matlab/AddLumpedElement.m new file mode 100644 index 0000000..cefa4c7 --- /dev/null +++ b/CSXCAD/matlab/AddLumpedElement.m @@ -0,0 +1,29 @@ +function CSX = AddLumpedElement(CSX, name, direction, varargin) +% function CSX = AddLumpedElement(CSX, name, direction, varargin) +% +% Add a lumped element property to CSX with the given name. +% Remember to add at least one or more box primitive to any +% property. +% +% arguments: +% direction: 0,1,2 or 'x','y','z' for the orientation of the lumped element +% +% optional arguments: +% 'R', 'C', 'L': definition of the lumped element properties +% 'Caps': 0 or 1 to (de)activate lumped element caps (1=default) +% If Caps are enable, a small PEC plate is added to each +% end of the lumped element to ensure contact to e.g +% a microstrip line +% +% examples: +% lumped element capacitor in y-direction with 1pF +% CSX = AddLumpedElement( CSX, 'Capacitor', 1, 'Caps', 1, 'C', 1e-12 ); +% CSX = AddBox( CSX, 'Capacitor', 0, [0 0 0], [10 10 10] ); +% +% See also AddMaterial, AddMetal, AddExcitation, AddProbe, AddDump, AddBox +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +CSX = AddProperty(CSX, 'LumpedElement', name, 'Direction', DirChar2Int(direction), varargin{:}); diff --git a/CSXCAD/matlab/AddMaterial.m b/CSXCAD/matlab/AddMaterial.m new file mode 100644 index 0000000..d21a8c9 --- /dev/null +++ b/CSXCAD/matlab/AddMaterial.m @@ -0,0 +1,26 @@ +function CSX = AddMaterial(CSX, name, varargin) +% function CSX = AddMaterial(CSX, name, varargin) +% +% Add a material property to CSX with the given name. +% Remember to add at least one or more geometrical primitives to any +% property. +% +% Use SetMaterialProperty to define the material constants: +% 'Epsilon': relative electric permitivity +% 'Mue': relative magnetic permeability +% 'Kappa': electric conductivity +% 'Sigma': magnetc conductivity (non-physical property) +% +% examples: +% CSX = AddMaterial( CSX, 'RO3010' ); +% CSX = SetMaterialProperty( CSX, 'RO3010', 'Epsilon', 10.2, 'Mue', 1 ); +% CSX = AddBox( CSX, 'RO3010', 0, [0 0 0], [100 1000 1000] ); +% +% See also SetMaterialProperty, SetMaterialWeight, AddMetal, AddExcitation, +% AddProbe, AddDump, AddBox +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +CSX = AddProperty(CSX, 'Material', name, varargin{:}); diff --git a/CSXCAD/matlab/AddMetal.m b/CSXCAD/matlab/AddMetal.m new file mode 100644 index 0000000..98ce9d1 --- /dev/null +++ b/CSXCAD/matlab/AddMetal.m @@ -0,0 +1,26 @@ +function CSX = AddMetal(CSX, name) +%function CSX = AddMetal(CSX, name) +% +% Add a metal property (PEC) to CSX with the given name. +% Remember to add at least one or more geometrical primitives to any +% property. +% +% example: +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddBox(CSX,'metal',10,[0 0 0],[100 100 200]); %assign box +% +% See also AddMaterial, AddExcitation, AddProbe, AddDump, AddBox +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +% check if this property already exists +[type_found pos] = FindProperty(CSX, name); + +% since it has no varargs, accept already existing metal with this name +if ((pos>0) && strcmp(type_found,'Metal')) + return +end + +CSX = AddProperty(CSX, 'Metal', name); diff --git a/CSXCAD/matlab/AddPlaneWaveExcite.m b/CSXCAD/matlab/AddPlaneWaveExcite.m new file mode 100644 index 0000000..aa59c02 --- /dev/null +++ b/CSXCAD/matlab/AddPlaneWaveExcite.m @@ -0,0 +1,38 @@ +function CSX = AddPlaneWaveExcite(CSX, name, k_dir, E_dir, f0, varargin) +% function CSX = AddPlaneWaveExcite(CSX, name, k_dir, E_dir, <f0, varargin>) +% +% Creates a plane wave excitation in the sense of a total-field/scattered +% field approach. +% +% Note: A plane wave excitation must not intersect with any kind of +% material. This exctiation type can only be applies in air/vacuum and +% completely surrounding a structure!!! +% +% Note: Only a single Box can be applied to this property!! +% +% Arguments +% CSX: CSX-struct created by InitCSX +% name: property name for the excitation +% k_dir: unit vector of wave progation direction +% E_dir: electric field polarisation vector (must be orthogonal to k_dir) +% f0: frequency for numerical phase velocity compensation (optional) +% +% example: +% inc_angle = 0 /180*pi; %incident angle on the x-axis +% k_dir = [cos(inc_angle) sin(inc_angle) 0]; % plane wave direction +% E_dir = [0 0 1]; % plane wave polarization --> E_z +% f0 = 500e6; % frequency for numerical phase velocity compensation +% +% CSX = AddPlaneWaveExcite(CSX, 'plane_wave', k_dir, E_dir, f0); +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig +% +% See also AddExcitation, AddBox + +if (nargin<5) + f0 = 0; +end + +CSX = AddExcitation(CSX, name, 10, E_dir, 'PropDir', k_dir, 'Frequency', f0, varargin{:}); diff --git a/CSXCAD/matlab/AddPoint.m b/CSXCAD/matlab/AddPoint.m new file mode 100644 index 0000000..15af7d4 --- /dev/null +++ b/CSXCAD/matlab/AddPoint.m @@ -0,0 +1,14 @@ +function CSX = AddPoint(CSX, propName, prio, pos, varargin) +% CSX = AddPoint(CSX, propName, prio, pos, varargin) +% +% CSXCAD matlab interface + +point.ATTRIBUTE.Priority = prio; + +point.ATTRIBUTE.X=pos(1); +point.ATTRIBUTE.Y=pos(2); +point.ATTRIBUTE.Z=pos(3); + +point = AddPrimitiveArgs(point,varargin{:}); + +CSX = Add2Property(CSX,propName, point, 'Point'); diff --git a/CSXCAD/matlab/AddPolygon.m b/CSXCAD/matlab/AddPolygon.m new file mode 100644 index 0000000..4590248 --- /dev/null +++ b/CSXCAD/matlab/AddPolygon.m @@ -0,0 +1,41 @@ +function CSX = AddPolygon( CSX, materialname, prio, normDir, elevation, points, varargin) +% CSX = AddPolygon( CSX, materialname, prio, normDir, elevation, points, varargin) +% +% CSX: CSX-object created by InitCSX() +% materialname: created by AddMetal() or AddMaterial() +% prio: priority +% normDir: normal direction of the polygon, +% e.g. 'x', 'y' or 'z', or numeric 0..2 +% elevation: position of the polygon plane +% points: two-dimensional coordinates +% +% Warning: Polygon has to be defined using Cartesian Coords +% for use with cylindrical mesh, set 'CoordSystem',0 +% +% example: +% p(1,1) = 0; % x-coord point 1 +% p(2,1) = 0; % y-coord point 1 +% p(1,2) = 10; % x-coord point 2 +% p(2,2) = 20; % y-coord point 2 +% % normal direction: z (2) +% CSX = AddPolygon( CSX, 'PEC', 1, 'z', 254, p, 'CoordSystem',0) +% +% +% (c) 2011 Thorsten Liebig <thorsten.liebig@gmx.de> +% (c) 2010 Sebastian Held <sebastian.held@gmx.de> +% +% See also InitCSX AddMetal AddMaterial AddLinPoly AddRotPoly + +polygon.ATTRIBUTE.Priority = prio; +polygon.ATTRIBUTE.Elevation = elevation; +polygon.ATTRIBUTE.NormDir = DirChar2Int(normDir); + +polygon.Vertex = {}; +for s=1:size(points,2) + polygon.Vertex{end+1}.ATTRIBUTE.X1 = points(1,s); + polygon.Vertex{end}.ATTRIBUTE.X2 = points(2,s); +end + +polygon = AddPrimitiveArgs(polygon,varargin{:}); + +CSX = Add2Property( CSX, materialname, polygon, 'Polygon'); diff --git a/CSXCAD/matlab/AddPolyhedron.m b/CSXCAD/matlab/AddPolyhedron.m new file mode 100644 index 0000000..7da8a86 --- /dev/null +++ b/CSXCAD/matlab/AddPolyhedron.m @@ -0,0 +1,51 @@ +function CSX = AddPolyhedron(CSX, propName, prio, vertices, faces, varargin) +% CSX = AddPolyhedron(CSX, propName, prio, vertices, faces, varargin) +% +% Add a polyhedron to CSX and assign to a property with name <propName>. +% +% prio: primitive priority +% vertices: cell array of all vertices +% faces: cell array of all faces +% +% Note: - The polyhedron must be a closed surface for 3D discretisation +% - All faces must contain the vertices in a right-handed order with +% the normal direction for each face pointing out of the solid +% +% optional: +% Transformation: perform a transformation on a primitive by adding +% e.g.: 'Transform', {'Scale','1,1,2','Rotate_X',pi/4,'Translate','0,0,100'} +% Note: This will only affect the 3D material/metal discretisation +% +% example: +% % example tetrahedron +% vertices{1}=[0 0 0]; +% vertices{2}=[1 0 0]; +% vertices{3}=[0 1 0]; +% vertices{4}=[0 0 1]; +% faces{1}=[0 2 1]; +% faces{2}=[0 1 3]; +% faces{3}=[0 3 2]; +% faces{4}=[1 2 3]; +% CSX = AddMetal( CSX, 'metal' ); +% CSX = AddPolyhedron(CSX, 'metal', 0, vertices, faces); +% +% +% See also AddBox, AddCylinder, AddCylindricalShell, AddSphere, AddSphericalShell, +% AddCurve, AddWire, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +polyhedron.ATTRIBUTE.Priority = prio; + +for n=1:numel(vertices) + polyhedron.Vertex{n}=vector2str(vertices{n}); +end +for n=1:numel(faces) + polyhedron.Face{n}=vector2str(faces{n}); +end + +polyhedron = AddPrimitiveArgs(polyhedron,varargin{:}); + +CSX = Add2Property(CSX,propName, polyhedron, 'Polyhedron'); diff --git a/CSXCAD/matlab/AddProbe.m b/CSXCAD/matlab/AddProbe.m new file mode 100644 index 0000000..2a6692e --- /dev/null +++ b/CSXCAD/matlab/AddProbe.m @@ -0,0 +1,79 @@ +function CSX = AddProbe(CSX, name, type, varargin) +% function CSX = AddProbe(CSX, name, type, varargin) +% +% Add a probe property to CSX with the given name. +% Remember to add a geometrical primitive to any property. +% +% name: name of the property and probe file +% +% type: 0 for voltage probing +% 1 for current probing +% 2 for E-field probing +% 3 for H-field probing +% +% 10 for waveguide voltage mode matching +% 11 for waveguide current mode matching +% +% all following parameter are optional key/value parameter: +% +% weight: weighting factor (default is 1) +% frequency: dump in the frequency domain at the given samples (in Hz) +% ModeFunction: A mode function (used only with type 3/4) +% NormDir: necessary for current probing box with dimension~=2 +% StartTime/StopTime: Define a start and/or stop time (in seconds) +% for this probe to be active. + +% examples: +% CSX = AddProbe(CSX,'ut1',0); %voltate probe +% CSX = AddProbe(CSX,'it1',1); %current probe +% +% See also ReadUI in the openEMS matlab interface, AddDump, +% AddExcitation, AddMaterial, AddExcitation, AddProbe, AddBox +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +FD_samples = []; +ModeFunction = {}; + +if ~ischar(name) + error('CSXCAD::AddProbe: name must be a string'); +end + +prop_args = {'Type', type}; + +for n=1:2:numel(varargin) + if (strcmpi(varargin{n},'weight')==1); + prop_args{end+1} = 'Weight'; + prop_args{end+1} = varargin{n+1}; + elseif (strcmpi(varargin{n},'Frequency')==1); + FD_samples = varargin{n+1}; + elseif (strcmpi(varargin{n},'ModeFunction')==1); + ModeFunction = varargin{n+1}; + elseif (strcmpi(varargin{n},'NormDir')==1); + prop_args{end+1} = 'NormDir'; + prop_args{end+1} = varargin{n+1}; + elseif (strcmpi(varargin{n},'StartTime')==1); + prop_args{end+1} = 'StartTime'; + prop_args{end+1} = varargin{n+1}; + elseif (strcmpi(varargin{n},'StopTime')==1); + prop_args{end+1} = 'StopTime'; + prop_args{end+1} = varargin{n+1}; + else + warning('CSXCAD:AddProbe',['variable argument key: "' varargin{n+1} '" unknown']); + end +end + +[CSX pos] = AddProperty(CSX, 'ProbeBox', name, prop_args{:}); + +if (numel(FD_samples)>0) + CSX.Properties.ProbeBox{pos}.FD_Samples=FD_samples; +end + +if (numel(ModeFunction)>0) + CSX.Properties.ProbeBox{pos}.Attributes.ATTRIBUTE.ModeFunctionX = ModeFunction{1}; + CSX.Properties.ProbeBox{pos}.Attributes.ATTRIBUTE.ModeFunctionY = ModeFunction{2}; + CSX.Properties.ProbeBox{pos}.Attributes.ATTRIBUTE.ModeFunctionZ = ModeFunction{3}; +end + diff --git a/CSXCAD/matlab/AddPropAttribute.m b/CSXCAD/matlab/AddPropAttribute.m new file mode 100644 index 0000000..05408f9 --- /dev/null +++ b/CSXCAD/matlab/AddPropAttribute.m @@ -0,0 +1,24 @@ +function CSX = AddPropAttribute(CSX, name, att_name, att_value) +% CSX = AddPropAttribute(CSX, name, att_name, att_value) +% +% Add a given attribute (name and value) to the given property +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (c) 2013 + +type = GetPropertyType(CSX, name); + +pos=0; +for n=1:numel(CSX.Properties.(type)) + if strcmp(CSX.Properties.(type){n}.ATTRIBUTE.Name, name) + pos=n; + end +end + +if (pos==0) + error('CSXCAD::AddPropAttribute: property not found'); + return; +end + +CSX.Properties.(type){pos}.Attributes.ATTRIBUTE.(att_name) = att_value; diff --git a/CSXCAD/matlab/AddRotPoly.m b/CSXCAD/matlab/AddRotPoly.m new file mode 100644 index 0000000..e694e98 --- /dev/null +++ b/CSXCAD/matlab/AddRotPoly.m @@ -0,0 +1,48 @@ +function CSX = AddRotPoly( CSX, materialname, prio, normDir, points, RotAxisDir, angle, varargin) +% function CSX = AddRotPoly( CSX, materialname, prio, normDir, points, RotAxisDir, angle, varargin) +% +% CSX: CSX-object created by InitCSX() +% materialname: created by AddMetal() or AddMaterial() +% prio: priority +% normDir: normal direction of the polygon, +% e.g. 'x', 'y' or 'z', or numeric 0..2 +% points: two-dimensional coordinates +% RotAxisDir: direction of the rotational axis +% e.g. 'x', 'y' or 'z', or numeric 0..2 +% angle (optional): rotational start/stop angle, default is [0 2pi] +% +% Warning: Polygon has to be defined using Cartesian Coords +% for use with cylindrical mesh, set 'CoordSystem',0 +% +% example: +% p(1,1) = 0; % x-coord point 1 +% p(2,1) = 0; % y-coord point 1 +% p(1,2) = 10; % x-coord point 2 +% p(2,2) = 20; % y-coord point 2 +% % normal direction: z +% CSX = AddRotPoly( CSX, 'PEC', 1, 'z', p , 'y'); +% +% 2011, Thorsten Liebig <thorsten.liebig@gmx.de> +% +% See also InitCSX AddMetal AddMaterial AddPolygon + +if (nargin<7) + angle = [0 2*pi]; +end + +polygon.ATTRIBUTE.Priority = prio; +polygon.ATTRIBUTE.NormDir = DirChar2Int(normDir); +polygon.ATTRIBUTE.RotAxisDir = DirChar2Int(RotAxisDir); + +polygon.Angles.ATTRIBUTE.Start = angle(1); +polygon.Angles.ATTRIBUTE.Stop = angle(2); + +polygon.Vertex = {}; +for s=1:size(points,2) + polygon.Vertex{end+1}.ATTRIBUTE.X1 = points(1,s); + polygon.Vertex{end}.ATTRIBUTE.X2 = points(2,s); +end + +polygon = AddPrimitiveArgs(polygon,varargin{:}); + +CSX = Add2Property( CSX, materialname, polygon, 'RotPoly'); diff --git a/CSXCAD/matlab/AddSphere.m b/CSXCAD/matlab/AddSphere.m new file mode 100644 index 0000000..457142f --- /dev/null +++ b/CSXCAD/matlab/AddSphere.m @@ -0,0 +1,30 @@ +function CSX = AddSphere(CSX, propName, prio, center, rad, varargin) +% function CSX = AddSphere(CSX, propName, prio, center, rad, varargin) +% +% Add a sphere to CSX and assign to a property with name <propName>. +% +% center: sphere center coordinates +% rad : sphere radius +% prio : primitive priority +% +% example: +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddSphere(CSX,'metal',10,[0 0 0],50); +% +% See also AddBox, AddCylindricalShell, AddCylinder, AddSphericalShell, +% AddCurve, AddWire, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +sphere.ATTRIBUTE.Priority = prio; +sphere.ATTRIBUTE.Radius = rad; + +sphere.Center.ATTRIBUTE.X=center(1); +sphere.Center.ATTRIBUTE.Y=center(2); +sphere.Center.ATTRIBUTE.Z=center(3); + +sphere = AddPrimitiveArgs(sphere,varargin{:}); + +CSX = Add2Property(CSX,propName, sphere, 'Sphere');
\ No newline at end of file diff --git a/CSXCAD/matlab/AddSphericalShell.m b/CSXCAD/matlab/AddSphericalShell.m new file mode 100644 index 0000000..3e76cef --- /dev/null +++ b/CSXCAD/matlab/AddSphericalShell.m @@ -0,0 +1,36 @@ +function CSX = AddSphericalShell(CSX, propName, prio, center, rad, shell_width, varargin) +% function CSX = AddSphericalShell(CSX, propName, prio, center, rad, shell_width, varargin) +% +% Add a sphere shell to CSX and assign to a property with name <propName>. +% +% center: sphere center coordinates +% rad : sphere radius +% shell_width: sphere shell width +% prio : primitive priority +% +% Note: +% the inner radius of this shell is rad-shell_width/2 +% the outer radius of this shell is rad+shell_width/2 +% +% example: +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddSphericalShell(CSX,'metal',10,[0 0 0],50,10); +% +% See also AddBox, AddCylindricalShell, AddCylinder, AddSphere, +% AddCurve, AddWire, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +sphere.ATTRIBUTE.Priority = prio; +sphere.ATTRIBUTE.Radius = rad; +sphere.ATTRIBUTE.ShellWidth = shell_width; + +sphere.Center.ATTRIBUTE.X=center(1); +sphere.Center.ATTRIBUTE.Y=center(2); +sphere.Center.ATTRIBUTE.Z=center(3); + +sphere = AddPrimitiveArgs(sphere,varargin{:}); + +CSX = Add2Property(CSX,propName, sphere, 'SphericalShell');
\ No newline at end of file diff --git a/CSXCAD/matlab/AddWire.m b/CSXCAD/matlab/AddWire.m new file mode 100644 index 0000000..f90843c --- /dev/null +++ b/CSXCAD/matlab/AddWire.m @@ -0,0 +1,46 @@ +function CSX = AddWire(CSX, propName, prio, points, wire_rad, varargin) +% function CSX = AddWire(CSX, propName, prio, points, wire_rad, varargin) +% +% Add a wire to CSX and assign to a property with name <propName>. +% +% Warning: This is a 1D object, not all properties may be compatible with a +% 1D object, e.g. a material property. +% +% points: curve coordinates array +% prio : primitive priority +% wire_rad: wire radius +% +% example: +% %first point +% points(1,1) = 0; +% points(2,1) = 5; +% points(3,1) = 10; +% %second point +% points(1,2) = 0; +% points(2,2) = 10; +% points(3,2) = 10; +% %third point ... +% % create a metal wire with finite radius... +% CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +% CSX = AddCurve(CSX,'metal',10, points, 2); +% +% See also AddBox, AddCylindricalShell, AddCylinder, AddSphere, +% AddSphericalShell, AddCurve, AddMetal +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +wire.ATTRIBUTE.Priority = prio; +wire.ATTRIBUTE.WireRadius = wire_rad; + +wire.Vertex={}; +for s=1:size(points,2) + wire.Vertex{end+1}.ATTRIBUTE.X = points(1,s); + wire.Vertex{end}.ATTRIBUTE.Y = points(2,s); + wire.Vertex{end}.ATTRIBUTE.Z = points(3,s); +end + +wire = AddPrimitiveArgs(wire,varargin{:}); + +CSX = Add2Property(CSX,propName, wire, 'Wire');
\ No newline at end of file diff --git a/CSXCAD/matlab/AnalyseMesh.m b/CSXCAD/matlab/AnalyseMesh.m new file mode 100644 index 0000000..93295c3 --- /dev/null +++ b/CSXCAD/matlab/AnalyseMesh.m @@ -0,0 +1,42 @@ +function results = AnalyseMesh(lines) +% function results = AnalyseMesh(lines) +% +% Analyse a given mesh line vector +% +% output structure: +% results.numLines: number of lines +% results.max_res: max. resolution found +% results.min_res: min. resolution found +% results.max_ratio: max. grading ratio found +% results.homogeneous: true/false for homogeneous mesh +% results.symmetric: true/false for symmetric mesh +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (C) 2012 + +results = []; + +lines = sort(unique(lines)); + +if (numel(lines)<=1) + warning('CSXCAD:AnalyseMesh', 'more than one line needed to analyse mesh'); +end + +diff_lines = diff(lines); + +results.numLines = numel(lines); + +results.max_res = max(diff_lines); +results.min_res = min(diff_lines); +if (results.max_res==results.min_res) + results.homogeneous = 1; +else + results.homogeneous = 0; +end + +results.symmetric = CheckSymmtricLines(lines); + +ratio_lines = diff_lines(1:end-1)./diff_lines(2:end); + +results.max_ratio = max([ratio_lines 1./ratio_lines]); diff --git a/CSXCAD/matlab/AutoSmoothMeshLines.m b/CSXCAD/matlab/AutoSmoothMeshLines.m new file mode 100644 index 0000000..e41e261 --- /dev/null +++ b/CSXCAD/matlab/AutoSmoothMeshLines.m @@ -0,0 +1,159 @@ +function [lines quality] = AutoSmoothMeshLines( lines, max_res, ratio, varargin) +% function [lines quality] = AutoSmoothMeshLines( lines, max_res, ratio, varargin) +% +% Generate smooth mesh lines by choosing an appropriate algorithm. +% +% Currently supported algorithm: +% SmoothMeshLines, SmoothMeshLines2 and RecursiveSmoothMesh +% +% arguments: +% lines: given fixed lines to create a smooth mesh in between +% max_res: desired max. resolution +% ratio: grading ratio: desired neighboring line-delta ratio +% - default is 1.5 +% - see also 'allowed_max_ratio' argument +% +% variable arguments ('keyword',value): +% algorithm: define subset of tried algorihm, e.g. [1 3] +% symmetric: 0/1 force symmetric mesh (default is input symmetry) +% homogeneous: 0/1 force homogeneous mesh +% allowed_min_res: allow a given min resolution only +% allowed_max_ratio: allow only a given max. grading ratio +% (default --> ratio*1.25) +% debug: 0/1 off/on +% +% example: +% lines = AutoSmoothMeshLines([-100 -10 10 100], 20, 1.5, 'algorihm', ... +% 1:3); +% +% See also InitCSX, DefineRectGrid +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (C) 2012 + +if (nargin<2) + error('CSXCAD:AutoSmoothMeshLines','lines and max_res are a required parameter'); +end + +if (nargin<3) + ratio = 1.5; +end + +lines = sort(unique(lines)); + +range = lines(end)-lines(1); +if (~isempty(find(diff(lines)<range*1e-6))) + warning('CSXCAD:AutoSmoothMeshLines','some lines found with very small distance which may cause smoothing failure!'); +end + +methods = {}; +methods{end+1} = @SmoothMeshLines; +methods{end+1} = @SmoothMeshLines2; +methods{end+1} = @RecursiveSmoothMesh; + +requires_homogen = 0; +requires_symmetric = CheckSymmtricLines(lines); +allowed_min_res = 0; +debug = 0; +max_ratio = ratio*1.25; + +algorithm = 1:numel(methods); + +for vn=1:2:numel(varargin) + if (strcmpi(varargin{vn},'algorithm')) + algorithm = intersect(varargin{vn+1},algorithm); + end + if (strcmpi(varargin{vn},'symmetric')) + requires_symmetric = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'homogeneous')) + requires_homogen = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'force_min_res')) + requires_homogen = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'allowed_min_res')) + allowed_min_res = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'allowed_max_ratio')) + max_ratio = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'debug')) + debug = varargin{vn+1}; + end +end + +for m=algorithm + if (debug>0) + disp(['AutoSmoothMeshLines: trying method: ' func2str(methods{m})]); + tic + end + out_lines{m} = methods{m}(lines, max_res, ratio, 'CheckMesh', false); + if (debug>0) + toc + end + quality(m) = eval_mesh(out_lines{m}, max_res, max_ratio, requires_homogen, requires_symmetric, allowed_min_res, 1); + if (quality(m)==100) % uncomment to release! + lines = out_lines{m}; + if (debug>0) + disp(['AutoSmoothMeshLines: The winner with 100% is ' func2str(methods{m})]); % remove to release + end + return + end +end +winner = find(quality==max(quality),1); +lines = out_lines{winner}; +if (debug>0) + disp(['New_SmoothMeshLines: The winner with ' num2str(quality(winner)) '% is ' func2str(methods{winner})]); % remove to release +end + +% show mesh problems +eval_mesh(lines, max_res, ratio, requires_homogen, requires_symmetric, allowed_min_res, 0); +return + +error('CSXCAD:AutoSmoothMeshLines',['unknown algorithm requested: ' algorithm ]); + +end + +function quality = eval_mesh(lines, max_res, ratio, requires_homogen, requires_symmetric, allowed_min_res, silent) + +quality = 100; + +results = AnalyseMesh(lines); +if ((requires_homogen==1) && (results.homogeneous~=1)) + if (silent==0) + warning('CSXCAD:AutoSmoothMeshLines','method failed to create homogenous mesh'); + end + quality = -1; + return +end +if ((requires_symmetric==1) && (results.symmetric~=1)) + if (silent==0) + warning('CSXCAD:AutoSmoothMeshLines','method failed to create symmetric mesh'); + end + quality = -1; + return +end + +if ((allowed_min_res>0) && (results.min_res<allowed_min_res)) + if (silent==0) + warning('CSXCAD:AutoSmoothMeshLines','method failed to obey allowed min res!'); + end + quality = -1; + return +end + +if (results.max_res>max_res*1.01) + if (silent==0) + warning('CSXCAD:AutoSmoothMeshLines',['method failed to fulfill max. res: ' num2str(results.max_res) ' > ' num2str(max_res)]); + end + quality = quality*(max_res/results.max_res); +end +if (results.max_ratio>ratio*1.01) + if (silent==0) + warning('CSXCAD:AutoSmoothMeshLines',['method failed to fulfill the max. ratio: ' num2str(results.max_ratio) ' > ' num2str(ratio)]'); + end + quality = quality*(ratio/results.max_ratio); +end +end
\ No newline at end of file diff --git a/CSXCAD/matlab/CSXGeomPlot.m b/CSXCAD/matlab/CSXGeomPlot.m new file mode 100644 index 0000000..df33c94 --- /dev/null +++ b/CSXCAD/matlab/CSXGeomPlot.m @@ -0,0 +1,43 @@ +function CSXGeomPlot(CSX_filename, args_string) +% function CSXGeomPlot(CSX_filename<, args_string>) +% +% Show the geometry stored in the CSX file using AppCSXCAD +% +% Optional AppCSXCAD arguments (args_string): +% '--RenderDiscMaterial', enable material rendering +% +% exports: +% '--export-polydata-vtk=<path-for-export>' +% '--export-STL=<path-for-export>' +% +% See also InitCSX, DefineRectGrid +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +if nargin < 1 + error 'specify the xml file to open' +end + +if nargin < 2 + args_string = ''; +end + +filename = mfilename('fullpath'); +dir = fileparts( filename ); + +if isunix + AppCSXCAD_bin = searchBinary('AppCSXCAD.sh', ... + {[dir filesep '..' filesep '..' filesep 'AppCSXCAD' filesep], ... + [dir filesep '..' filesep '..' filesep '..' filesep 'bin' filesep]}); +else % assume windows + AppCSXCAD_bin = searchBinary('AppCSXCAD.exe',[dir filesep '..' filesep]); +end + +command = [AppCSXCAD_bin ' --disableEdit ' args_string ' ' CSX_filename]; +disp( ['invoking AppCSXCAD, exit to continue script...'] ); +if isOctave() + fflush(stdout); +end +system(command); diff --git a/CSXCAD/matlab/CalcDebyeMaterial.m b/CSXCAD/matlab/CalcDebyeMaterial.m new file mode 100644 index 0000000..a4a3a5e --- /dev/null +++ b/CSXCAD/matlab/CalcDebyeMaterial.m @@ -0,0 +1,29 @@ +function eps_debye = CalcDebyeMaterial(f, eps_r, kappa, eps_Delta, t_relax) +% eps_debye = CalcDebyeMaterial(f, eps_r, kappa, eps_Delta, t_relax) +% +% Calculate the Debye type dispersive material constant +% +% arguments: +% f: frequeny range of interest +% eps_r: eps_r infinity +% kappa: conductivity (losses) +% eps_Delta: (vector) delta of relative permitivity +% t_relax: (vector) relaxation time (losses) +% +% return: +% eps_debye: the complex relative permitivity +% +% See also: CalcLorentzMaterial +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (2013) + +EPS0 = 8.85418781762e-12; +eps_debye = ones(size(f))*eps_r - 1j*kappa./(2*pi*f)/EPS0; + +for n=1:numel(eps_Delta) + eps_debye = eps_debye + eps_Delta(n)./(1+2j*pi*f*t_relax(n)); +end + +end diff --git a/CSXCAD/matlab/CalcDrudeMaterial.m b/CSXCAD/matlab/CalcDrudeMaterial.m new file mode 100644 index 0000000..9901566 --- /dev/null +++ b/CSXCAD/matlab/CalcDrudeMaterial.m @@ -0,0 +1,37 @@ +function eps_drude = CalcDrudeMaterial(f, eps_r, kappa, plasmaFreq, t_relax) +% eps_drude = CalcDrudeMaterial(f, eps_r, kappa, plasmaFreq, t_relax) +% +% Calculate the Drude type dispersive material constant +% +% arguments: +% f: frequeny range of interest +% eps_r: eps_r infinity +% kappa: conductivity (losses) +% plasmaFreq: (vector) plasma frequencies +% t_relax: (vector) relaxation time (losses) +% +% return: +% eps_drude: the complex relative permitivity +% +% Example: +% % silver (AG) at optical frequencies (Drude model) [1, p. 201] +% f = linspace(300e12, 1100e12, 201); +% eps_model = CalcDrudeMaterial(f, 3.942, 7.97e3, 7e15/2/pi, 0, 1/2.3e13); +% +% figure +% plot(f,real(eps_model)) +% hold on; +% grid on; +% plot(f,imag(eps_model),'r--') +% +% See also: CalcDebyeMaterial, CalcLorentzMaterial +% +% [1] Rennings, Andre: "Elektromagnetische Zeitbereichssimulationen +% innovativer Antennen auf Basis von Metamaterialien." +% PhD Thesis, University of Duisburg-Essen, September 2008 +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (2013) + +eps_drude = CalcLorentzMaterial(f, eps_r, kappa, plasmaFreq, plasmaFreq*0, t_relax); diff --git a/CSXCAD/matlab/CalcLorentzMaterial.m b/CSXCAD/matlab/CalcLorentzMaterial.m new file mode 100644 index 0000000..2785747 --- /dev/null +++ b/CSXCAD/matlab/CalcLorentzMaterial.m @@ -0,0 +1,61 @@ +function eps_lorentz = CalcLorentzMaterial(f, eps_r, kappa, plasmaFreq, LorPoleFreq, t_relax) +% eps_lorentz = CalcLorentzMaterial(f, eps_r, kappa, plasmaFreq, LorPoleFreq, t_relax) +% +% Calculate the Lorentz type dispersive material constant +% +% arguments: +% f: frequeny range of interest +% eps_r: eps_r infinity +% kappa: conductivity (losses) +% plasmaFreq: (vector) plasma frequencies (Drude model) +% LorPoleFreq: (vector) Lorentz pole frequencies (zero for pure Drude model) +% t_relax: (vector) relaxation time (losses) +% +% return: +% eps_lorentz: the complex relative permitivity +% +% Example: +% % silver (AG) at optical frequencies (Drude model) [1, p. 201] +% f = linspace(300e12, 1100e12, 201); +% eps_model = CalcLorentzMaterial(f, 3.942, 7.97e3, 7e15/2/pi, 0, 1/2.3e13); +% +% figure +% plot(f,real(eps_model)) +% hold on; +% grid on; +% plot(f,imag(eps_model),'r--') +% +% % silver (AG) at optical frequencies (Drude+Lorentz model) [1, p. 201] +% f = linspace(300e12, 1100e12, 201); +% eps_model = CalcLorentzMaterial(f, 1.138, 4.04e3, [13e15 9.61e15]/2/pi, [0 7.5e15]/2/pi,[1/2.59e13 1/3e14]); +% +% figure +% plot(f,real(eps_model)) +% hold on; +% grid on; +% plot(f,imag(eps_model),'r--') +% +% See also: CalcDebyeMaterial +% +% [1] Rennings, Andre: "Elektromagnetische Zeitbereichssimulationen +% innovativer Antennen auf Basis von Metamaterialien." +% PhD Thesis, University of Duisburg-Essen, September 2008 +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (2013) + +EPS0 = 8.85418781762e-12; +eps_lorentz = ones(size(f))*eps_r - 1j*kappa./(2*pi*f)/EPS0; + +w = 2*pi*f; +for n=1:numel(plasmaFreq) + if (t_relax(n)>0) + w_r = 1/t_relax(n); + else + w_r = 0; + end + eps_lorentz = eps_lorentz - eps_r*(2*pi*plasmaFreq(n))^2./(w.^2 - (2*pi*LorPoleFreq(n))^2 - 2j*pi*f*w_r); +end + +end diff --git a/CSXCAD/matlab/CheckMesh.m b/CSXCAD/matlab/CheckMesh.m new file mode 100644 index 0000000..68cac0b --- /dev/null +++ b/CSXCAD/matlab/CheckMesh.m @@ -0,0 +1,65 @@ +function [EC pos E_type] = CheckMesh(lines, min_res, max_res, ratio, be_quiet) +% function [EC pos E_type] = CheckMesh(lines, min_res, max_res, ratio, be_quiet) +% +% Check if mesh lines are valid +% +% parameter: +% min_res: minimal allowed mesh-diff +% max_res: maximal allowed mesh-diff +% ratio: maximal allowed mesh-diff ratio +% be_quiet: disable warnings +% +% return: +% EC: error code (number of found errors) +% pos: line positions with error +% E_type: error type +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +if (nargin<5) + be_quiet = 0; +end + +diff_lines = diff(lines); +EC = 0; +E_type = []; + +pos = []; +max_err = find(diff_lines>max_res); +if (~isempty(max_err) && ~be_quiet) + warning('CSXCAD:CheckMesh','found resolution larger than max_res'); + pos = [pos max_err]; + EC = EC + numel(max_err); + E_type = [E_type 1]; +end + +min_err = find(diff_lines<min_res); +if (~isempty(min_err) && ~be_quiet) + warning('CSXCAD:CheckMesh','found resolution smaller than min_res'); + pos = [pos min_err]; + EC = EC + numel(min_err); + E_type = [E_type 2]; +end + +for n=1:numel(diff_lines)-1 + if (diff_lines(n+1)/diff_lines(n) > ratio*1.01) + str = ['lines: ' num2str(n) '@' num2str(lines(n)) ' ' num2str(n+1) '@' num2str(lines(n+1)) ' ' num2str(n+2) '@' num2str(lines(n+2))]; + if (~be_quiet) + warning('CSXCAD:CheckMesh', [str '\nfound resolution increase larger than ratio: ' num2str(diff_lines(n+1)/diff_lines(n)) ' > ' num2str(ratio) ]); + end + pos = [pos n+1]; + EC = EC + 1; + E_type = [E_type 3]; + end + if (diff_lines(n+1)/diff_lines(n) < 1/ratio/1.01) + str = ['lines: ' num2str(n) '@' num2str(lines(n)) ' ' num2str(n+1) '@' num2str(lines(n+1)) ' ' num2str(n+2) '@' num2str(lines(n+2))]; + if (~be_quiet) + warning('CSXCAD:SmoothRange', [str '\nfound resolution decrease smaller than ratio: ' num2str(diff_lines(n+1)/diff_lines(n)) ' < 1/' num2str(ratio) '=' num2str(1/ratio) ]); + end + pos = [pos n]; + EC = EC + 1; + E_type = [E_type 4]; + end +end diff --git a/CSXCAD/matlab/Convert_VF_DiscMaterial.m b/CSXCAD/matlab/Convert_VF_DiscMaterial.m new file mode 100644 index 0000000..5d78451 --- /dev/null +++ b/CSXCAD/matlab/Convert_VF_DiscMaterial.m @@ -0,0 +1,255 @@ +function Convert_VF_DiscMaterial(raw_filesuffix, mat_db_file, out_filename, varargin) +% function Convert_VF_DiscMaterial(raw_filesuffix, mat_db_file, out_filename, varargin) +% +% function to convert the virtual family raw voxel data to a disc material +% property required by CSXCAD/openEMS +% +% Note: The conversion is (currently) done, converting the broad-band data +% into a relative permittivity and conductivity for a given frequency of +% interest. Thus the converted model is only valid within a small frequency +% range around the used conversion frequency. +% +% required arguments: +% raw_filesuffix: suffix for the virtual family body model files +% the files: +% <raw_filesuffix>.txt and <raw_filesuffix>.raw +% must be found! +% example: '/tmp/Ella_26y_V2_1mm' +% mat_db_file: tissue database file (incl. full path if necessary) +% example: '/tmp/DB_h5_20120711_SEMCADv14.8.h5' +% out_filename: outfile name, e.g. 'Ella_298MHz.h5' +% +% variable arguments (key/value): +% 'Frequency': specifiy the frequency of interest (required!) +% 'Center': 0/1 make the model centered around (0,0,0) +% (default is off) +% +% Requirements: +% - Matlab, Octave is currently not supported due to missing hdf5 write +% functions +% - ~6GB of RAM to convert the largest voxel model +% - Virtual Family voxel models +% e.g. Duke_34y_V5_1mm.raw and .txt +% See http://www.itis.ethz.ch/itis-for-health/virtual-population/overview/ +% - Virtual Family tissue database +% e.g. DB_h5_20120711_SEMCADv14.8.h5 +% Download from: http://www.itis.ethz.ch/assets/Downloads/TissueDb/DBh5_20120711_SEMCADv14.8.zip +% +% example: +% Convert_VF_DiscMaterial('/tmp/Ella_26y_V2_1mm', ... +% '/tmp/DB_h5_20120711_SEMCADv14.8.h5', ... +% 'Ella_centered_298MHz.h5', ... +% 'Frequency', 2.4e9); +% +% (c) 2013 Thorsten Liebig + +if exist(out_filename,'file') + disp(['requested model file: "' out_filename '" already exists, skipping']); + return +end + +%% default settings +center_Model = 0; % make the model centered around (0,0,0) +freq = nan; % frequency of interest +range = {[],[],[]}; + +% internal +range_used = 0; + +for n=1:2:numel(varargin) + if (strcmp(varargin{n},'Frequency')) + freq = varargin{n+1}; + elseif (strcmp(varargin{n},'Center')) + center_Model = varargin{n+1}; + elseif (strcmp(varargin{n},'Range')) + range = varargin{n+1}; + range_used = 1; + else + warning('CSXCAD:Convert_VF_DiscMaterial',['unknown argument: ' varagin{n}]); + end +end + +if (isnan(freq)) + error('CSXCAD:Convert_VF_DiscMaterial','a frequency of interest must be specified'); +end + +%% end of settings (do not edit below) +physical_constants +w = 2*pi*freq; + +%% +disp(['reading raw body specs: ' raw_filesuffix '.txt']); +fid = fopen([raw_filesuffix '.txt']); +material_list = textscan(fid,'%d %f %f %f %s'); + +frewind(fid) +feof(fid); +while ~feof(fid) + line = fgetl(fid); + if strcmp(line,'Grid extent (number of cells)') + n_cells = textscan(fid,'n%c %f'); + end + if strcmp(line,'Spatial steps [m]') + d_cells = textscan(fid,'d%c %f'); + end +end + +fclose(fid); + +nx = n_cells{2}(1); +ny = n_cells{2}(2); +nz = n_cells{2}(3); + +disp(['body model contains ' num2str(nx*ny*nz) ' voxels']); + +dx = d_cells{2}(1); +dy = d_cells{2}(2); +dz = d_cells{2}(3); + +% +x = (0:nx)*dx; +y = (0:ny)*dy; +z = (0:nz)*dz; + +if (center_Model>0) + disp('centering model...'); + x = x - x(1) - (max(x)-min(x))/2; + y = y - y(1) - (max(y)-min(y))/2; + z = z - z(1) - (max(z)-min(z))/2; +end + +disp('Body model mesh range (original):'); +disp(['x: ' num2str(x(1)) ' --> ' num2str(x(end)) ', width: ' num2str(max(x)-min(x))]); +disp(['y: ' num2str(y(1)) ' --> ' num2str(y(end)) ', width: ' num2str(max(y)-min(y))]); +disp(['z: ' num2str(z(1)) ' --> ' num2str(z(end)) ', width: ' num2str(max(z)-min(z))]); + +%% map to range +if (isempty(range{1})) + x_idx = 1:nx; +else + x_idx = find(x>=range{1}(1) & x<=range{1}(2)); +end + +if (isempty(range{2})) + y_idx = 1:ny; +else + y_idx = find(y>=range{2}(1) & y<=range{2}(2)); +end + +if (isempty(range{3})) + z_idx = 1:nz; +else + z_idx = find(z>=range{3}(1) & z<=range{3}(2)); +end + +tic +disp(['reading raw body data: ' raw_filesuffix '.raw']); +fid = fopen([raw_filesuffix '.raw']); + +skip=(z(z_idx(1))-z(1))/dz*nx*ny; +fseek(fid,skip,'bof'); + +% read in line by line to save memory +for n=1:length(z_idx) + data_plane = reshape(fread(fid,nx*ny),nx,ny); + mat(1:numel(x_idx),1:numel(y_idx),n) = uint8(data_plane(x_idx,y_idx)); +end +fclose(fid); + +% resize to range +x = x(x_idx);x=[x x(end)+dx]; +y = y(y_idx);y=[y y(end)+dx]; +z = z(z_idx);z=[z z(end)+dx]; + +if (range_used) + disp('Body model mesh range (new):'); + disp(['x: ' num2str(x(1)) ' --> ' num2str(x(end)) ', width: ' num2str(max(x)-min(x))]); + disp(['y: ' num2str(y(1)) ' --> ' num2str(y(end)) ', width: ' num2str(max(y)-min(y))]); + disp(['z: ' num2str(z(1)) ' --> ' num2str(z(end)) ', width: ' num2str(max(z)-min(z))]); +end + +nx = numel(x_idx); +ny = numel(y_idx); +nz = numel(z_idx); + +mesh.x=x; +mesh.y=y; +mesh.z=z; + +toc + +%% +disp(['reading/analysing material database: ' mat_db_file]); + +mat_db_info = h5info(mat_db_file); + +%% +% there is no material in the list for number 0 +mat_db.epsR(1) = 1; +mat_db.kappa(1) = 0; +mat_db.density(1) = 0; +mat_db.Name{1} = 'Background'; + +for n=1:numel(material_list{end}) + mat_para = GetColeData(mat_db_info, material_list{end}{n}); + eps_colo_cole = mat_para.ef + mat_para.sig/(1j*w*EPS0); + if (mat_para.del1>0) + eps_colo_cole = eps_colo_cole + mat_para.del1/(1 + (1j*w*mat_para.tau1*1e-12)^(1-mat_para.alf1)); + end + if (mat_para.del2>0) + eps_colo_cole = eps_colo_cole + mat_para.del2/(1 + (1j*w*mat_para.tau2*1e-9)^(1-mat_para.alf2)); + end + if (mat_para.del3>0) + eps_colo_cole = eps_colo_cole + mat_para.del3/(1 + (1j*w*mat_para.tau3*1e-6)^(1-mat_para.alf3)); + end + if (mat_para.del4>0) + eps_colo_cole = eps_colo_cole + mat_para.del4/(1 + (1j*w*mat_para.tau4*1e-3)^(1-mat_para.alf4)); + end + mat_db.epsR(n+1) = real(eps_colo_cole); + mat_db.kappa(n+1) = -1*imag(eps_colo_cole)*w*EPS0; + mat_db.density(n+1) = mat_para.Dens; + mat_db.Name{n+1} = mat_para.Name; +end +toc +%% + +disp(['Creating DiscMaterial file: ' out_filename]); +CreateDiscMaterial(out_filename, mat, mat_db, mesh) +toc + +disp('done..'); + +end + +function [mat_para] = GetColeData(mat_db, mat_name) + +% search for mat_name in all "AlternNames" + +N_grp = numel(mat_db.Groups(1).Groups); + +mat_name = regexp(mat_name, '/','split'); +mat_name = mat_name{end}; + +for n=1:N_grp + found = false; + name = regexp(mat_db.Groups(1).Groups(n).Name, '/','split'); + alt_names = [name(end) regexp(mat_db.Groups(1).Groups(n).Attributes(end).Value, '/','split')]; + + alt_names = regexprep(alt_names,' ','_'); + + if (sum(strcmpi(alt_names, mat_name))>0) + % we found it, break + found = true; + break; + end +end + +if (found==false) + error(['property "' mat_name '" not found']); +end + +for a = 1:numel(mat_db.Groups(1).Groups(n).Attributes) + mat_para.(mat_db.Groups(1).Groups(n).Attributes(a).Name) = mat_db.Groups(1).Groups(n).Attributes(a).Value; +end + +end diff --git a/CSXCAD/matlab/CreateDiscMaterial.m b/CSXCAD/matlab/CreateDiscMaterial.m new file mode 100644 index 0000000..f874230 --- /dev/null +++ b/CSXCAD/matlab/CreateDiscMaterial.m @@ -0,0 +1,83 @@ +function CreateDiscMaterial(filename, data, mat_db, mesh) +% function CreateDiscMaterial(filename, data, mat_db, mesh) +% +% Create the discrete material hdf5 file (version 2) usable by AddDiscMaterial +% +% Note: This function currently requires Matlab. Octave is missing the +% necessary hdf5 write functions. +% +% arguments: +% filename: hdf5 file to create (must not exist) +% data: voxel based index data, index as used in mat_db-1 +% mat_db: material database +% mesh: used voxel mesh. Note size is size(data)+[1 1 1] +% +% example: +% mat_db.epsR = [1 3 4]; %relative permittivity +% mat_db.kappa = [0 0.2 0.4]; %electric conductivity (S/m) +% mat_db.density = [0 1000 1010]; %material density (kg/m³) +% mat_db.Name = {'Background','mat2','mat3'}; +% +% data = [0 1 0; 2 1 2; 0 1 0]; % 3x3x3 data +% mesh.x = [0 0.1 0.2 0.3]; % 4 mesh lines in x-dir (3 cells) +% mesh.y = [-0.1 0 0.1 0.2]; % 4 mesh lines in y-dir (3 cells) +% mesh.z = [0 0.4 0.8 1.2]; % 4 mesh lines in z-dir (3 cells) +% +% CreateDiscMaterial('test_mat.h5', data, mat_db, mesh); +% +% See also AddDiscMaterial +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (2013) + +if isOctave + error('CreateDiscMaterial currently does not support Octave, due to missing hdf5 functions!'); +end + +if (exist('h5create')==0) + error('CSXCAD:CreateDiscMaterial','Your matlab seems to be too old, h5create cannot be found!'); +end + +data_size = size(data); +mesh_size = [numel(mesh.x) numel(mesh.y) numel(mesh.z)]; + +if ((mesh_size-1)~=data_size) + error('data size and mesh size mismatch'); +end + +if (exist(filename,'file')) + error(['file "' filename '" already exist. Delete/rename first!']); +end + +h5create(filename, '/DiscData',data_size, 'Datatype', 'uint8', 'ChunkSize',data_size, 'Deflate',9); +h5write(filename, '/DiscData', data); + +clear data; + +h5writeatt(filename, '/DiscData', 'DB_Size', int32(numel(mat_db.epsR))); + +h5writeatt(filename, '/DiscData', 'epsR', mat_db.epsR); +h5writeatt(filename, '/DiscData', 'kappa', mat_db.kappa); +h5writeatt(filename, '/DiscData', 'density', mat_db.density); +h5writeatt(filename, '/DiscData', 'Name', strjoin(mat_db.Name,',')); + +h5create(filename, '/mesh/x', mesh_size(1)); +h5write(filename, '/mesh/x', mesh.x); + +h5create(filename, '/mesh/y', mesh_size(2)); +h5write(filename, '/mesh/y', mesh.y); + +h5create(filename, '/mesh/z', mesh_size(3)); +h5write(filename, '/mesh/z', mesh.z); + +h5writeatt(filename, '/', 'Version', 2); + + + +function out = strjoin(names, delimitier) + +out = names{1}; +for n=2:numel(names) + out = [out delimitier names{n}]; +end diff --git a/CSXCAD/matlab/DefineRectGrid.m b/CSXCAD/matlab/DefineRectGrid.m new file mode 100644 index 0000000..efc8a41 --- /dev/null +++ b/CSXCAD/matlab/DefineRectGrid.m @@ -0,0 +1,50 @@ +function CSX = DefineRectGrid(CSX, deltaUnit, mesh) +% function CSX = DefineRectGrid(CSX, deltaUnit, mesh); +% +% Create a rectiliniear grid. +% +% example Cartesian mesh: +% CSX = InitCSX(); +% mesh.x = AutoSmoothMeshLines([0 a], 10); +% mesh.y = AutoSmoothMeshLines([0 b], 10); +% mesh.z = AutoSmoothMeshLines([0 length], 15); +% CSX = DefineRectGrid(CSX, unit,mesh); +% +% example Cylindrical mesh: +% CSX = InitCSX('CoordSystem',1); +% mesh.r = AutoSmoothMeshLines([0 a], 10); +% mesh.a = AutoSmoothMeshLines([0 2*pi], pi/30); +% mesh.z = AutoSmoothMeshLines([-length 0 length], 15); +% CSX = DefineRectGrid(CSX, unit,mesh); +% +% See also InitCSX, SmoothMesh, AutoSmoothMeshLines, DetectEdges +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit = deltaUnit; + +if (isfield(CSX,'ATTRIBUTE')) + if (isfield(CSX.ATTRIBUTE,'CoordSystem')) + CSX.RectilinearGrid.ATTRIBUTE.CoordSystem = CSX.ATTRIBUTE.CoordSystem; + end +end + +if (isfield(mesh,'x')) + CSX.RectilinearGrid.XLines = mesh.x; +elseif ( (isfield(mesh,'r')) && (CSX.ATTRIBUTE.CoordSystem==1)) + CSX.RectilinearGrid.XLines = mesh.r; +else + error 'x/(r) direction not found' +end + +if (isfield(mesh,'y')) + CSX.RectilinearGrid.YLines = mesh.y; +elseif ((isfield(mesh,'a')) && (CSX.ATTRIBUTE.CoordSystem==1)) + CSX.RectilinearGrid.YLines = mesh.a; +else + error 'y/(a) direction not found' +end + +CSX.RectilinearGrid.ZLines = mesh.z; diff --git a/CSXCAD/matlab/DetectEdges.m b/CSXCAD/matlab/DetectEdges.m new file mode 100644 index 0000000..ea01909 --- /dev/null +++ b/CSXCAD/matlab/DetectEdges.m @@ -0,0 +1,278 @@ +function mesh = DetectEdges(CSX, mesh, varargin) +% mesh = DetectEdges(CSX <, mesh, varargin>) +% +% Returns a mesh with the edges added of certain (see below) primitives +% found in the CSX structure. +% +% optional arguments: +% mesh: add edges to an existing (predefined) mesh +% +% optional: 'keyword', value +% 'Debug' enable debug mode (default is 1) +% 'AddPropertyType' add a list of additional property types to detect +% e.g. 'DumpBox' or {'DumpBox','ProbeBox'} +% 'SetPropertyType' set the list of property types to detect (override default) +% e.g. 'Metal' or {'Metal','ConductingSheet'} +% 'ExcludeProperty' give a list of property names to exclude from +% detection +% 'SetProperty' give a list of property names to handly exlusively for detection +% +% advanced options: 'keyword', value +% '2D_Metal_Edge_Res' define a one-third/two-third metal edge resolution +% +% example: +% CSX = InitCSX(); +% % define all properties and primitives to detect +% % ... +% % ... +% mesh = DetectEdges(CSX); +% mesh = SmoothMesh(mesh, lambda/20/unit); +% CSX = DefineRectGrid(CSX, unit, mesh); +% +% Note: +% - Only primitives contained in Metal, Material, Excitation, LumpedElement +% or ConductingSheet properties are processed +% - Currently this function handles only Box, Polygons, LinPoly and Cylinder. +% +% CSXCAD matlab interface +% ----------------------- +% author: Koen De Vleeschauwer (c) 2012 +% modified by: Thorsten Liebig (c) 2012, 2013 +% +% See also InitCSX, SmoothMesh, DefineRectGrid + +supported_properties = {}; +supported_properties{end+1}='Metal'; +supported_properties{end+1}='Material'; +supported_properties{end+1}='Excitation'; +supported_properties{end+1}='LumpedElement'; +supported_properties{end+1}='ConductingSheet'; + +exclude_list = {}; +prop_list_only = {}; + +debug = 1; +mesh_2D_metal_edges = 0; + +for n=1:2:numel(varargin) + if (strcmpi(varargin{n},'debug')==1); + debug = varargin{n+1}; + elseif (strcmpi(varargin{n},'AddPropertyType')==1); + if iscell(varargin{n+1}) + supported_properties(end+1) = varargin{n+1}; + elseif ischar(varargin{n+1}) + supported_properties{end+1} = varargin{n+1}; + else + error('CSXCAD:DetectEdges','unknown property definition'); + end + elseif (strcmpi(varargin{n},'SetPropertyType')==1); + if iscell(varargin{n+1}) + supported_properties = varargin{n+1}; + elseif ischar(varargin{n+1}) + supported_properties = {varargin{n+1}}; + else + error('CSXCAD:DetectEdges','unknown property definition'); + end + elseif (strcmpi(varargin{n},'ExcludeProperty')==1); + exclude_list = varargin{n+1}; + elseif (strcmpi(varargin{n},'SetProperty')==1); + prop_list_only = varargin{n+1}; + elseif (strcmpi(varargin{n},'2D_Metal_Edge_Res')==1); + mesh_2D_metal_edges = varargin{n+1}*[-1 2]/3; + else + warning('CSXCAD:DetectEdges',['unknown argument: ' varargin{n}]); + end +end + +if (nargin<2) + mesh = []; +end + +if isempty(mesh) + if (CSX.ATTRIBUTE.CoordSystem==0) + mesh.x = []; + mesh.y = []; + mesh.z = []; + elseif (CSX.ATTRIBUTE.CoordSystem==1) + mesh.r = []; + mesh.a = []; + mesh.z = []; + else + error('CSXCAD:DetectEdges','unknown coordinate system used'); + end +end + +edges = {}; +edges.x = [ ]; +edges.y = [ ]; +edges.z = [ ]; + +if (~isstruct(CSX)) + error('expected a CSX structure'); +end + +CoordSystem = CSX.ATTRIBUTE.CoordSystem; + +if (isfield(CSX, 'Properties')) + prop_fn = fieldnames(CSX.Properties); + for p = 1:numel(prop_fn) + if (sum(strcmpi(prop_fn{p}, supported_properties))==0) + continue; + end + isMetal = sum(strcmpi(prop_fn{p},{'Metal','ConductingSheet'})); + property_group = CSX.Properties.(prop_fn{p}); + for m = 1:numel(property_group) + property=property_group{m}; + if ~isfield(property, 'Primitives') + continue; + end + if (sum(strcmpi(property.ATTRIBUTE.Name,exclude_list))) + continue; + end + if (~isempty(prop_list_only) && (sum(strcmpi(property.ATTRIBUTE.Name,prop_list_only))==0)) + continue; + end + primitives = property.Primitives; + prim_fn = fieldnames(primitives); + for n_prim = 1:numel(prim_fn) + if (strcmp(prim_fn{n_prim}, 'Box')) + for b = 1:length(primitives.Box) + box = primitives.Box{b}; + x1 = box.P1.ATTRIBUTE.X; + y1 = box.P1.ATTRIBUTE.Y; + z1 = box.P1.ATTRIBUTE.Z; + x2 = box.P2.ATTRIBUTE.X; + y2 = box.P2.ATTRIBUTE.Y; + z2 = box.P2.ATTRIBUTE.Z; + % check dimension + dim = (x1~=x2) + (y1~=y2) + (z1~=z2); + if ((dim==2) && isMetal) + % add to global list of edges, with a given 2D + % edge resolution + edges = AddEdge (edges, box, x1+sign(x1-x2)*mesh_2D_metal_edges, y1+sign(y1-y2)*mesh_2D_metal_edges, ... + z1+sign(z1-z2)*mesh_2D_metal_edges, CoordSystem, debug); + edges = AddEdge (edges, box, x2+sign(x2-x1)*mesh_2D_metal_edges, y2+sign(y2-y1)*mesh_2D_metal_edges, ... + z2+sign(z2-z1)*mesh_2D_metal_edges, CoordSystem, debug); + else + % add to global list of edges + edges = AddEdge (edges, box, x1, y1, z1, CoordSystem, debug); + edges = AddEdge (edges, box, x2, y2, z2, CoordSystem, debug); + end + end + elseif (strcmp(prim_fn{n_prim}, 'LinPoly') || strcmp(prim_fn{n_prim}, 'Polygon')) + for l = 1:length(primitives.(prim_fn{n_prim})) + poly = primitives.(prim_fn{n_prim}){l}; + dir = poly.ATTRIBUTE.NormDir + 1; + dirP = mod(poly.ATTRIBUTE.NormDir+1,3) + 1; + dirPP = mod(poly.ATTRIBUTE.NormDir+2,3) + 1; + lin_length = 0; + if (strcmp(prim_fn{n_prim}, 'LinPoly')) + lin_length = poly.ATTRIBUTE.Length; + end + if (isfield(poly, 'Vertex')) + for v = 1:length(poly.Vertex) + vertex = poly.Vertex{v}; + edge(dir) = poly.ATTRIBUTE.Elevation; + edge(dirP) = vertex.ATTRIBUTE.X1; + edge(dirPP) = vertex.ATTRIBUTE.X2; + edges = AddEdge (edges, poly, edge(1), edge(2), edge(3), CoordSystem, debug); + if (lin_length~=0) + edge(dir) = edge(dir) + lin_length; + edges = AddEdge (edges, poly, edge(1), edge(2), edge(3), CoordSystem, debug); + end + end + end + end + elseif (strcmp(prim_fn{n_prim}, 'Cylinder')) + for c = 1:length(primitives.Cylinder) + cylinder = primitives.Cylinder{c}; + r = cylinder.ATTRIBUTE.Radius; + x1 = cylinder.P1.ATTRIBUTE.X; + y1 = cylinder.P1.ATTRIBUTE.Y; + z1 = cylinder.P1.ATTRIBUTE.Z; + x2 = cylinder.P2.ATTRIBUTE.X; + y2 = cylinder.P2.ATTRIBUTE.Y; + z2 = cylinder.P2.ATTRIBUTE.Z; + if ((x1 == x2) && (y1 == y2) && (z1 ~= z2)) + % cylinder parallel with z axis + edges = AddEdge (edges, cylinder, x1 - r, y1 - r, z1, CoordSystem, debug); + edges = AddEdge (edges, cylinder, x2 + r, y2 + r, z2, CoordSystem, debug); + elseif ((x1 == x2) && (y1 ~= y2) && (z1 == z2)) + % cylinder parallel with y axis + edges = AddEdge (edges, cylinder, x1 - r, y1, z1 - r, CoordSystem, debug); + edges = AddEdge (edges, cylinder, x2 + r, y2, z2 + r, CoordSystem, debug); + elseif ((x1 ~= x2) && (y1 == y2) && (z1 == z2)) + % cylinder parallel with x axis + edges = AddEdge (edges, cylinder, x1, y1 - r, z1 - r, CoordSystem, debug); + edges = AddEdge (edges, cylinder, x2, y2 + r, z2 + r, CoordSystem, debug); + elseif (debug > 0) + warning('CSXCAD:DetectEdges',['unsupported primitive of type: "' prim_fn{n_prim} '" found, skipping edges']); + end + end + else + if (debug>0) + warning('CSXCAD:DetectEdges',['unsupported primitive of type: "' prim_fn{n_prim} '" found, skipping edges']); + end + end + end + end + end +end + +if (CSX.ATTRIBUTE.CoordSystem==0) + mesh.x = sort(unique([mesh.x edges.x])); + mesh.y = sort(unique([mesh.y edges.y])); + mesh.z = sort(unique([mesh.z edges.z])); +elseif (CSX.ATTRIBUTE.CoordSystem==1) + mesh.r = sort(unique([mesh.r edges.x])); + mesh.a = sort(unique([mesh.a edges.y])); + mesh.z = sort(unique([mesh.z edges.z])); +else + error('CSXCAD:DetectEdges','unknown coordinate system used'); +end + +end + + +function edges = AddEdge(edges, csx_prim, x, y, z, CoordSystem, debug) +% Add edges of CSX primitives including some transformations + +xt = unique(x); +yt = unique(y); +zt = unique(z); + +if isfield(csx_prim.ATTRIBUTE,'CoordSystem') + if (csx_prim.ATTRIBUTE.CoordSystem~=CoordSystem) + if (debug>2) + warning('CSXCAD:DetectEdges','different coordinate systems not supported, skipping edges'); + end + return + end +end + +if (isfield(csx_prim, 'Transformation')) + + transformation = csx_prim.Transformation; + trans_fn = fieldnames(transformation); + + for t=1:numel(trans_fn) + if (strcmp(trans_fn{t}, 'Translate')) + xt = xt + transformation.Translate.ATTRIBUTE.Argument(1); + yt = yt + transformation.Translate.ATTRIBUTE.Argument(2); + zt = zt + transformation.Translate.ATTRIBUTE.Argument(3); + else + if (debug>0) + warning('CSXCAD:DetectEdges','unsupported transformation found in primitive, skipping edges'); + end + return + end + end + +end + +% add to global list of edges +edges.x = [edges.x xt]; +edges.y = [edges.y yt]; +edges.z = [edges.z zt]; +end + diff --git a/CSXCAD/matlab/DirChar2Int.m b/CSXCAD/matlab/DirChar2Int.m new file mode 100644 index 0000000..0454448 --- /dev/null +++ b/CSXCAD/matlab/DirChar2Int.m @@ -0,0 +1,31 @@ +function n = DirChar2Int(dir_char) +% function n = DirChar2Int(dir_char) +% +% internal function to convert a character like 'x','y','z' into a numeric +% direction: 0..2! +% If input already is a numeric value from 0..2, it will just be copied! +% Everything else will raise an error! +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (c) 2013 + +if (numel(dir_char)>1) + error('CSXCAD:DirChar2Int','invalid normal direction') +end + +if (ischar(dir_char)) + if (strcmp(dir_char,'x') || strcmp(dir_char,'r')) + n = 0; return; + elseif (strcmp(dir_char,'y') || strcmp(dir_char,'a')) + n = 1; return; + elseif strcmp(dir_char,'z') + n = 2; return; + else + error('CSXCAD:DirChar2Int','invalid normal direction') + end +elseif (isnumeric(dir_char) && ((dir_char==0) || (dir_char==1) || (dir_char==2))) + n = dir_char; +else + error('CSXCAD:DirChar2Int','invalid normal direction') +end diff --git a/CSXCAD/matlab/ImportPLY.m b/CSXCAD/matlab/ImportPLY.m new file mode 100644 index 0000000..a221b38 --- /dev/null +++ b/CSXCAD/matlab/ImportPLY.m @@ -0,0 +1,23 @@ +function CSX = ImportPLY(CSX, propName, prio, filename, varargin) +% function CSX = ImportPLY(CSX, propName, prio, filename, varargin) +% +% example: +% CSX = AddMetal( CSX, 'cad_model' ); % create a perfect electric conductor (PEC) +% CSX = ImportPLY(CSX, 'cad_model',10, 'sphere.ply','Transform',{'Scale',1/unit}); +% +% Note: make sure the file 'sphere.ply' is in the working directory +% +% See also AddBox, AddCylinder, AddCylindricalShell, AddSphere, AddSphericalShell, +% AddCurve, AddWire, AddMetal, ImportSTL +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +plyfile.ATTRIBUTE.Priority = prio; +plyfile.ATTRIBUTE.FileName = filename; +plyfile.ATTRIBUTE.FileType = 'PLY'; + +plyfile = AddPrimitiveArgs(plyfile,varargin{:}); + +CSX = Add2Property(CSX,propName, plyfile, 'PolyhedronReader'); diff --git a/CSXCAD/matlab/ImportSTL.m b/CSXCAD/matlab/ImportSTL.m new file mode 100644 index 0000000..a3641f1 --- /dev/null +++ b/CSXCAD/matlab/ImportSTL.m @@ -0,0 +1,23 @@ +function CSX = ImportSTL(CSX, propName, prio, filename, varargin) +% function CSX = ImportSTL(CSX, propName, prio, filename, varargin) +% +% example: +% CSX = AddMetal( CSX, 'cad_model' ); % create a perfect electric conductor (PEC) +% CSX = ImportSTL(CSX, 'cad_model',10, 'sphere.stl','Transform',{'Scale',1/unit}); +% +% Note: make sure the file 'sphere.stl' is in the working directory +% +% See also AddBox, AddCylinder, AddCylindricalShell, AddSphere, AddSphericalShell, +% AddCurve, AddWire, AddMetal, ImportPLY +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +stlfile.ATTRIBUTE.Priority = prio; +stlfile.ATTRIBUTE.FileName = filename; +stlfile.ATTRIBUTE.FileType = 'STL'; + +stlfile = AddPrimitiveArgs(stlfile,varargin{:}); + +CSX = Add2Property(CSX,propName, stlfile, 'PolyhedronReader'); diff --git a/CSXCAD/matlab/InitCSX.m b/CSXCAD/matlab/InitCSX.m new file mode 100644 index 0000000..2411d85 --- /dev/null +++ b/CSXCAD/matlab/InitCSX.m @@ -0,0 +1,31 @@ +function CSX = InitCSX(varargin) +% function CSX = InitCSX() +% +% Inititalize the CSX data-structure. +% +% variable arguments: +% 'CoordSystem' : define the default coordinate system +% 0 -> Cartesian +% 1 -> Cylindircal +% 2 -> Sphercial (not yet implemented) +% +% example: +% CSX = InitCSX(); %for a default cartesian mesh +% or +% CSX = InitCSX('CoordSystem','1'); % for a cylindrical mesh definition +% +% See also DefineRectGrid, SmoothMeshLines, SmoothMeshLines2, +% SetBackgroundMaterial +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +CSX.Properties = []; + +%Cartesian mesh as default coordinate system +CSX.ATTRIBUTE.CoordSystem = 0; + +for n=1:2:numel(varargin) + CSX.ATTRIBUTE.(varargin{n}) = (varargin{n+1}); +end diff --git a/CSXCAD/matlab/RecursiveSmoothMesh.m b/CSXCAD/matlab/RecursiveSmoothMesh.m new file mode 100644 index 0000000..0ce791c --- /dev/null +++ b/CSXCAD/matlab/RecursiveSmoothMesh.m @@ -0,0 +1,176 @@ +function resultMeshPoints = RecursiveSmoothMesh(fixedMeshPoints, max_res, ratio, varargin)
+% function resultMeshPoints = RecursiveSmoothMesh(fixedMeshPoints, max_res, ratio, varargin)
+%
+% RecursiveSmoothMesh: This function tries to generate an optimal mesh
+% between the given fixedMeshPoints. The space between the fixed points is
+% filled with the largest possible mesh step size considdering mesh max
+% resolution and mesh increase/decrease ratio.
+%
+% Algorithm: Mesh creation is done via a simple try and error approach:
+% The next mesh point could be between lastDistance*ratio and
+% lastDistance/ratio away from the last mesh point. If a fixed mesh point
+% exist in this distance, this point is used. If there is a fixed mesh
+% point below this distance, the last (determined) mesh point is wrong, go
+% back and try with a mesh point a few mm below. (If that is not possible
+% because the distance to the mesh point behind the last point is to low,
+% also update this mesh point).
+% The algorithm seems to work very well, except if the ratio is to low.
+% At a value of 1.2, the calculation takes a very long time.
+%
+%
+% Parameter: fixedMeshPoints: List containing points at which a mesh
+% line should appear
+% max_res: Maximum distance between two mesh lines
+% ratio: Maximum grading ratio between two
+% neighbour mesh lines
+%
+% optional variable arguments ('key', value):
+% CheckMesh: Do a final mesh check (default is true)
+% allowed_max_ratio: allow only a given max. grading ratio
+% (default --> ratio*1.25)
+%
+% Returns: resultMeshPoints: List containing the positions of all
+% mesh lines. If a empty list is
+% returned, no resolution could be
+% found.
+%
+% CSXCAD matlab interface
+% -----------------------
+% Author: Florian Pfanner
+% Date: 28.09.2012
+
+
+if (nargin < 3)
+ ratio = 1.4;
+end
+
+if (ratio < 1.2)
+ warning('ratio must be > 1.2, set it to 1.2');
+ ratio = 1.2;
+end
+
+check_mesh = true;
+max_ratio = ratio*1.25;
+
+for vn=1:2:numel(varargin)
+ if (strcmpi(varargin{vn},'CheckMesh'))
+ check_mesh = varargin{vn+1};
+ end
+ if (strcmpi(varargin{vn},'allowed_max_ratio'))
+ max_ratio = varargin{vn+1};
+ end
+end
+
+fixedMeshPoints = sort(unique(fixedMeshPoints(:)'));
+if (length(fixedMeshPoints) < 2)
+ error('Not enoughts points to create mesh!');
+end
+
+
+% special behaviour to generate a symetric mesh (If the fixedMeshPoints are
+% given in a symetric way.)
+symMesh = CheckSymmtricLines(fixedMeshPoints);
+if (symMesh)
+ nPoints = length(fixedMeshPoints);
+
+ symMeshEven = mod(nPoints,2)==0;
+
+ if (symMeshEven)
+ center = mean(fixedMeshPoints);
+ if (center-fixedMeshPoints(nPoints/2) > 0.5*max_res)
+ fixedMeshPoints = [fixedMeshPoints(1:nPoints/2) center];
+ else
+ % this is not working!!!
+ % fixedMeshPoints = fixedMeshPoints(1:nPoints/2+1);
+ fixedMeshPoints = [fixedMeshPoints(1:nPoints/2) center];
+ end
+ else
+ % For symetric mesh, only the first half of the fixed mesh points is
+ % needed. But add center of the symetric mesh to the fixed mesh points
+ % to ensure the center is included.
+ fixedMeshPoints = fixedMeshPoints(1:(nPoints+1)/2);
+ end
+end
+
+
+minDistance = min(fixedMeshPoints(2:end)-fixedMeshPoints(1:end-1));
+
+% two cases:
+% * minDistance is smaller than max_res
+% -> try spaces for the first mesh between max_res and minDistance/2
+% * max_res is smaller than minDistance
+% -> try spaces for the first mesh between max_res and max_res/10
+if (minDistance < max_res)
+ trySpaces = linspace(max_res, minDistance/2, 10);
+else
+ trySpaces = linspace(max_res, max_res/10, 10);
+end
+for k=1:length(trySpaces)
+ [resultMeshPoints, success] = recursiveMeshSearch(fixedMeshPoints(1), ...
+ trySpaces(k), fixedMeshPoints(2:end), max_res, ratio);
+ if (success)
+ resultMeshPoints = [fixedMeshPoints(1) resultMeshPoints];
+ break;
+ end
+end
+
+
+if (symMesh)
+ resultMeshPoints = [resultMeshPoints(1:end-1) max(resultMeshPoints)*2-resultMeshPoints(end:-1:1)];
+end
+
+% check result
+if (check_mesh)
+ CheckMesh(resultMeshPoints,0,max_res,ratio,0);
+end
+
+end % main function
+
+function [meshPoints, success] = recursiveMeshSearch(lastMeshPoint, ...
+ lastSpace, fixedMeshPoints, max_res, ratio)
+
+nextMeshPointMax = lastMeshPoint + min(lastSpace*ratio, max_res);
+nextMeshPointMin = lastMeshPoint + lastSpace/ratio;
+
+% check if below ratio is a fixed mesh point -> abbort
+if (fixedMeshPoints(1) < nextMeshPointMin)
+ meshPoints = [];
+ success = false;
+ return;
+end
+
+% check if in range of ratio is a fixed mesh point:
+if (fixedMeshPoints(1) >= nextMeshPointMin && fixedMeshPoints(1) <= nextMeshPointMax)
+ % yes, use this fixed mesh point
+ nextMeshPoint = fixedMeshPoints(1);
+ if (length(fixedMeshPoints) > 1)
+ [meshPointsCall, success] = recursiveMeshSearch(nextMeshPoint, ...
+ nextMeshPoint-lastMeshPoint, fixedMeshPoints(2:end), ...
+ max_res, ratio);
+ if (success) % backtracking was successful, return result
+ meshPoints = [nextMeshPoint meshPointsCall];
+ else % no possible resulution found, last step has to be repeated!
+ meshPoints = [];
+ end
+ else % this was the last mesh point, finish!
+ meshPoints = nextMeshPoint;
+ success = true;
+ end
+ return;
+end
+
+% try to set next mesh point, begin with highest width
+trySpace = linspace(nextMeshPointMax, nextMeshPointMin, 10);
+for i=1:length(trySpace)
+ [meshPointsCall, success] = recursiveMeshSearch(trySpace(i), ...
+ trySpace(i)-lastMeshPoint, fixedMeshPoints, max_res, ...
+ ratio);
+ if (success) % backtracking was successful, return
+ meshPoints = [trySpace(i) meshPointsCall];
+ return;
+ end
+end
+% no sucessful points found, return false
+meshPoints = [];
+success = false;
+end
diff --git a/CSXCAD/matlab/SetBackgroundMaterial.m b/CSXCAD/matlab/SetBackgroundMaterial.m new file mode 100644 index 0000000..5aaa6df --- /dev/null +++ b/CSXCAD/matlab/SetBackgroundMaterial.m @@ -0,0 +1,23 @@ +function CSX = SetBackgroundMaterial(CSX, varargin) +% function CSX = SetBackgroundMaterial(CSX, varargin)) +% +% Set the background material properties +% +% variable arguments: +% 'Epsilon' : background rel. electric permittivity (default 1) +% 'Kappa' : background electric conductivity (default 0) +% 'Mue' : background rel. magnetic permeability (default 1) +% +% example: +% CSX = InitCSX(); +% CSX = SetBackgroundMaterial(CSX, 'Epsilon', 4) +% +% See also InitCSX +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (c) 2013 + +for n=1:2:numel(varargin) + CSX.BackgroundMaterial.ATTRIBUTE.(varargin{n}) = (varargin{n+1}); +end diff --git a/CSXCAD/matlab/SetExcitationWeight.m b/CSXCAD/matlab/SetExcitationWeight.m new file mode 100644 index 0000000..67a50ee --- /dev/null +++ b/CSXCAD/matlab/SetExcitationWeight.m @@ -0,0 +1,53 @@ +function CSX = SetExcitationWeight(CSX, name, weight) +% function CSX = SetExcitationWeight(CSX, name, weight) +% +% Define weighting functions for x-, y- and z-direction of excitation +% +% The functions can use the variables: +% x,y,z +% rho for the distance to z-axis +% r for the distance to origin +% a for alpha (as in cylindircal and spherical coord systems) +% t for theta (as in the spherical coord system +% +% all these variables are not weighted with the drawing unit defined by +% the grid +% +% example: +% start=[0 0 0]; +% stop=[width height 0]; +% CSX = AddExcitation(CSX,'excite',0,[1 1 0]); +% weight{1} = '2*cos(0.0031416*x)*sin(0.0062832*y)'; +% weight{2} = '1*sin(0.0031416*x)*cos(0.0062832*y)'; +% weight{3} = 0; +% CSX = SetExcitationWeight(CSX,'excite',weight); +% CSX = AddBox(CSX,'excite',0 ,start,stop); +% +% See also AddExcitation, InitCSX, DefineRectGrid +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +if ~isfield(CSX,'Properties') + error('CSXCAD::SetExcitationWeight: no properties not found'); +end +if ~isfield(CSX.Properties,'Excitation') + error('CSXCAD::SetExcitationWeight: no excitation properties found'); +end + +pos=0; +for n=1:numel(CSX.Properties.Excitation) + if strcmp(CSX.Properties.Excitation{n}.ATTRIBUTE.Name, name) + pos=n; + end +end + +if (pos==0) + error('CSXCAD::SetExcitationWeight: property not found'); + return; +end + +CSX.Properties.Excitation{pos}.Weight.ATTRIBUTE.X = weight{1}; +CSX.Properties.Excitation{pos}.Weight.ATTRIBUTE.Y = weight{2}; +CSX.Properties.Excitation{pos}.Weight.ATTRIBUTE.Z = weight{3}; diff --git a/CSXCAD/matlab/SetMaterialProperty.m b/CSXCAD/matlab/SetMaterialProperty.m new file mode 100644 index 0000000..463d22e --- /dev/null +++ b/CSXCAD/matlab/SetMaterialProperty.m @@ -0,0 +1,29 @@ +function CSX = SetMaterialProperty(CSX, name, varargin) +% function CSX = SetMaterialProperty(CSX, name, varargin) +% +% Use this function to define the material constants: +% 'Epsilon': relative electric permittivity: [Epsilon] = 1 +% 'Mue': relative magnetic permeability: [Mue} = 1 +% 'Kappa': electric conductivity: [Kappa] = S/m +% 'Sigma': magnetic conductivity (non-physical property): [Sigma] = Ohm/m +% 'Density': material mass density: [Density] = kg/m^3, e.g. water: 1000 +% necessary for SAR calculations +% +% examples: +% CSX = AddMaterial( CSX, 'RO3010' ); +% CSX = SetMaterialProperty( CSX, 'RO3010', 'Epsilon', 10.2, 'Mue', 1 ); +% CSX = AddBox( CSX, 'RO3010', 0, [0 0 0], [100 1000 1000] ); +% +% % anisotropic material +% CSX = AddMaterial( CSX, 'sheet','Isotropy',0); +% CSX = SetMaterialProperty(CSX, 'sheet', 'Kappa', [0 0 kappa]); +% CSX = AddBox( CSX, 'sheet', 0, [0 0 0], [10 1000 1000] ); +% +% See also AddMaterial, SetMaterialWeight +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +CSX = SetPropertyArgs(CSX, GetPropertyType(CSX,name), name, 'Property', varargin{:}); + diff --git a/CSXCAD/matlab/SetMaterialWeight.m b/CSXCAD/matlab/SetMaterialWeight.m new file mode 100644 index 0000000..825892e --- /dev/null +++ b/CSXCAD/matlab/SetMaterialWeight.m @@ -0,0 +1,31 @@ +function CSX = SetMaterialWeight(CSX, name, varargin) +% function CSX = SetMaterialWeight(CSX, name, varargin) +% +% Define the material weighting function(s) +% +% The functions can use the variables: +% x,y,z +% rho for the distance to z-axis +% r for the distance to origin +% a for alpha (as in cylindircal and spherical coord systems) +% t for theta (as in the spherical coord system +% +% all these variables are not weighted with the drawing unit defined by +% the grid +% +% example: +% %material distribution as a rect-function with 4 periods +% start=[-500 -100 -500]; +% stop =[ 500 100 500]; +% CSX = AddMaterial(CSX, 'material'); +% CSX = SetMaterialProperty(CSX, 'material', 'Epsilon', 1); +% CSX = SetMaterialWeight(CSX, 'material', 'Epsilon', ['(sin(4*z / 1000 *2*pi)>0)+1']); +% CSX = AddBox(CSX, 'material' ,10 , start, stop); +% +% See also AddMaterial, SetMaterialProperty, InitCSX, DefineRectGrid +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +CSX = SetPropertyArgs(CSX, GetPropertyType(CSX,name), name, 'Weight', varargin{:}); diff --git a/CSXCAD/matlab/SmoothMesh.m b/CSXCAD/matlab/SmoothMesh.m new file mode 100644 index 0000000..7ba37bc --- /dev/null +++ b/CSXCAD/matlab/SmoothMesh.m @@ -0,0 +1,73 @@ +function [mesh] = SmoothMesh( mesh, max_res, ratio, varargin) +% function [mesh] = SmoothMesh( mesh, max_res, <ratio, varargin>) +% +% Convienent function to create a smooth mesh in all directions. +% Generate smooth mesh by choosing an appropriate algorithm in each direction. +% +% Currently supported smoothing algorithm: +% SmoothMeshLines, SmoothMeshLines2 and RecursiveSmoothMesh +% +% arguments: +% lines: given fixed lines to create a smooth mesh in between +% max_res: scalar or vector of desired max. allowed resolution +% ratio: grading ratio: scalar or vector of desired neighboring +% line-delta ratio (optional, default is 1.5) +% - see also 'allowed_max_ratio' argument +% +% variable arguments ('keyword',value): +% algorithm: define subset of tried algorihm, e.g. [1 3] +% symmetric: 0/1 force symmetric mesh (default is input symmetry) +% homogeneous: 0/1 force homogeneous mesh +% allowed_min_res: allow a given min resolution only +% allowed_max_ratio: allow only a given max. grading ratio +% (default --> ratio*1.25) +% debug: 0/1 off/on +% +% example: +% mesh.x = [-BoundBox 0 BoundBox]; +% mesh.y = [-BoundBox 0 BoundBox]; +% mesh.z = [0 BoundBox]; +% mesh = SmoothMesh(mesh, lambda/20/unit); +% CSX = DefineRectGrid(CSX, unit, mesh); +% +% See also AutoSmoothMeshLines, InitCSX, DefineRectGrid, DetectEdges +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (C) 2012 + +if (nargin<3) + ratio = 1.5; +end + +if (numel(max_res)==1) + max_res = [max_res max_res max_res]; +end +if (numel(ratio)==1) + ratio = [ratio ratio ratio]; +end +if (numel(max_res)~=1) + max_res = [max_res max_res max_res]; +end + +if isfield(mesh,'x') + mesh.x = AutoSmoothMeshLines(mesh.x, max_res(1), ratio(1), varargin{:}); +elseif isfield(mesh,'r') + mesh.r = AutoSmoothMeshLines(mesh.r, max_res(1), ratio(1), varargin{:}); +else + error 'x/(r) direction not found' +end + +if isfield(mesh,'y') + mesh.y = AutoSmoothMeshLines(mesh.y, max_res(2), ratio(2), varargin{:}); +elseif isfield(mesh,'a') + mesh.a = AutoSmoothMeshLines(mesh.a, max_res(2), ratio(2), varargin{:}); +else + error 'y/(a) direction not found' +end + +if isfield(mesh,'z') + mesh.z = AutoSmoothMeshLines(mesh.z, max_res(3), ratio(3), varargin{:}); +else + error 'z direction not found' +end diff --git a/CSXCAD/matlab/SmoothMeshLines.m b/CSXCAD/matlab/SmoothMeshLines.m new file mode 100644 index 0000000..bc1160f --- /dev/null +++ b/CSXCAD/matlab/SmoothMeshLines.m @@ -0,0 +1,122 @@ +function lines = SmoothMeshLines( lines, max_res, ratio, varargin) +%function lines = SmoothMeshLines( lines, max_res, ratio, varargin) +% +% create smooth mesh lines +% +% Warning: This function may not always produce a desired output. +% +% lines: given fixed lines to create a smooth mesh in between +% max_res: desired max. resolution +% ratio: max. neighboring line-delta ratio, (optional, default is 1.3) +% +% optional variable arguments ('key', value) +% recursive: SmoothMeshLines a couple of times recursivly (default is 0) +% CheckMesh: Do a final mesh check (default is true) +% allowed_max_ratio: allow only a given max. grading ratio +% (default --> ratio*1.25) +% +% example: +% % create a x-mesh with lines at 0, 50 and 200 an a desired mesh +% resolution of 5 +% mesh.x = SmoothMeshLines([0 50 200],5,1.3); +% +% See also InitCSX, DefineRectGrid +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +if (numel(lines)<2) + return +end + +if (nargin<3) + ratio = 1.3; +end + +recursive = 0; +check_mesh = true; +max_ratio = ratio*1.25; + +for vn=1:2:numel(varargin) + if (strcmpi(varargin{vn},'recursive')) + recursive = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'CheckMesh')) + check_mesh = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'allowed_max_ratio')) + max_ratio = varargin{vn+1}; + end +end + +lines = unique(sort(lines)); + +diff_Lines = diff(lines); + +index = find(diff_Lines>(1.001*max_res)); + +% for n=1:numel(diff_Lines)-1 +% if ( (diff_Lines(n+1)/diff_Lines(n) > ratio) ) +% index = [index n+1]; +% end +% end + +index = unique(index); + +addLines = []; + +for n=1:numel(index) + if (index(n)==1) + start_res = max_res; + else + start_res = lines(index(n)) - lines(index(n)-1); + end + + if ((index(n)+1)==numel(lines)) + stop_res = max_res; + else + stop_res = lines(index(n)+2) - lines(index(n)+1); + end + + addLines = [addLines SmoothRange(lines(index(n)),lines(index(n)+1),start_res,stop_res,max_res,ratio)]; +end + + +lines = unique(sort([lines addLines])); + +addLines = []; +% relax ratio for test +ratio_relax = ratio + (ratio-1) * 1; + +[EC pos E_type] = CheckMesh(lines,0,max_res,ratio_relax,1); +diff_Lines = diff(lines); + +for n=1:EC + if pos(n)>1 + start_res = diff_Lines(pos(n)-1); + else + start_res = diff_Lines(pos(n)); + end + if pos(n) >= numel(diff_Lines) + stop_res = diff_Lines(end); + else + stop_res = diff_Lines(pos(n)+1); + end + max_res_R = max([start_res stop_res])/2/ratio; + addLines = [addLines SmoothRange(lines(pos(n)),lines(pos(n)+1),start_res,stop_res,max_res_R,ratio)]; +end + +lines = unique(sort([lines addLines])); + +for n=1:recursive + old_numL = numel(lines); + lines = SmoothMeshLines( lines, max_res, ratio, 0); + if numel(lines) == old_numL + return + end +end + +if (check_mesh) + CheckMesh(lines,0,max_res,max_ratio,0); +end diff --git a/CSXCAD/matlab/SmoothMeshLines2.m b/CSXCAD/matlab/SmoothMeshLines2.m new file mode 100644 index 0000000..5d347c0 --- /dev/null +++ b/CSXCAD/matlab/SmoothMeshLines2.m @@ -0,0 +1,318 @@ +function lines = SmoothMeshLines2(lines, max_res, ratio, varargin) +%lines = SmoothMeshLines2( lines, max_res [, ratio, varargin] ) +% +% Create smooth mesh lines. +% +% input: +% lines: 1xn vector of (fixed) mesh lines +% max_res: maximum distance between any two lines (e.g. lambda/10) +% ratio: (optional) default: 1.3 +% +% optional variable arguments ('key', value) +% CheckMesh: Do a final mesh check (default is true) +% allowed_max_ratio: allow only a given max. grading ratio +% (default --> ratio*1.25) +% output: +% lines: 1xn vector of (smoothed) mesh lines +% +% example: +% mesh.x = [0 100 2300 2400]; +% mesh.x = SmoothMeshLines2( mesh.x, 43 ); +% +% todo: +% - if gaps cannot be completely filled because of the ratio restriction, +% the spacing is not optimal. SmoothRange() needs to be optimized. +% - SmoothRange() has special handling if stop_taper is to long; can this +% also happen for start_taper? +% +% CSXCAD matlab interface +% ----------------------- +% (C) 2010 Sebastian Held <sebastian.held@uni-due.de> +% See also SmoothMeshLines + +if (numel(lines)<2) + return +end + +if (nargin<3) + ratio = 1.3; +end + +check_mesh = true; +max_ratio = ratio*1.25; + +for vn=1:2:numel(varargin) + if (strcmpi(varargin{vn},'CheckMesh')) + check_mesh = varargin{vn+1}; + end + if (strcmpi(varargin{vn},'allowed_max_ratio')) + max_ratio = varargin{vn+1}; + end +end + +% +% | | | | +% | gap(1) | gap(2) | ... | +% | | | | +% +% lines(1) lines(2) ... lines(end) + +lines = unique(sort(lines)); + +for n=1:(numel(lines)-1) + old_gap(n).start_res = -1; + old_gap(n).stop_res = -1; + old_gap(n).lines = []; +end + +if numel(lines) == 2 + % special case + addLines = SmoothRange( lines(1), lines(2), [], [], max_res, ratio ); + lines = sort(unique([lines addLines])); + [EC pos] = CheckMesh(lines,0,max_res,ratio); + return +end + +while 1 + gap = calc_gaps(lines,old_gap,max_res); + + % determine gap to process + index = find( [old_gap(2:end).start_res] ~= [gap(2:end).start_res], 1, 'first' ) + 1; + if isempty(index) + index = find( [old_gap(1:end-1).stop_res] ~= [gap(1:end-1).stop_res], 1, 'first' ); + end + if isempty(index) + break; % done + end + + start_res = min( max_res, calc_start_res(index,lines,old_gap) ); + stop_res = min( max_res, calc_stop_res(index,lines,old_gap) ); + + % create new lines + old_gap(index).lines = SmoothRange( lines(index), lines(index+1), start_res, stop_res, max_res, ratio ); + + % remember gap setting + old_gap(index).start_res = start_res; + old_gap(index).stop_res = stop_res; + + % debug +% plot_lines +end + +% merge lines +for n=1:numel(old_gap) + lines = [lines old_gap(n).lines]; +end +lines = sort(unique(lines)); + +if (check_mesh) + CheckMesh(lines,0,max_res,max_ratio,0); +end + +end % SmoothMeshLines2() + +%% ------------------------------------------------------------------------ +function gap = calc_gaps(lines,old_gap,max_res) + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + if numel(temp_lines) == 2 + gap(1).start_res = inf; % resolution not fixed + gap(1).stop_res = max_res; + return + end + + gap(1).start_res = inf; % resolution not fixed + idx = interp1( temp_lines, 1:numel(temp_lines), lines(2), 'nearest' ); + gap(1).stop_res = min( max_res, temp_lines(idx+1) - temp_lines(idx) ); + for n=2:numel(lines)-2 + idx = interp1( temp_lines, 1:numel(temp_lines), lines(n), 'nearest' ); + gap(n).start_res = min( max_res, temp_lines(idx) - temp_lines(idx-1) ); + idx = interp1( temp_lines, 1:numel(temp_lines), lines(n+1), 'nearest' ); + gap(n).stop_res = min( max_res, temp_lines(idx+1) - temp_lines(idx) ); + gap(n).lines = old_gap(n).lines; + end + idx = interp1( temp_lines, 1:numel(temp_lines), lines(end-1), 'nearest' ); + gap(numel(lines)-1).start_res = min( max_res, temp_lines(idx) - temp_lines(idx-1) ); + gap(numel(lines)-1).stop_res = inf; % resolution not fixed +end % calc_gaps() + + +%% ------------------------------------------------------------------------ +function plot_lines + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + plot( temp_lines, ones(size(temp_lines)), 'r+' ); + hold on + plot( lines, ones(size(lines)), 'bo' ); + hold off +end % plot_lines + + +%% ------------------------------------------------------------------------ +function res = calc_start_res(pos,lines,old_gap) + if (pos < 2) || (pos > numel(lines)) + res = []; + return + end + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + idx = interp1( temp_lines, 1:numel(temp_lines), lines(pos), 'nearest' ); + res = temp_lines(idx) - temp_lines(idx-1); +end % calc_res() + + +%% ------------------------------------------------------------------------ +function res = calc_stop_res(pos,lines,old_gap) + if (pos < 1) || (pos >= numel(lines)-1) + res = []; + return + end + temp_lines = lines; + for n=1:numel(old_gap) + temp_lines = [temp_lines old_gap(n).lines]; + end + temp_lines = sort(unique(temp_lines)); + + idx = interp1( temp_lines, 1:numel(temp_lines), lines(pos+1), 'nearest' ); + res = temp_lines(idx+1) - temp_lines(idx); +end % calc_res() + + +%% ------------------------------------------------------------------------ +function lines = SmoothRange(start, stop, start_res, stop_res, max_res, ratio) + +if (nargin<6) + ratio = 1.3; +end + +if isempty(start_res) && isempty(stop_res) + % special case: fill entire range with max_res + n1 = ceil( (stop-start) / max_res ) + 1; + lines = linspace( start, stop, n1 ); + return +end + +if isempty(start_res) + % special case: taper from stop_res at stop to max_res at start (if + % possible) + taper = stop_res*ratio; + stop_taper = stop; + while (taper*ratio < max_res) && (stop_taper(1)-taper > start) + stop_taper = [stop_taper(1)-taper stop_taper]; + taper = taper*ratio; + end + if (stop_taper(1) - start) >= max_res + % fill with equidistant lines + n1 = ceil( (stop_taper(1)-start) / max_res ) + 1; + stop_taper = [linspace(start,stop_taper(1),n1) stop_taper]; + else + % not enough space for entire taper + stop_taper(1) = []; % likely too near to start + if numel(stop_taper) > 0 + stop_taper = [(start+stop_taper(1))/2 stop_taper]; % create a centered line + end + end + lines = sort(unique(stop_taper(stop_taper>=start))); + return +end + +if isempty(stop_res) + % special case: taper from stop_res at stop to max_res at start (if + % possible) + taper = start_res*ratio; + start_taper = start; + while (taper*ratio < max_res) && (start_taper(end)+taper < stop) + start_taper = [start_taper start_taper(end)+taper]; + taper = taper*ratio; + end + if (stop - start_taper(end)) >= max_res + % fill with equidistant lines + n1 = ceil( (stop-start_taper(end)) / max_res ) + 1; + start_taper = [start_taper linspace(start_taper(end),stop,n1)]; + else + % not enough space for entire taper + start_taper(end) = []; % likely too near to stop + if numel(start_taper) > 0 + start_taper = [start_taper (stop+start_taper(end))/2]; % create a centered line + end + end + lines = sort(unique(start_taper(start_taper<=stop))); + return +end + +taper = start_res*ratio; +start_taper = start; +while (taper*ratio<max_res) && (start_taper(end)+taper < stop) + start_taper = [start_taper start_taper(end)+taper]; + taper = taper*ratio; +end +if (numel(start_taper) > 1) && (start_taper(end) - start_taper(end-1) > stop - start_taper(end)) + % not enough space for entire taper + start_taper(end) = []; % likely too near to stop + start_taper = [start_taper (stop+start_taper(end))/2]; % create a centered line +end + +taper = stop_res*ratio; +stop_taper = stop; +while (taper*ratio<max_res) && (stop_taper(1)-taper > start) + stop_taper = [stop_taper(1)-taper stop_taper]; + taper = taper*ratio; +end +if (numel(stop_taper) > 1) && (stop_taper(2) - stop_taper(1) > start - stop_taper(1)) + % not enough space for entire taper + stop_taper(1) = []; % likely too near to start + stop_taper = [(start+stop_taper(1))/2 stop_taper]; % create a centered line +end + +while ~isempty(start_taper) && ~isempty(stop_taper) && (stop_taper(1) < start_taper(end)) + + diff_start = diff(start_taper); + if isempty(diff_start) + diff_start = start_res; + end + diff_stop = diff(stop_taper); + if isempty(diff_stop) + diff_stop = stop_res; + end + + if (diff_start(end)>diff_stop(1)) + start_taper(end) = []; + else + stop_taper(1) = []; + end +end + +% it may happen, that the space between start_taper(end) and stop_taper(1) +% is very small; correct it: +if numel(start_taper)>=2 && numel(stop_taper)>=2 + d = diff( [start_taper(end-1) start_taper(end) stop_taper(1) stop_taper(2) ]); + if (d(1)/d(2) > ratio) || (d(3)/d(2) > ratio) + addLines = linspace( start_taper(end-1), stop_taper(2), 4 ); + lines = unique( [start_taper(1:end-1) addLines stop_taper(2:end)] ); + return + end +end + +if isempty(start_taper) + start_taper = start; +end +if isempty(stop_taper) + stop_taper = stop; +end + +% fill remaining space with equidistant lines +numL = ceil((stop_taper(1) - start_taper(end))/max_res)+1; +lines = unique([start_taper linspace(start_taper(end),stop_taper(1),numL) stop_taper]); +end % SmoothRange() diff --git a/CSXCAD/matlab/export_empire.m b/CSXCAD/matlab/export_empire.m new file mode 100644 index 0000000..104a0bf --- /dev/null +++ b/CSXCAD/matlab/export_empire.m @@ -0,0 +1,305 @@ +function export_empire( CSX, FDTD, filename, options ) +% export_empire( CSX, FDTD, filename, options ) +% +% Exports the geometry defined in CSX to filename as an Empire python file. +% +% CSX: CSX-object created by InitCSX() +% FDTD: FDTD-object created by InitFDTD() +% filename: export filename (e.g. '/tmp/export.py') +% options (optional): struct +% options.ignore : cell array with property names to ignore +% +% CSXCAD matlab interface +% ----------------------- +% author: Sebastian Held <sebastian.held@gmx.de> +% 17. Jun 2010: initial version +% 28. Sep 2010: rewritten using Empire python scripting +% See also InitCSX InitFDTD + +if nargin < 3 + options = []; +end + +fid = fopen( filename, 'w' ); + +% write header +fprintf( fid, ['# Empire python script\n\n' ... + '# start Empire and select File/"Run Script" and select this file\n\n'] ); +fprintf( fid, 'from scripting import *\n\n' ); +export_empire_meshlines( fid, CSX ); +export_empire_settings( fid, CSX, FDTD ); + +if isfield(CSX.Properties,'Material') + % process material + Material = CSX.Properties.Material; + options.Property = 'Material'; + process_primitives( fid, Material, options ); +end +if isfield(CSX.Properties,'Metal') + % process PEC + Metal = CSX.Properties.Metal; + options.Property = 'Metal'; + process_primitives( fid, Metal, options ); +end + +fclose( fid ); + +% ----------------------------------------------------------------------------- +function str = gym_vect( start,stop ) +str = ['[' num2str(start(1)) ',' num2str(stop(1)) ',' num2str(start(2)) ',' num2str(stop(2)) ',' num2str(start(3)) ',' num2str(stop(3)) ']']; +function str = gym_vector( v ) +str = ['[' num2str(v(1))]; +for a=2:numel(v), str = [str ', ' num2str(v(a))]; end +str = [str ']']; + +% ----------------------------------------------------------------------------- +function primitive_box( fid, CSX_box, layertype ) +%primitive_box( fid, CSX_box, layertype ) +% layertype may be: 'dielectric Dielectric' or 'conductor Conductor' +global current_layer_prio +properties = ''; +if current_layer_prio ~= CSX_box.ATTRIBUTE.Priority + properties = [properties ' prio ' num2str(CSX_box.ATTRIBUTE.Priority)]; +end +if ~isempty(properties) + fprintf( fid, 'g.objecttext("%s")\n', [layertype properties] ); +end +start = [CSX_box.P1.ATTRIBUTE.X CSX_box.P1.ATTRIBUTE.Y CSX_box.P1.ATTRIBUTE.Z]; +stop = [CSX_box.P2.ATTRIBUTE.X CSX_box.P2.ATTRIBUTE.Y CSX_box.P2.ATTRIBUTE.Z]; +fprintf( fid, 'g.layerextrude("z",%g,%g)\n', start(3), stop(3) ); +fprintf( fid, 'g.box(%g,%g,%g,%g)\n', start(1), start(2), stop(1), stop(2) ); + +% ----------------------------------------------------------------------------- +function primitive_cylinder( fid, CSX_cylinder, layertype ) +%primitive_cylinder( fid, CSX_cylinder, layertype ) +% layertype may be: 'dielectric Dielectric' or 'conductor Conductor' +global current_layer_prio +properties = ''; +if current_layer_prio ~= CSX_cylinder.ATTRIBUTE.Priority + properties = [properties ' prio ' num2str(CSX_cylinder.ATTRIBUTE.Priority)]; +end +if ~isempty(properties) + fprintf( fid, 'g.objecttext("%s")\n', [layertype properties] ); +end + +start = [CSX_cylinder.P1.ATTRIBUTE.X CSX_cylinder.P1.ATTRIBUTE.Y CSX_cylinder.P1.ATTRIBUTE.Z]; +stop = [CSX_cylinder.P2.ATTRIBUTE.X CSX_cylinder.P2.ATTRIBUTE.Y CSX_cylinder.P2.ATTRIBUTE.Z]; + +radius = CSX_cylinder.ATTRIBUTE.Radius; +if start([1 2]) == stop([1 2]) + fprintf( fid, 'g.layerextrude("z",%g,%g)\n', start(3), stop(3) ); + fprintf( fid, 'g.circle(%g,%g,%g)\n', start(1), start(2), radius ); +elseif start([1 3]) == stop([1 3]) + fprintf( fid, 'g.layerextrude("y",%g,%g)\n', start(2), stop(2) ); + fprintf( fid, 'g.circle(%g,%g,%g)\n', start(3), start(1), radius ); +elseif start([2 3]) == stop([2 3]) + fprintf( fid, 'g.layerextrude("x",%g,%g)\n', start(1), stop(1) ); + fprintf( fid, 'g.circle(%g,%g,%g)\n', start(2), start(3), radius ); +else + disp( ['cylinder coordinates: (' num2str(start) ') -> (' num2str(stop) ')'] ); + error 'cylinder orientation not supported' +end + +% ----------------------------------------------------------------------------- +function str = primitive_wire( CSX_wire, options ) +radius = CSX_wire.ATTRIBUTE.WireRadius; +str = ['sphere_sweep { linear_spline, ' num2str(numel(CSX_wire.Vertex))]; +for a=1:numel(CSX_wire.Vertex) + % iterate over all vertices + v = [CSX_wire.Vertex{a}.ATTRIBUTE.X CSX_wire.Vertex{a}.ATTRIBUTE.Y CSX_wire.Vertex{a}.ATTRIBUTE.Z]; + str = [str ', ' pov_vect(v) ', ' num2str(radius)]; +end +str = [str ' ' options '}']; + +% ----------------------------------------------------------------------------- +function primitive_polygon( fid, CSX_polygon, layertype ) +%primitive_polygon( fid, CSX_polygon, layertype ) +% layertype may be: 'dielectric Dielectric' or 'conductor Conductor' +global current_layer_prio +properties = ''; +if current_layer_prio ~= CSX_polygon.ATTRIBUTE.Priority + properties = [properties ' prio ' num2str(CSX_polygon.ATTRIBUTE.Priority)]; +end +if ~isempty(properties) + fprintf( fid, 'g.objecttext("%s")\n', [layertype properties] ); +end +Elevation = CSX_polygon.ATTRIBUTE.Elevation; +if (CSX_polygon.NormDir.ATTRIBUTE.X ~= 0) + NormDir = 'x'; +elseif (CSX_polygon.NormDir.ATTRIBUTE.Y ~= 0) + NormDir = 'y'; +elseif (CSX_polygon.NormDir.ATTRIBUTE.Z ~= 0) + NormDir = 'z'; +end +coords = []; +for a=1:numel(CSX_polygon.Vertex) + % iterate over all vertices + coords = [coords CSX_polygon.Vertex{a}.ATTRIBUTE.X1 CSX_polygon.Vertex{a}.ATTRIBUTE.X2]; +end +fprintf( fid, 'g.layerextrude("%s",%g,%g)\n', NormDir, Elevation, Elevation+1 ); +fprintf( fid, 'g.poly(' ); +for n=1:numel(coords) + if n > 1 + fprintf( fid, ', %g', coords(n) ); + else + fprintf( fid, '%g', coords(n) ); + end +end +fprintf( fid, ')\n' ); + + + + +% ----------------------------------------------------------------------------- +function process_primitives( fid, prop, options ) +global current_layer_prio +ignore = {}; +if isfield(options,'ignore'), ignore = options.ignore; end + +% iterate over all properties and extract the priority +prio = []; +for num=1:numel(prop) + if isfield(prop{num}.Primitives,'Box') + for a=1:numel(prop{num}.Primitives.Box) + % iterate over all boxes + prio(end+1) = prop{num}.Primitives.Box{a}.ATTRIBUTE.Priority; + end + end + if isfield(prop{num}.Primitives,'Cylinder') + for a=1:numel(prop{num}.Primitives.Cylinder) + % iterate over all cylinders + prio(end+1) = prop{num}.Primitives.Cylinder{a}.ATTRIBUTE.Priority; + end + end + if isfield(prop{num}.Primitives,'Polygon') + for a=1:numel(prop{num}.Primitives.Polygon) + % iterate over all polygons + prio(end+1) = prop{num}.Primitives.Polygon{a}.ATTRIBUTE.Priority; + end + end +end +% Now prio is a vector full of priorities. Extract the priority, which is +% most often used. +u_prios = unique(prio); +count = histc(prio,u_prios); +[~,idx] = max(count); +current_layer_prio = u_prios(idx); + +% now export all properties/primitives +for num=1:numel(prop) + Name = prop{num}.ATTRIBUTE.Name; + if any( strcmp(Name,ignore) ) + disp( ['omitting ' Name '...'] ); + continue + end + + disp( ['processing ' prop{num}.ATTRIBUTE.Name '...'] ); + + if strcmp(options.Property,'Metal') + layertype = 'conductor'; + properties = [layertype ' sigma infinite sides auto rough 0 prio ' num2str(current_layer_prio) ' automodel auto']; + else + layertype = 'dielectric'; + epsr = 1; + sigma = 0; + if isfield(prop{num},'PropertyX') + if isfield(prop{num}.PropertyX.ATTRIBUTE,'Epsilon'), epsr = prop{num}.PropertyX.ATTRIBUTE.Epsilon; end + if isfield(prop{num}.PropertyX.ATTRIBUTE,'Kappa'), sigma = prop{num}.PropertyX.ATTRIBUTE.Kappa; end + end + properties = [layertype ' epsr ' num2str(epsr) ' sigma ' num2str(sigma) ' prio ' num2str(current_layer_prio) ' tand 0 density 0']; + end + + export_empire_layer( fid, Name, properties ) + + if isfield(prop{num}.Primitives,'Box') + for a=1:numel(prop{num}.Primitives.Box) + % iterate over all boxes + Box = prop{num}.Primitives.Box{a}; + primitive_box( fid, Box, layertype ); + end + end + if isfield(prop{num}.Primitives,'Cylinder') + for a=1:numel(prop{num}.Primitives.Cylinder) + % iterate over all cylinders + Cylinder = prop{num}.Primitives.Cylinder{a}; + primitive_cylinder( fid, Cylinder, layertype ); + end + end + if isfield(prop{num}.Primitives,'Wire') + for a=1:numel(prop{num}.Primitives.Wire) + % iterate over all wires + Wire = prop{num}.Primitives.Wire{a}; +% str = primitive_wire( Wire, obj_modifier ); +% fprintf( fid, '%s\n', str ); + end + warning( 'CSXCAD:export_empire', 'primitive wire currently unsupported' ); + end + if isfield(prop{num}.Primitives,'Polygon') + for a=1:numel(prop{num}.Primitives.Polygon) + % iterate over all polygons + Polygon = prop{num}.Primitives.Polygon{a}; + primitive_polygon( fid, Polygon, layertype ); + end + end + + fprintf( fid, 'g.load()\n\n' ); +end + + + + +function export_empire_layer( fid, name, layertext ) +fprintf( fid, '# create a new layer\n' ); +fprintf( fid, 'g = gmf()\n' ); +fprintf( fid, 'g.layer("%s")\n', name ); +fprintf( fid, 'g.layerextrude("z",0,0)\n' ); +fprintf( fid, 'g.layertext("%s")\n', layertext ); + +function export_empire_meshlines(fid,CSX) +fprintf( fid, '# create the mesh\n' ); +fprintf( fid, 'empty_disc()\n' ); +fprintf( fid, 'g = gmf()\n' ); +fprintf( fid, 'g.raw_write("DISC X %g")\n', CSX.RectilinearGrid.XLines ); +fprintf( fid, 'g.raw_write("DISC Y %g")\n', CSX.RectilinearGrid.YLines ); +fprintf( fid, 'g.raw_write("DISC Z %g")\n', CSX.RectilinearGrid.ZLines ); +fprintf( fid, 'g.load()\n\n' ); + +function empire_BC = convert_BC( openEMS_BC ) +if ischar(openEMS_BC) + if strcmp(openEMS_BC,'PEC') + empire_BC = 'electric'; + elseif strcmp(openEMS_BC,'PMC') + empire_BC = 'magnetic'; + elseif strncmp(openEMS_BC,'MUR',3) || strncmp(openEMS_BC,'PML',3) + empire_BC = 'pml 6'; + else + empire_BC = 'UNKNOWN'; + end +else + if openEMS_BC == 0 + empire_BC = 'electric'; + elseif openEMS_BC == 1 + empire_BC = 'magnetic'; + elseif openEMS_BC == 2 || openEMS_BC == 3 + empire_BC = 'pml 6'; + else + empire_BC = 'UNKNOWN'; + end +end + +function export_empire_settings(fid,CSX,FDTD) +fprintf( fid, '# settings\n' ); +fprintf( fid, 'set_fs_parameter( ''dx'', %g )\n', CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit ); +fprintf( fid, 'set_fs_parameter( ''xmin'', ''%s'' )\n', convert_BC(FDTD.BoundaryCond.ATTRIBUTE.xmin) ); +fprintf( fid, 'set_fs_parameter( ''xmax'', ''%s'' )\n', convert_BC(FDTD.BoundaryCond.ATTRIBUTE.xmax) ); +fprintf( fid, 'set_fs_parameter( ''ymin'', ''%s'' )\n', convert_BC(FDTD.BoundaryCond.ATTRIBUTE.ymin) ); +fprintf( fid, 'set_fs_parameter( ''ymax'', ''%s'' )\n', convert_BC(FDTD.BoundaryCond.ATTRIBUTE.ymax) ); +fprintf( fid, 'set_fs_parameter( ''zmin'', ''%s'' )\n', convert_BC(FDTD.BoundaryCond.ATTRIBUTE.zmin) ); +fprintf( fid, 'set_fs_parameter( ''zmax'', ''%s'' )\n', convert_BC(FDTD.BoundaryCond.ATTRIBUTE.zmax) ); +fprintf( fid, 'set_fs_parameter( ''f_0'', %g )\n', FDTD.Excitation.ATTRIBUTE.f0 ); +fprintf( fid, 'set_fs_parameter( ''f_c'', %g )\n', FDTD.Excitation.ATTRIBUTE.fc ); +f_start = max(0,FDTD.Excitation.ATTRIBUTE.f0 - FDTD.Excitation.ATTRIBUTE.fc); +f_end = max(0,FDTD.Excitation.ATTRIBUTE.f0 + FDTD.Excitation.ATTRIBUTE.fc); +fprintf( fid, 'set_fs_parameter( ''F_START'', %g )\n', f_start ); +fprintf( fid, 'set_fs_parameter( ''F_END'', %g )\n', f_end ); +fprintf( fid, 'set_fs_parameter( ''F_MID'', %g )\n', (f_start+f_end)/2 ); diff --git a/CSXCAD/matlab/export_excellon.m b/CSXCAD/matlab/export_excellon.m new file mode 100644 index 0000000..8877fcd --- /dev/null +++ b/CSXCAD/matlab/export_excellon.m @@ -0,0 +1,144 @@ +function export_excellon( CSX, filename, options ) +% export_excellon( CSX, filename, options ) +% +% Exports the geometry defined in CSX to filename as an excellon drill file. +% Only cylinders will be considered for export. +% The xy-plane is exported. +% +% CSX: CSX-object created by InitCSX() +% filename: export filename (e.g. '/tmp/export.gbr') +% options (optional): struct +% .header: (string) add this to the header of the file +% comment line must have ';' as the very first character +% .exclude: (cell array) ignore these CSX-properties +% .include: (cell array) include only these CSX-properties +% if neither .exclude or .include is specified, process all properties +% +% See also InitCSX +% CSXCAD matlab interface +% ----------------------- +% author: Sebastian Held <sebastian.held@gmx.de> +% 26. Aug 2010: initial version + + +if nargin < 3 + options = []; +end + +fid = fopen( filename, 'w' ); + +% create header +header = []; +header = [header 'M48\n']; +header = [header '; EXCELLON file written by CSXCAD\n' ]; +header = [header '; ' datestr(now) '\n' ]; +header = [header 'METRIC,LZ\n']; +header = [header 'VER,1\n']; +header = [header 'FMAT,1\n']; + +if isfield(options,'header') + header = [header options.header]; +end + +% determine the drawing unit +options.drawingunit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit; + +if isfield(CSX.Properties,'Material') + % process material + Material = CSX.Properties.Material; + options.Property = 'Material'; + [tools,drill] = process_primitives( Material, options ); +end +if isfield(CSX.Properties,'Metal') + % process PEC + Metal = CSX.Properties.Metal; + options.Property = 'Metal'; + [tools_,drill_] = process_primitives( Metal, options ); + tools = unique( [tools tools_] ); + drill = [drill drill_]; +end + +% check the generated tools +[~,m] = unique( cellfun( @(x) x(1:3), tools, 'UniformOutput', 0 ) ); +if length(m) ~= numel(tools) + disp( 'the tool list is not unique:' ); + disp( tools ); + error( 'the tool list is not unique!' ); +end +a = unique( cellfun( @(x) x(4), tools, 'UniformOutput', 0 ) ); +if (length(a) > 1) || (a ~= 'C') + disp( 'the tool list has errors:' ); + disp( tools ); + error( 'the tool list has errors!' ); +end + +% convert cell array of strings to string +tools = char(cellfun( @(x)[x '\n'], tools )); +drill = char(cellfun( @(x)[x '\n'], drill )); + +% save the file +fprintf( fid, header ); +fprintf( fid, tools ); +fprintf( fid, '%s\n', 'M95' ); % end of program header +fprintf( fid, drill ); +fprintf( fid, '%s\n', 'M00' ); % end +fclose( fid ); + + +% ----------------------------------------------------------------------------- +function str = coord(v) +% str = coord(v) +% 2D-vector v takes coordinates in unit m +x = sprintf( '%.3f', v(1)*1e3 ); % mm +y = sprintf( '%.3f', v(2)*1e3 ); % mm +str = ['X' x 'Y' y]; + +% ----------------------------------------------------------------------------- +function [tool,drill] = primitive_cylinder( CSX_cylinder, drawingunit ) +start = [CSX_cylinder.P1.ATTRIBUTE.X CSX_cylinder.P1.ATTRIBUTE.Y] * drawingunit; +stop = [CSX_cylinder.P2.ATTRIBUTE.X CSX_cylinder.P2.ATTRIBUTE.Y] * drawingunit; +radius = CSX_cylinder.ATTRIBUTE.Radius * drawingunit; +if start(1:2) == stop(1:2) + % via in z-plane + dia_mm = radius * 2 * 1e3; % mm + tool = sprintf( 'T%02iC1.3f', round(dia_mm*10), dia_mm ); + drill = sprintf( 'T%02i\n%s', round(dia_mm*10), coord(start(1:2)) ); +else + disp( 'omitting primitive cylinder, because the projection onto the z-plane is not a circle' ); +end + + +% ----------------------------------------------------------------------------- +function [tools,drill] = process_primitives( prop, options ) +exclude = {}; +if isfield(options,'exclude'), exclude = options.exclude; end +for num=1:numel(prop) + Name = prop{num}.ATTRIBUTE.Name; + if any( strcmp(Name,exclude) ) + disp( ['omitting ' Name '...'] ); + continue + end + if isfield(options,'include') + include = options.include; + if ~any( strcmp(Name,include) ) + disp( ['omitting ' Name '...'] ); + continue + end + end + + tools = {}; + drill = {}; + + disp( ['processing ' prop{num}.ATTRIBUTE.Name '...'] ); + fprintf( fid, '%s\n', ['%LN' Name '*%'] ); + if isfield(prop{num}.Primitives,'Cylinder') + for a=1:numel(prop{num}.Primitives.Cylinder) + % iterate over all cylinders + Cylinder = prop{num}.Primitives.Cylinder{a}; + [tool,drill_] = primitive_cylinder( fid, Cylinder ); + tools = [tools tool]; + drill = [drill drill_]; + end + end +end + diff --git a/CSXCAD/matlab/export_gerber.m b/CSXCAD/matlab/export_gerber.m new file mode 100644 index 0000000..8a50f6f --- /dev/null +++ b/CSXCAD/matlab/export_gerber.m @@ -0,0 +1,198 @@ +function export_gerber( CSX, filename, options ) +% export_gerber( CSX, filename, options ) +% +% Exports the geometry defined in CSX to filename as a gerber RS-274X file. +% The xy-plane is exported. +% +% CSX: CSX-object created by InitCSX() +% filename: export filename (e.g. '/tmp/export.gbr') +% options (optional): struct +% .header: (string) add this to the header of the file +% .ignore: (cell array) ignore these CSX-properties +% +% See also InitCSX +% CSXCAD matlab interface +% ----------------------- +% author: Sebastian Held <sebastian.held@gmx.de> +% 13. Jun 2010: initial version + +%TODO: +% .include: (cell array) include only these CSX-properties +% if neither .ignore or .include is specified, process all properties + +if nargin < 3 + options = []; +end + +fid = fopen( filename, 'w' ); + +% write header +fprintf( fid, '%s\n', 'G04 gerber RS274X-file exported by openEMS*' ); +fprintf( fid, '%s\n', '%FSLAX66Y66*%' ); % leading zeros omitted; absolute coordinate values; 123456.123456 +fprintf( fid, '%s\n', '%MOMM*%' ); % set units to mm +fprintf( fid, '%s\n', '%INopenEMS export*%' ); % image name +fprintf( fid, '%s\n', '%ADD10C,0.00100*%' ); % aperture description +if isfield(options,'header') + fprintf( fid, '%s\n', header ); +end + +global drawingunit +drawingunit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit; +global CoordSystem +CoordSystem = CSX.ATTRIBUTE.CoordSystem; + +if isfield(CSX.Properties,'Material') + % process material + Material = CSX.Properties.Material; + options.Property = 'Material'; + process_primitives( fid, Material, options ); +end +if isfield(CSX.Properties,'Metal') + % process PEC + Metal = CSX.Properties.Metal; + options.Property = 'Metal'; + process_primitives( fid, Metal, options ); +end + +fprintf( fid, '%s\n', 'M02*' ); % end of program +fclose( fid ); + + +% ----------------------------------------------------------------------------- +function str = gerber_coord(v, CoordSystem) +global drawingunit + +if (CoordSystem==1) + r = v(1); + a = v(2); + v(1) = r*cos(a) * drawingunit; + v(2) = r*sin(a) * drawingunit; +else + v(1) = v(1) * drawingunit; + v(2) = v(2) * drawingunit; +end +x = sprintf( '%+013i', round(v(1)*1e9) ); % mm +y = sprintf( '%+013i', round(v(2)*1e9) ); % mm +if (numel(x) ~= 13) || (numel(y) ~= 13) + error( ['gerber_coord(): coordinate transformation failed: x=' x ' y=' y] ); +end +str = ['X' x 'Y' y]; + +% ----------------------------------------------------------------------------- +function primitive_box( fid, CSX_box, CoordSystem ) +start = [CSX_box.P1.ATTRIBUTE.X CSX_box.P1.ATTRIBUTE.Y]; +stop = [CSX_box.P2.ATTRIBUTE.X CSX_box.P2.ATTRIBUTE.Y]; + +if (CoordSystem==0) + fprintf( fid, '%s\n', 'G36*' ); + fprintf( fid, '%s\n', [gerber_coord(start, CoordSystem) 'D02*'] ); + fprintf( fid, '%s\n', [gerber_coord([stop(1) start(2)], CoordSystem) 'D01*'] ); + fprintf( fid, '%s\n', [gerber_coord([stop], CoordSystem) 'D01*'] ); + fprintf( fid, '%s\n', [gerber_coord([start(1) stop(2)], CoordSystem) 'D01*'] ); + fprintf( fid, '%s\n', [gerber_coord(start, CoordSystem) 'D01*'] ); + fprintf( fid, '%s\n', 'G37*' ); +elseif (CoordSystem==1) + r = sort([start(1) stop(1)]); + a = sort([start(2) stop(2)]); + max_arc = 0.5*pi/180; + a = linspace(a(1),a(2),ceil((a(2)-a(1))/max_arc)); + + fprintf( fid, '%s\n', 'G36*' ); + fprintf( fid, '%s\n', [gerber_coord([r(1) a(1)], CoordSystem) 'D02*'] ); + + for ang = a(2:end) + fprintf( fid, '%s\n', [gerber_coord([r(1) ang], CoordSystem) 'D01*'] ); + end + + for ang = fliplr(a) + fprintf( fid, '%s\n', [gerber_coord([r(2) ang], CoordSystem) 'D01*'] ); + end + + fprintf( fid, '%s\n', [gerber_coord([r(1) a(1)], CoordSystem) 'D01*'] ); + fprintf( fid, '%s\n', 'G37*' ); +else + error 'unknown coordinate system' +end + +% ----------------------------------------------------------------------------- +function primitive_cylinder( fid, CSX_cylinder, CoordSystem ) +global drawingunit +start = [CSX_cylinder.P1.ATTRIBUTE.X CSX_cylinder.P1.ATTRIBUTE.Y]; +stop = [CSX_cylinder.P2.ATTRIBUTE.X CSX_cylinder.P2.ATTRIBUTE.Y]; +radius = CSX_cylinder.ATTRIBUTE.Radius * drawingunit; +if start(1:2) == stop(1:2) + % via => draw circle + fprintf( fid, '%%ADD10C,%f*%%\n', radius*2*1e3 ); % aperture definition (mm) + fprintf( fid, '%s\n', 'G54D10*' ); % select aperture D10 + fprintf( fid, '%s\n', [gerber_coord(start, CoordSystem) 'D03*'] ); % flash +else + disp( ['omitting primitive cylinder, because the projection onto the z-plane is not a circle'] ); +end + + +% ----------------------------------------------------------------------------- +function primitive_polygon( fid, CSX_polygon, CoordSystem ) +if CSX_polygon.ATTRIBUTE.NormDir ~= 2 + disp( ['omitting primitive polygon, because the normal direction is not 2'] ); +end +fprintf( fid, '%s\n', 'G36*' ); +v = [CSX_polygon.Vertex{1}.ATTRIBUTE.X1 CSX_polygon.Vertex{1}.ATTRIBUTE.X2]; +fprintf( fid, '%s\n', [gerber_coord(v, CoordSystem) 'D02*'] ); +for a=2:numel(CSX_polygon.Vertex) + % iterate over all vertices + v = [CSX_polygon.Vertex{a}.ATTRIBUTE.X1 CSX_polygon.Vertex{a}.ATTRIBUTE.X2]; + fprintf( fid, '%s\n', [gerber_coord(v, CoordSystem) 'D01*'] ); +end +fprintf( fid, '%s\n', 'G37*' ); + +% ----------------------------------------------------------------------------- +function primCoordSystem = GetCoordSystem(prim) +global CoordSystem +primCoordSystem = CoordSystem; +if (isfield(prim.ATTRIBUTE,'CoordSystem')) + primCoordSystem = prim.ATTRIBUTE.CoordSystem; +end + +% ----------------------------------------------------------------------------- +function process_primitives( fid, prop, options ) +ignore = {}; +if isfield(options,'ignore'), ignore = options.ignore; end +for num=1:numel(prop) + Name = prop{num}.ATTRIBUTE.Name; + if any( strcmp(Name,ignore) ) + disp( ['omitting ' Name '...'] ); + continue + end + if strcmp(options.Property,'Material') && ~isfield(prop{num}.Property.ATTRIBUTE,'Kappa') + disp( ['omitting ' Name ' (no Kappa-value)...'] ); + continue + end + + disp( ['processing ' prop{num}.ATTRIBUTE.Name '...'] ); + fprintf( fid, '%s\n', ['%LN' Name '*%'] ); + if ~isfield(prop{num},'Primitives') + continue + end + if isfield(prop{num}.Primitives,'Box') + for a=1:numel(prop{num}.Primitives.Box) + % iterate over all boxes + Box = prop{num}.Primitives.Box{a}; + primitive_box( fid, Box, GetCoordSystem(Box) ); + end + end + if isfield(prop{num}.Primitives,'Cylinder') + for a=1:numel(prop{num}.Primitives.Cylinder) + % iterate over all cylinders + Cylinder = prop{num}.Primitives.Cylinder{a}; + primitive_cylinder( fid, Cylinder, GetCoordSystem(Cylinder) ); + end + end + if isfield(prop{num}.Primitives,'Polygon') + for a=1:numel(prop{num}.Primitives.Polygon) + % iterate over all polygons + Polygon = prop{num}.Primitives.Polygon{a}; + primitive_polygon( fid, Polygon, GetCoordSystem(Polygon) ); + end + end +end + diff --git a/CSXCAD/matlab/export_povray.m b/CSXCAD/matlab/export_povray.m new file mode 100644 index 0000000..a53cb13 --- /dev/null +++ b/CSXCAD/matlab/export_povray.m @@ -0,0 +1,237 @@ +function export_povray( CSX, filename, options ) +% export_povray( CSX, filename, options ) +% +% Exports the geometry defined in CSX to filename as a povray file. +% +% CSX: CSX-object created by InitCSX() +% filename: export filename (e.g. '/tmp/export.pov') +% options (optional): struct +% .camera: (string) use this as the camera definition line +% 1 camera on positive x-axis +% 2 camera on positive y-axis +% 3 camera on positive z-axis +% .light: (string) use this as the light definition line +% 1 point light at camera position +% .header: (string) add this to the header of the file +% .ignore: (cell array) ignore these CSX-properties +% .obj_modifier: struct +% .<propname>: (string) povray object modifier for corresponding primitives +% example: options.obj_modifier.copper = 'pigment { color rgbt <0.8,0.5,0,0> }'; +% +% See also InitCSX +% CSXCAD matlab interface +% ----------------------- +% author: Sebastian Held <sebastian.held@gmx.de> +% 12. Jun 2010: initial version + +if nargin < 3 + options = []; +end + +fid = fopen( filename, 'w' ); + +% write header +fprintf( fid, '%s\n', '// povray-file exported by openEMS' ); +fprintf( fid, '%s\n', '#include "colors.inc"' ); +fprintf( fid, '%s\n', '#include "metals.inc"' ); +fprintf( fid, '%s\n', '#include "textures.inc"' ); +fprintf( fid, '%s\n', '#include "transforms.inc"' ); +fprintf( fid, '%s\n', 'background { color rgb<1.000000,1.000000,1.000000> }' ); +if isfield(options,'header') + fprintf( fid, '%s\n', header ); +end + +if isfield(CSX.Properties,'Material') + % process material + Material = CSX.Properties.Material; + process_primitives( fid, Material, 'pigment { color rgbt <0.000, 0.533, 0.800,0.0> } finish { diffuse 0.6 }', options ); +end +if isfield(CSX.Properties,'Metal') + % process PEC + Metal = CSX.Properties.Metal; + process_primitives( fid, Metal, 'texture { Copper_Metal }', options ); +end + +% create coordinate system vectors (for debugging) +debug_coords = 0; +if debug_coords + fprintf( fid, [... + '#macro Axis_( AxisLen, RedTexture,WhiteTexture) \n' ... + 'union{\n' ... + ' cylinder {<0,-AxisLen,0>,<0,AxisLen,0>,0.5\n' ... + ' texture{checker texture{RedTexture} \n' ... + ' texture{WhiteTexture}\n' ... + ' translate<0.1,0,0.1>}}\n' ... + ' cone{<0,AxisLen,0>,2,<0,AxisLen+0.7,0>,0\n' ... + ' texture{RedTexture}}\n' ... + ' } // end of union \n' ... + '#end // of macro "Axis( )"\n' ... + '\n' ... + '#macro AxisXYZ( AxisLenX, AxisLenY, AxisLenZ, TexRed, TexWhite)\n' ... + 'union{\n' ... + ' object{Axis_(AxisLenX, TexRed, TexWhite) rotate< 0,0,-90>} // x axis\n' ... + ' object{Axis_(AxisLenY, TexRed, TexWhite) rotate< 0,0, 0>} // y axis \n' ... + ' object{Axis_(AxisLenZ, TexRed, TexWhite) rotate<90,0, 0>} // z axis\n' ... + ' text{ ttf"timrom.ttf", "x", 0.15, 0 texture{TexRed} \n' ... + ' scale 10 translate <AxisLenX+0.05,0.4,-0.10>}\n' ... + ' text{ ttf"timrom.ttf", "y", 0.15, 0 texture{TexRed} \n' ... + ' scale 10 translate <-0.75,AxisLenY+0.50,-0.00>}\n' ... + ' text{ ttf"timrom.ttf", "z", 0.15, 0 texture{TexRed} \n' ... + ' scale 10 translate <-0.75,0.2,AxisLenZ+0.50>}\n' ... + ' scale <1000,1000,1000>\n' ... + '} // end of union\n' ... + '#end// of macro\n' ... + '\n' ... + 'object{AxisXYZ( 10, 10, 10, texture{ pigment{rgb<1,0,0>} finish{ phong 1}}, texture{ pigment{rgb<1,1,1>} finish{ phong 1}} )}\n' ... + ''] ); +end + +% create camera and light +xmin = min(CSX.RectilinearGrid.XLines); +xmax = max(CSX.RectilinearGrid.XLines); +ymin = min(CSX.RectilinearGrid.YLines); +ymax = max(CSX.RectilinearGrid.YLines); +zmin = min(CSX.RectilinearGrid.ZLines); +zmax = max(CSX.RectilinearGrid.ZLines); +center = [(xmin+xmax)/2 (ymin+ymax)/2 (zmin+zmax)/2]; +% default camera +camera_pos = [center(1), ymin, zmax+sqrt((xmax-xmin)^2+(ymax-ymin)^2)]; +camera = ['camera {location ' pov_vect(camera_pos) ' right <-1.33,0,0> look_at ' pov_vect(center) ' angle 40}']; % right-handed coordinate system +if isfield(options,'camera') + if ischar(options.camera) + camera = options.camera; + elseif options.camera == 1 + % looking from positive x-axis + camera_pos = [xmax+sqrt((ymax-ymin)^2+(zmax-zmin)^2), 0, 0]; + camera = ['camera {location ' pov_vect(camera_pos) ' right <-1.33,0,0> look_at ' pov_vect(center) ' angle 40}']; + elseif options.camera == 2 + % looking from positive y-axis + camera_pos = [0, ymax+sqrt((xmax-xmin)^2+(zmax-zmin)^2), 0]; + camera = ['camera {location ' pov_vect(camera_pos) ' right <-1.33,0,0> look_at ' pov_vect(center) ' angle 40}']; + elseif options.camera == 3 + % looking from positive z-axis + camera_pos = [0, 0, zmax+sqrt((xmax-xmin)^2+(ymax-ymin)^2)]; + camera = ['camera {location ' pov_vect(camera_pos) ' right <-1.33,0,0> look_at ' pov_vect(center) ' angle 40}']; + end +end +fprintf( fid, '%s\n', camera ); +% default light +light = 'light_source { <3500,-3500,10000> White area_light <10000, 0, 0>, <0, 10000, 0>, 2, 2 adaptive 1 }'; +if isfield(options,'light') + if ischar(options.light) + light = options.light; + elseif options.light == 1 + % point light at position of camera + light = ['light_source { ' pov_vect(camera_pos) ', rgb <1,1,1> }']; + end +end +fprintf( fid, '%s\n', light ); + +fclose( fid ); + +% ----------------------------------------------------------------------------- +function str = pov_vect( vec ) +if numel(vec) == 3 + str = ['<' num2str(vec(1)) ',' num2str(vec(2)) ',' num2str(vec(3)) '>']; +else + str = ['<' num2str(vec(1)) ',' num2str(vec(2)) '>']; +end + +% ----------------------------------------------------------------------------- +function str = primitive_box( CSX_box, options ) +start = [CSX_box.P1.ATTRIBUTE.X CSX_box.P1.ATTRIBUTE.Y CSX_box.P1.ATTRIBUTE.Z]; +stop = [CSX_box.P2.ATTRIBUTE.X CSX_box.P2.ATTRIBUTE.Y CSX_box.P2.ATTRIBUTE.Z]; +if any( start == stop ) + % 2D box + % povray supports 2D polygons, but has no priority concept, therefore use the box primitive + epsilon = 1; % FIXME this should be small compared to any other linear dimension of any object in the scene + epsilon = (start == stop) * epsilon; % select identical components + start = start - epsilon; + stop = stop + epsilon; +end +str = ['box { ' pov_vect(start) ', ' pov_vect(stop) ' ' options '}']; + +% ----------------------------------------------------------------------------- +function str = primitive_cylinder( CSX_cylinder, options ) +start = [CSX_cylinder.P1.ATTRIBUTE.X CSX_cylinder.P0.ATTRIBUTE.Y CSX_cylinder.P1.ATTRIBUTE.Z]; +stop = [CSX_cylinder.P2.ATTRIBUTE.X CSX_cylinder.P1.ATTRIBUTE.Y CSX_cylinder.P2.ATTRIBUTE.Z]; +radius = CSX_cylinder.ATTRIBUTE.Radius; +str = ['cylinder { ' pov_vect(start) ', ' pov_vect(stop) ', ' num2str(radius) ' ' options '}']; + +% ----------------------------------------------------------------------------- +function str = primitive_wire( CSX_wire, options ) +radius = CSX_wire.ATTRIBUTE.WireRadius; +str = ['sphere_sweep { linear_spline, ' num2str(numel(CSX_wire.Vertex))]; +for a=1:numel(CSX_wire.Vertex) + % iterate over all vertices + v = [CSX_wire.Vertex{a}.ATTRIBUTE.X CSX_wire.Vertex{a}.ATTRIBUTE.Y CSX_wire.Vertex{a}.ATTRIBUTE.Z]; + str = [str ', ' pov_vect(v) ', ' num2str(radius)]; +end +str = [str ' ' options '}']; + +% ----------------------------------------------------------------------------- +function str = primitive_polygon( CSX_polygon, options ) +Elevation = -CSX_polygon.ATTRIBUTE.Elevation; +% NormDir = CSX_polygon.ATTRIBUTE.NormDir; +epsilon = 1; % FIXME this should be small compared to any other linear dimension of any object in the scene +str = ['prism { linear_spline linear_sweep ' num2str(Elevation - epsilon) ', ' num2str(Elevation + epsilon) ', ' num2str(numel(CSX_polygon.Vertex)+1)]; +for a=1:numel(CSX_polygon.Vertex) + % iterate over all vertices + v = [CSX_polygon.Vertex{a}.ATTRIBUTE.X1 CSX_polygon.Vertex{a}.ATTRIBUTE.X2]; + str = [str ', ' pov_vect(v)]; +end +v = [CSX_polygon.Vertex{1}.ATTRIBUTE.X1 CSX_polygon.Vertex{1}.ATTRIBUTE.X2]; % close prism +str = [str ', ' pov_vect(v)]; +% str = [str ' ' options ' Point_At_Trans(' pov_vect(NormDir) ')}']; % needs transforms.inc +str = [str ' ' options ' rotate<-90,0,0> }']; + +% ----------------------------------------------------------------------------- +function process_primitives( fid, prop, default_obj_modifier, options ) +ignore = {}; +if isfield(options,'ignore'), ignore = options.ignore; end +for num=1:numel(prop) + Name = prop{num}.ATTRIBUTE.Name; + if any( strcmp(Name,ignore) ) + disp( ['omitting ' Name '...'] ); + continue + end + obj_modifier = default_obj_modifier; + if isfield(options,'obj_modifier') && isfield(options.obj_modifier,Name) + obj_modifier = options.obj_modifier.(Name); + end + disp( ['processing ' prop{num}.ATTRIBUTE.Name '...'] ); + fprintf( fid, '%s\n', ['// ' Name] ); + if isfield(prop{num}.Primitives,'Box') + for a=1:numel(prop{num}.Primitives.Box) + % iterate over all boxes + Box = prop{num}.Primitives.Box{a}; + str = primitive_box( Box, obj_modifier ); + fprintf( fid, '%s\n', str ); + end + end + if isfield(prop{num}.Primitives,'Cylinder') + for a=1:numel(prop{num}.Primitives.Cylinder) + % iterate over all cylinders + Cylinder = prop{num}.Primitives.Cylinder{a}; + str = primitive_cylinder( Cylinder, obj_modifier ); + fprintf( fid, '%s\n', str ); + end + end + if isfield(prop{num}.Primitives,'Wire') + for a=1:numel(prop{num}.Primitives.Wire) + % iterate over all wires + Wire = prop{num}.Primitives.Wire{a}; + str = primitive_wire( Wire, obj_modifier ); + fprintf( fid, '%s\n', str ); + end + end + if isfield(prop{num}.Primitives,'Polygon') + for a=1:numel(prop{num}.Primitives.Polygon) + % iterate over all polygons + Polygon = prop{num}.Primitives.Polygon{a}; + str = primitive_polygon( Polygon, obj_modifier ); + fprintf( fid, '%s\n', str ); + end + end +end + diff --git a/CSXCAD/matlab/isOctave.m b/CSXCAD/matlab/isOctave.m new file mode 100644 index 0000000..b536e48 --- /dev/null +++ b/CSXCAD/matlab/isOctave.m @@ -0,0 +1,46 @@ +function [isOct,ver] = isOctave() +% [isOct,version] = isOctave() +% +% Function to test if matlab or octave is used. +% +% Output: +% isOct: bool; true, if interpreter is Octave +% version: struct +% .major: number; major version +% .minor: number; minor version +% .release: number; release version +% +% example (Octave 3.6.1): +% isOct = true +% version.major = 3 +% version.minor = 6 +% version.release = 1 +% +% example (Matlab 7.8.0.347 (R2009a)): +% isOct = false +% version.major = 7 +% version.minor = 8 +% version.release = 0 +% +% openEMS matlab/octave interface +% ----------------------- +% (C) 2011 Thorsten Liebig <thorsten.leibig@gmx.de> +% (C) 2012 Sebastian Held <sebastian.held@gmx.de> + +isOct = exist('OCTAVE_VERSION','builtin') ~= 0; + +if (isOct) + ver_cell = strsplit( OCTAVE_VERSION, '.' ); +else + remain = version(); + ver_cell = {}; + while ~isempty(remain) + [str, remain] = strtok(remain, '.'); + ver_cell{end+1} = str; + end +end + +ver = []; +ver.major = str2double(ver_cell{1}); +ver.minor = str2double(ver_cell{2}); +ver.release = str2double(ver_cell{3}); diff --git a/CSXCAD/matlab/private/Add2Property.m b/CSXCAD/matlab/private/Add2Property.m new file mode 100644 index 0000000..e38ddc9 --- /dev/null +++ b/CSXCAD/matlab/private/Add2Property.m @@ -0,0 +1,26 @@ +function CSX = Add2Property(CSX, propName, newPrim, primName) +% function CSX = Add2Property(CSX, propName, newPrim, primName) +% +% meant for internal use!! +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +type = GetPropertyType(CSX, propName); +if isempty(type) + error('CSXCAD:Add2Property',['the type for the property "' propName '" cannot be found']); +end + +pos = GetPropertyPosition(CSX, type, propName); +if (pos==0) + error('CSXCAD:Add2Property',['property "' propName '" of type "' type '" not found!']); +end + +if ~isfield(CSX.Properties.(type){pos}, 'Primitives') + CSX.Properties.(type){pos}.Primitives.(primName){1}=newPrim; +elseif ~isfield(CSX.Properties.(type){pos}.Primitives, primName) + CSX.Properties.(type){pos}.Primitives.(primName){1}=newPrim; +else + CSX.Properties.(type){pos}.Primitives.(primName){end+1}=newPrim; +end diff --git a/CSXCAD/matlab/private/AddPrimitiveArgs.m b/CSXCAD/matlab/private/AddPrimitiveArgs.m new file mode 100644 index 0000000..df3cf28 --- /dev/null +++ b/CSXCAD/matlab/private/AddPrimitiveArgs.m @@ -0,0 +1,28 @@ +function primitive = AddPrimitiveArgs(primitive, varargin) +% prim = AddPrimitiveArgs(primitive, varargin) +% +% meant for internal use only +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +transform = []; + +for n=1:2:(nargin-2) + if (strcmp(varargin{n},'Transform')) + transform = varargin([n:n+1]); + varargin([n:n+1]) = []; + break + end +end + +for n=1:2:numel(varargin) + primitive.ATTRIBUTE.(varargin{n}) = varargin{n+1}; +end + +if ~isempty(transform) + for n=1:2:numel(transform{2}) + primitive.Transformation.(transform{2}{n}).ATTRIBUTE.Argument=transform{2}{n+1}; + end +end diff --git a/CSXCAD/matlab/private/AddProperty.m b/CSXCAD/matlab/private/AddProperty.m new file mode 100644 index 0000000..ec4c2ed --- /dev/null +++ b/CSXCAD/matlab/private/AddProperty.m @@ -0,0 +1,29 @@ +function [CSX pos] = AddProperty(CSX, type, name, varargin) +% function [CSX pos] = AddProperty(CSX, type, name, varargin) +% +% internal function to add a property to CSX.Properties +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +% check if this property already exists +[type_found pos] = FindProperty(CSX, name); +if (pos>0) + error('CSXCAD:AddProperty',['property with name "' name '" with type "' type_found '" already exists! Choose a different name!']); +end + +if isfield(CSX.Properties,type) + pos = numel(CSX.Properties.(type))+1; +else + CSX.Properties.(type) = {}; % create cell array + pos = 1; +end + +CSX.Properties.(type){pos}.ATTRIBUTE.Name=name; +for n=1:numel(varargin)/2 + if ~ischar(varargin{2*n-1}) + error(['CSXCAD::AddProperty: not an attribute: ' varargin{2*n-1}]); + end + CSX.Properties.(type){pos}.ATTRIBUTE.(varargin{2*n-1})=varargin{2*n}; +end diff --git a/CSXCAD/matlab/private/CheckSymmtricLines.m b/CSXCAD/matlab/private/CheckSymmtricLines.m new file mode 100644 index 0000000..6106d63 --- /dev/null +++ b/CSXCAD/matlab/private/CheckSymmtricLines.m @@ -0,0 +1,34 @@ +function result = CheckSymmtricLines(lines) +% function result = CheckSymmtricLines(lines) +% +% check mesh lines for symmetry +% +% Note: make sure lines are sorted and unique +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (C) 2012 + +result = false; + +tolerance = 1e-10; +NP = numel(lines); +range = lines(end)-lines(1); +center = 0.5*(lines(end)+lines(1)); + +% check all lines for symmetry +for n=1:NP/2 + if (abs((center-lines(n))-(lines(end-n+1)-center)) > range*tolerance/NP) + return; + end +end + +% check central point to be symmetry-center +if (mod(NP,2)) + if (abs(lines((NP+1)/2)-center) > range*tolerance/NP) + return; + end +end + +% if all checks pass, return true +result = true; diff --git a/CSXCAD/matlab/private/FindProperty.m b/CSXCAD/matlab/private/FindProperty.m new file mode 100644 index 0000000..1df5ef6 --- /dev/null +++ b/CSXCAD/matlab/private/FindProperty.m @@ -0,0 +1,17 @@ +function [type pos] = FindProperty(CSX, name) +% function [type pos] = FindProperty(CSX, name) +% +% internal function to find a given property +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (c) 2013 + +pos = 0; +type = GetPropertyType(CSX, name); + +if isempty(type) + return; +end + +pos = GetPropertyPosition(CSX, type, name); diff --git a/CSXCAD/matlab/private/GetPropertyPosition.m b/CSXCAD/matlab/private/GetPropertyPosition.m new file mode 100644 index 0000000..8305876 --- /dev/null +++ b/CSXCAD/matlab/private/GetPropertyPosition.m @@ -0,0 +1,41 @@ +function pos = GetPropertyPosition(CSX, type, name) +% function pos = GetPropertyPosition(CSX, type, name) +% +% - internal function to get the position of property with name: <name> +% inside a given type +% - function will perform a series of validitiy tests +% - will return 0 if not found +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (c) 2013 + +pos = 0; + +if ~ischar(name) + error('CSXCAD::GetPropertyPosition: name must be a string'); +end + +if ~ischar(type) + error('CSXCAD::GetPropertyPosition: type name must be a string'); +end + +if ~isfield(CSX,'Properties') + error('CSXCAD:GetPropertyPosition','CSX.Properties is not defined'); +end + +if isempty(type) + error('CSXCAD:GetPropertyPosition','type is empty, maybe the property you requested is undefined'); +end + +% type not (yet) defined, thus <name> not found +if ~isfield(CSX.Properties,type) + return +end + +for n=1:numel(CSX.Properties.(type)) + if strcmp(CSX.Properties.(type){n}.ATTRIBUTE.Name, name) + pos=n; + return + end +end diff --git a/CSXCAD/matlab/private/GetPropertyType.m b/CSXCAD/matlab/private/GetPropertyType.m new file mode 100644 index 0000000..5a63dff --- /dev/null +++ b/CSXCAD/matlab/private/GetPropertyType.m @@ -0,0 +1,30 @@ +function type_name = GetPropertyType(CSX, name) +% function type_name = GetPropertyType(CSX, name) +% +% internal function to get the type of a given property +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig (c) 2010-2013 + +if ~ischar(name) + error('CSXCAD::GetPropertyType: name must be a string'); +end +if ~isfield(CSX,'Properties') + error('CSXCAD:GetPropertyPosition','CSX.Properties is not defined'); +end + +type_name = ''; +if isempty(CSX.Properties) + return +end + +prop_types = fieldnames(CSX.Properties); +for n=1:numel(prop_types) + for p = 1:numel(CSX.Properties.(prop_types{n})) + if (strcmp(CSX.Properties.(prop_types{n}){p}.ATTRIBUTE.Name,name)) + type_name = prop_types{n}; + return; + end + end +end diff --git a/CSXCAD/matlab/private/SetPropertyArgs.m b/CSXCAD/matlab/private/SetPropertyArgs.m new file mode 100644 index 0000000..893a642 --- /dev/null +++ b/CSXCAD/matlab/private/SetPropertyArgs.m @@ -0,0 +1,16 @@ +function CSX = SetPropertyArgs(CSX, type, name, property, varargin) +% CSX = SetPropertyArgs(CSX, type, name, property, varargin) +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +pos = GetPropertyPosition(CSX, type, name); + +if (pos==0) + error('CSXCAD:SetPropertyArgs',['property "' name '" of type "' type '" not found!']); +end + +for n=1:numel(varargin)/2 + CSX.Properties.(type){pos}.(property).ATTRIBUTE.(varargin{2*n-1}) = varargin{2*n}; +end diff --git a/CSXCAD/matlab/private/SmoothRange.m b/CSXCAD/matlab/private/SmoothRange.m new file mode 100644 index 0000000..d36db70 --- /dev/null +++ b/CSXCAD/matlab/private/SmoothRange.m @@ -0,0 +1,58 @@ +function [ lines ] = SmoothRange(start, stop, start_res, stop_res, max_res, ratio) +%function [ lines ] = SmoothRange(start, stop, start_res, stop_res, max_res, ratio) +% +% internal function only, use SmoothMeshLines instead +% +% See also SmoothMeshLines +% +% CSXCAD matlab interface +% ----------------------- +% author: Thorsten Liebig + +if (nargin<6) + ratio = 1.3; +end + +taper = start_res*ratio; +start_taper = start; +while (taper<max_res) + start_taper = [start_taper start_taper(end)+taper]; + taper = taper*ratio; +end + +taper = stop_res*ratio; +stop_taper = stop; +while (taper<max_res) + stop_taper = [stop_taper stop_taper(end)-taper]; + taper = taper*ratio; +end +stop_taper = sort(stop_taper); + +while ( (abs(stop_taper(1) - start_taper(end)) < max_res) || (stop_taper(1) < start_taper(end)) ) + + diff_start = diff(start_taper); + if isempty(diff_start) + diff_start = 0; + end + diff_stop = diff([stop_taper]); + if isempty(diff_stop) + diff_stop = 0; + end + + if (diff_start(end)>diff_stop(1)) + start_taper = start_taper(1:end-1); + else + stop_taper = stop_taper(2:end); + end + + if (numel(stop_taper)==0) || (numel(start_taper)==0) + break + end +end + +if (numel(stop_taper)==0) || (numel(start_taper)==0) + lines = unique([start_taper stop_taper]); +else + numL = ceil((stop_taper(1) - start_taper(end))/max_res)+1; + lines = unique([start_taper linspace(start_taper(end),stop_taper(1),numL) stop_taper]); +end
\ No newline at end of file diff --git a/CSXCAD/matlab/private/octave_struct2xml_2.m b/CSXCAD/matlab/private/octave_struct2xml_2.m new file mode 100644 index 0000000..3ebc549 --- /dev/null +++ b/CSXCAD/matlab/private/octave_struct2xml_2.m @@ -0,0 +1,57 @@ +function out = octave_struct2xml_2( in, rootName, indent ) + +out = [indent '<' rootName]; + +float_accuracy = 15; + +% process attributes +if isstruct( in ) + fnames = fieldnames( in ); + for n=1:numel(fnames) + current_field = fnames{n}; + if strcmp( current_field, 'ATTRIBUTE' ) + attributes = fieldnames( in.ATTRIBUTE ); + for m=1:numel( attributes ) + temp = in.ATTRIBUTE.(attributes{m}); + if ~ischar( temp ) + temp = vector2str( temp, float_accuracy ); + end + out = [out ' ' attributes{m} '="' temp '"']; + end + break + end + end +end + +out = [out '>']; + + +% process content +if ~isstruct( in ) + if ~ischar( in ) + temp = vector2str(in, float_accuracy); + else + temp = in; + end + out = [out temp '</' rootName '>\n']; + return +end + +out = [out '\n']; + +fnames = fieldnames( in ); +for n=1:numel(fnames) + current_field = fnames{n}; + if strcmp( current_field, 'ATTRIBUTE' ) + continue + end + if iscell( in.(current_field) ) + for m=1:numel( in.(current_field) ) + out = [out octave_struct2xml_2( in.(current_field){m}, current_field, [indent ' '] )]; + end + else + out = [out octave_struct2xml_2( in.(current_field), current_field, [indent ' '] )]; + end +end + +out = [out indent '</' rootName '>\n']; diff --git a/CSXCAD/matlab/private/struct_2_xmlNode.m b/CSXCAD/matlab/private/struct_2_xmlNode.m new file mode 100644 index 0000000..207632f --- /dev/null +++ b/CSXCAD/matlab/private/struct_2_xmlNode.m @@ -0,0 +1,48 @@ +function docElem = struct_2_xmlNode(docNode, docElem, mat_struct) + +if (isfield(mat_struct,'ATTRIBUTE')) + names = fieldnames(mat_struct.ATTRIBUTE); + for n=1:numel(names) + if isnumeric(getfield(mat_struct.ATTRIBUTE,names{n})) || islogical(getfield(mat_struct.ATTRIBUTE,names{n})) + docElem.setAttribute(names{n},vector2str(getfield(mat_struct.ATTRIBUTE,names{n}),15)); + else + docElem.setAttribute(names{n},getfield(mat_struct.ATTRIBUTE,names{n})); + end + end + mat_struct = rmfield(mat_struct,'ATTRIBUTE'); +end + +names = fieldnames(mat_struct); + +for n=1:numel(names) + if isstruct(getfield(mat_struct,names{n})) + docNewElem = docNode.createElement(names{n}); + docNewElem = struct_2_xmlNode(docNode, docNewElem, getfield(mat_struct,names{n})); + docElem.appendChild(docNewElem); + elseif iscell(getfield(mat_struct,names{n})) + cellfield = getfield(mat_struct,names{n}); + for m=1:numel(cellfield) + if ischar(cellfield{m}) + docNewElem = docNode.createElement(names{n}); + docNewElem.appendChild(docNode.createTextNode(cellfield{m})); + docElem.appendChild(docNewElem); + else + docNewElem = docNode.createElement(names{n}); + docNewElem = struct_2_xmlNode(docNode, docNewElem, cellfield{m}); + docElem.appendChild(docNewElem); + end + end + elseif isempty(getfield(mat_struct,names{n})) + %do nothing... + elseif isnumeric(getfield(mat_struct,names{n})) + number = getfield(mat_struct,names{n}); + str = vector2str(number,15); + docNewElem = docNode.createElement(names{n}); + docNewElem.appendChild(docNode.createTextNode(str)); + docElem.appendChild(docNewElem); + elseif ischar(getfield(mat_struct,names{n})) + docNewElem = docNode.createElement(names{n}); + docNewElem.appendChild(docNode.createTextNode(getfield(mat_struct,names{n}))); + docElem.appendChild(docNewElem); + end +end diff --git a/CSXCAD/matlab/private/vector2str.m b/CSXCAD/matlab/private/vector2str.m new file mode 100644 index 0000000..5e3685b --- /dev/null +++ b/CSXCAD/matlab/private/vector2str.m @@ -0,0 +1,35 @@ +function str = vector2str(vec, acc) +% str = vector2str( vec [, acc] ) +% +% internal function for CSXCAD +% +% (c) 2012-2015 Thorsten Liebig <thorsten.liebig@gmx.de> +% (C) 2012 Sebastian Held <sebastian.held@gmx.de> + +str = ''; + +if (numel(vec)==0) + return +end + +if (nargin<2) + acc = 9; +end + +% Octave < 3.8.0 has a bug in num2str() (only 64bit versions) +% see bug report https://savannah.gnu.org/bugs/index.php?36121 +[isOct,version] = isOctave(); +if (isOct && version.major <= 3 && version.minor < 8) + % affected Octave versions + for n = 1:numel(vec) + str = [str sprintf('%.*g', acc, vec(n)) ',']; + end +else + % Matlab and non affected Octave + for n = 1:numel(vec) + str = [str num2str(vec(n),acc) ',']; + end +end + +% remove the last ',' +str = str(1:end-1); diff --git a/CSXCAD/matlab/searchBinary.m b/CSXCAD/matlab/searchBinary.m new file mode 100644 index 0000000..6310dde --- /dev/null +++ b/CSXCAD/matlab/searchBinary.m @@ -0,0 +1,47 @@ +function [binary_location] = searchBinary(name, searchpath, err_fail) +% [binary_location] = searchBinary(name, searchpath<, err_fail>) +% +% Function to search for executable. If the executable isn't found +% in searchpath, look in environment search path. +% +% parameter: +% name: name of the binary to search +% searchpath: (list of) search paths +% err_fail: 0/1, throw an error if binary is not found (default is 1) +% +% Output: +% binary_location: +% if found in searchpath: full path of binary +% if found in environment search path: name of binary +% if not found: empty string +% +% +% openEMS matlab/octave interface +% ----------------------- +% (C) 2013 Stefan Mahr <dac922@gmx.de> + +if (nargin<3) + err_fail = 1; +end + +if ischar(searchpath) + searchpath = {searchpath}; +end + +% append PATH search paths +searchpath = [searchpath regexp(getenv('PATH'), pathsep, 'split')]; + +% try all search paths +for n=1:numel(searchpath) + binary_location = [searchpath{n} name]; + if exist(binary_location, 'file') + return + end +end + +% binary not found +binary_location = ''; + +if (err_fail) + error('CSXCAD:binary_location', [name ' binary not found!']); +end diff --git a/CSXCAD/matlab/struct_2_xml.m b/CSXCAD/matlab/struct_2_xml.m new file mode 100644 index 0000000..73c0a30 --- /dev/null +++ b/CSXCAD/matlab/struct_2_xml.m @@ -0,0 +1,23 @@ +function struct_2_xml(filename, xml_struct, rootName); + +if ~isOctave() + docNode = com.mathworks.xml.XMLUtils.createDocument(rootName); + docElem = docNode.getDocumentElement; + + docElem = struct_2_xmlNode(docNode, docElem, xml_struct); + + % Save the sample XML document. + xmlFileName = [filename]; + xmlwrite(xmlFileName,docNode); +else +% % for octave you need the octave_xmltoolbox (C) 2007 Thomas Geiger +% % http://wiki.octave.org/wiki.pl?XMLToolboxPort +% xml_struct = octave_struct2xml( xml_struct, rootName ); +% xml_save( filename, xml_struct, 'any' ); + +% xml_toolbox is buggy (sequence of elements is not preserved) + fid = fopen( filename, 'w' ); + fprintf( fid, '<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>\n' ); + fprintf( fid, octave_struct2xml_2(xml_struct,rootName,'') ); + fclose( fid ); +end |