diff --git a/lib/AddBlocks.m b/lib/AddBlocks.m
new file mode 100644
index 0000000..d8a6acb
--- /dev/null
+++ b/lib/AddBlocks.m
@@ -0,0 +1,33 @@
+function [object_blk,blk_pos_corner ]=AddBlocks(model,blk)
+%Add blocks, used in model expansion and loading rectangular sources on the topographic surface
+
+blk_pos_corner = [];
+ins0 = 0;
+for n = 1:size(blk.Lx,2)
+ %Size unit：m
+ lenx = blk.Lx(n);
+ leny = blk.Ly(n);
+ lenz = blk.Lz(n);
+ %Center position coordinate
+ xc = blk.CentCoord(n,1);
+ yc = blk.CentCoord(n,2);
+ zc = blk.CentCoord(n,3);
+
+ blk_position = [xc yc zc];
+ blk_size = [lenx leny lenz];
+ blkLabel = ['blk' num2str(n+ins0)];
+
+ model.geom('geom1').feature.create(blkLabel,'Block');
+ model.geom('geom1').feature(blkLabel).set('size',blk_size);
+ model.geom('geom1').feature(blkLabel).set('pos',blk_position);
+ model.geom("geom1").feature(blkLabel).set("rot", blk.angle);
+ model.component("mod1").geom("geom1").feature( blkLabel).set("base", "center");
+
+ object_blk{n} = blkLabel;
+ temp = [];
+ blk_pos_corner = cat(1, blk_pos_corner, temp);
+end
+
+model.component("mod1").geom("geom1").run();
+
+
diff --git a/lib/AddRecParametricSurface.m b/lib/AddRecParametricSurface.m
new file mode 100644
index 0000000..5b546dc
--- /dev/null
+++ b/lib/AddRecParametricSurface.m
@@ -0,0 +1,34 @@
+function [object_rec,ps_pos_rec ]=AddRecParametricSurface(model,rec,lengthcurve,heightcurve)
+%
+%A vertical auxiliary parametric surfaces are added for constructing the
+% receives on the terrain surface
+% ins0 = 0;
+% xrec=[2 10];
+% yrec=[0 0];
+% zrec=[0 0];
+% lengthcurve=[5 5];
+% heightcurve=5;
+
+x = rec(:,1);
+y = rec(:,2);
+z = rec(:,3);
+nrec = length(x);
+object_rec=cell(1,nrec);%{};
+ps_pos_rec=[];
+ins0 = 0;
+for n = 1:nrec
+ psLabel = ['ps' num2str(n+ins0)];
+ ps(n+ins0)= model.geom('geom1').feature.create( psLabel ,'ParametricSurface');
+ model.geom('geom1').feature( psLabel ).set('parmin1',num2str(x(n)-lengthcurve/2));
+ model.geom('geom1').feature( psLabel ).set('parmin2',num2str(z(n)-heightcurve));
+ model.geom('geom1').feature( psLabel ).set('parmax1',num2str(x(n)+lengthcurve/2));
+ model.geom('geom1').feature( psLabel ).set('parmax2',num2str(z(n)+heightcurve));
+ model.geom('geom1').feature( psLabel ).set('coord',{'s1',num2str(y(n)),'s2'});
+ model.geom('geom1').feature( psLabel ).set('maxknots',{'4'});
+ object_rec{n} = psLabel;
+ ps_pos_temp = [x(n)-lengthcurve/2 y(n) z(n)-heightcurve x(n)-lengthcurve/2 y(n) z(n)+heightcurve];
+ ps_pos_rec = cat(1, ps_pos_rec, ps_pos_temp);
+end
+model.component("mod1").geom("geom1").run();
+
+
diff --git a/lib/AndInterpolationCurve.m b/lib/AndInterpolationCurve.m
new file mode 100644
index 0000000..4f601cd
--- /dev/null
+++ b/lib/AndInterpolationCurve.m
@@ -0,0 +1,13 @@
+function objectIPC = AndInterpolationCurve(model,filenames)
+%Interpolating curves are added, which can be used to construct 3D irregular volumes
+for n = 1:size(filenames,1)
+ IPCname =['ipc' num2str(n)];
+ model.geom('geom1').create(IPCname, 'InterpolationCurve');
+ model.geom('geom1').feature(IPCname).set('type', 'closed');
+ model.geom('geom1').feature(IPCname).set('source', 'file');
+ model.geom('geom1').feature(IPCname).set('filename', filenames{n});
+ model.geom('geom1').feature(IPCname).set('struct', 'sectionwise');
+ objectIPC{n} = IPCname;
+end
+model.geom('geom1').run;
+
diff --git a/lib/AndPoint.m b/lib/AndPoint.m
new file mode 100644
index 0000000..f7fb786
--- /dev/null
+++ b/lib/AndPoint.m
@@ -0,0 +1,13 @@
+
+function objectIPC = AndPoint(model,p)
+% Add points that can be used for the endpoints of 3D irregular volumes
+for np = 1:size(p,1)
+ PTname =['pt' num2str(n)];
+ model.geom("geom1").create(PTname, "Point");
+ model.geom("geom1").feature(PTname).setIndex("p", p(np,1), 0);
+ model.geom("geom1").feature(PTname).setIndex("p", p(np,2), 1);
+ model.geom("geom1").feature(PTname).setIndex("p", p(np,3), 2);
+ objectPT{n} = PTname;
+end
+model.geom('geom1').run;
+end
diff --git a/lib/Comsol_with_Matlab_Start.m b/lib/Comsol_with_Matlab_Start.m
new file mode 100644
index 0000000..0b59ab6
--- /dev/null
+++ b/lib/Comsol_with_Matlab_Start.m
@@ -0,0 +1,19 @@
+
+
+% Start the 'COMSOL Multiphysics with MATLAB' interfaces
+% Required: '...\comsolmphserver.exe' and '...\Multiphysics\mli' file address
+% You can also manually launch the executable COMSOL Multiphysics with MATLAB. exe
+
+path = pwd;
+
+try
+ mphtags -show
+ warning('Already connected to a server!');
+catch
+ winopen('D:\Software\COMSOL60\Multiphysics\bin\win64\comsolmphserver.exe');
+ cd 'D:\Software\COMSOL60\Multiphysics\mli';
+ mphstart;
+end
+cd(path);
+
+
diff --git a/lib/Variables 1.txt b/lib/Variables 1.txt
new file mode 100644
index 0000000..8e7d0a8
--- /dev/null
+++ b/lib/Variables 1.txt
@@ -0,0 +1,9 @@
+rho_xy ((abs((Ex2_G*mf.Hx-Ex1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, xy"
+rho_yx ((abs((Ey1_G*mf2.Hy-Ey2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, yx"
+rho_xx ((abs((Ex1_G*mf2.Hy-Ex2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, xx"
+rho_yy ((abs((Ey2_G*mf.Hx-Ey1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy)))^2/(2*pi*freq*mu0_const)) "Apparent resistivity, yy"
+phi_xy arg(1[S]*(Ex2_G*mf.Hx-Ex1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, xy"
+phi_yx arg(1[S]*(Ey1_G*mf2.Hy-Ey2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, yx"
+phi_xx arg(1[S]*(Ex1_G*mf2.Hy-Ex2_G*mf.Hy)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, xx"
+phi_yy arg(1[S]*(Ey2_G*mf.Hx-Ey1_G*mf2.Hx)/(mf.Hx*mf2.Hy-mf2.Hx*mf.Hy))[rad] "Apparent resistivity phase, yy"
+S abs((mf2.Ex/mf2.Hx+mf.Ey/mf.Hy)/(mf.Ex/mf.Hy-mf2.Ey/mf2.Hx)) ""
diff --git a/lib/Variables 2.txt b/lib/Variables 2.txt
new file mode 100644
index 0000000..9a2f371
--- /dev/null
+++ b/lib/Variables 2.txt
@@ -0,0 +1,12 @@
+Ex1_G real(mf.Ex)+mf.omega*1[s]*imag(d(mf.psi0,x))/1[S/m]+i*(imag(mf.Ex)+mf.omega*1[s]*real(d(mf.psi0,x))/1[S/m]) ""
+Ey1_G real(mf.Ey)+mf.omega*1[s]*imag(d(mf.psi0,y))/1[S/m]+i*(imag(mf.Ey)+mf.omega*1[s]*real(d(mf.psi0,y))/1[S/m]) ""
+Ez1_G real(mf.Ez)+mf.omega*1[s]*imag(d(mf.psi0,z))/1[S/m]+i*(imag(mf.Ez)+mf.omega*1[s]*real(d(mf.psi0,z))/1[S/m]) ""
+Ex2_G real(mf2.Ex)+mf2.omega*1[s]*imag(d(mf2.psi0,x))/1[S/m]+i*(imag(mf2.Ex)+mf2.omega*1[s]*real(d(mf2.psi0,x))/1[S/m]) ""
+Ey2_G real(mf2.Ey)+mf2.omega*1[s]*imag(d(mf2.psi0,y))/1[S/m]+i*(imag(mf2.Ey)+mf2.omega*1[s]*real(d(mf2.psi0,y))/1[S/m]) ""
+Ez2_G real(mf2.Ez)+mf2.omega*1[s]*imag(d(mf2.psi0,z))/1[S/m]+i*(imag(mf2.Ez)+mf2.omega*1[s]*real(d(mf2.psi0,z))/1[S/m]) ""
+normE1_G sqrt(Ex1_G^2+Ey1_G^2+Ez1_G^2) ""
+normE2_G sqrt(Ex2_G^2+Ey2_G^2+Ez2_G^2) ""
+normE2_G2 sqrt(Ey2_G^2+Ez2_G^2) ""
+Eyz_r (sqrt(real(Ey2_G)^2+real(Ez2_G)^2)) ""
+Eyz_i (sqrt(imag(Ey2_G)^2+imag(Ez2_G)^2)) ""
+Exy_i (sqrt(imag(Ex2_G)^2+imag(Ey2_G)^2)) ""
diff --git a/lib/get_curveTxtFile.m b/lib/get_curveTxtFile.m
new file mode 100644
index 0000000..da0acad
--- /dev/null
+++ b/lib/get_curveTxtFile.m
@@ -0,0 +1,17 @@
+
+%Gets the path of a file with a specific character in the folder
+function filenames = get_curveTxtFile(data_dir,id,str)
+% data_dir：folder path
+% id：The position of the character
+% Finding characters
+% example: curve_01.txt;curve_02.txt
+% CurveFiles = get_curveTxtFile(data_dir,[1:5],'curve');
+
+D = dir(data_dir );
+nf = 0;
+for i=3:numel(D)
+ if strcmp(D(i).name(id),str)
+ nf = nf+1;
+ filenames{nf,1} = [data_dir '\' D(i).name] ;
+ end
+end
\ No newline at end of file
diff --git a/lib/interpclosed.m b/lib/interpclosed.m
new file mode 100644
index 0000000..6e5d6a1
--- /dev/null
+++ b/lib/interpclosed.m
@@ -0,0 +1,420 @@
+function varargout = interpclosed(x,y,varargin)
+% INTERPCLOSED Arc-length interpolation, perimeter and area of 2D closed curves defined by points
+%
+% xyq = INTERPCLOSED(x,y,tq) Interpolates new data points xyq at given
+% positions defined by an arc-length parametrization and the query points
+% tq, along the closed curve defined by the points specified by x and y.
+% The default method 'spline' is used. tq is a n-by-1 array with its
+% elements constrained within [0,1], with 0 being the first point of the
+% curve and 1 the last one.
+%
+% [len,area] = INTERPCLOSED(x,y) If tq is not specified and two output
+% variables are requested, then only the perimeter and area of the
+% interpolated curve are returned. Both outputs are obviously scalars.
+%
+% [len,area,c] = INTERPCLOSED(x,y) If tq is not specified and three output
+% variables are requested, then the perimeter, area and centroid of the
+% interpolated curve are returned. The centroid is a vector with the
+% position as (x,y).
+%
+% [len,area,c,I] = INTERPCLOSED(x,y) If tq is not specified and four output
+% variables are requested, then the perimeter, area, centroid, and second
+% moments of area of the interpolated curve are returned. The second moment
+% of area is a vector with three elements: (Ixx, Iyy, Ixy).
+%
+% pp = INTERPCLOSED(x,y,'pp') If only one output argument is defined and
+% the string 'pp' is given as input, the returned variable is the piesewise
+% polynomial pp, regardless of the definition of tq.
+%
+% [xyq,len,area] = INTERPCLOSED(x,y,tq) If tq is provided and there are
+% three output variables, the perimeter and area are additionally returned.
+%
+% [xyq,len,area,c,I] = INTERPCLOSED(x,y,tq) If tq is provided and there are
+% five output variables, the perimeter, the area, the centroid, and the
+% second moments of area are additionally returned.
+%
+% [___] = INTERPCLOSED(___,method) By specifying the string method it
+% is possible to change how the interpolated points are obtained.
+% 'linear': Linear interpolation. The interpolated value at a query point
+% is based on linear interpolation of the values at neighboring
+% points in each respective dimension. This is the fastest
+% method.
+% 'spline': Spline interpolation using periodic end conditions. The
+% interpolated value at a query point is based on a cubic
+% interpolation of the values at neighboring points in each
+% respective dimension. This is the default method.
+% 'pchip': Shape-preserving piecewise cubic interpolation. The
+% interpolated value at a query point is based on a shape-
+% preserving piecewise cubic periodic interpolation of the
+% values at neighboring points.
+%
+% [___] = INTERPCLOSED(___,print) By setting the boolean variable print to
+% true, more details about the interpolation can be obtained.
+%
+% Examples:
+% 1) Use the spline method to interpolate at 64 query points within the
+% 8 points used to sample the original circle:
+%
+% t = linspace(0,2*pi,9);
+% x = sin(t) + 0.2; y = cos(t) - 0.2;
+% [len,area,c,I] = interpclosed(x,y);
+% fprintf(['Perimeter: %4.5f, Area: %4.5f,\n',...
+% 'Centroid: [%4.5f %4.5f], Iz: %4.5f\n',...
+% 'To be compared to\n',...
+% '2*pi: %4.5f, pi: %4.5f,\n',...
+% 'Centroid: [%4.5f %4.5f], pi/2: %4.5f\n'],...
+% len,area,c,(I(1)-area*c(2)^2+I(2)-area*c(1)^2),...
+% 2*pi,pi,[0.2,-0.2],pi/2)
+%
+% 2) Get the piecewise polyonomial of a linear interpolation of a given
+% set of points and use the polynomial to make a plot:
+%
+% x = [0 .82 .92 0 -.92 -.82]; y = [.66 .9 0 -.83 0 .9];
+% pp = interpclosed(x,y,'pp','linear');
+% tq = min(pp.breaks):0.001:max(pp.breaks);
+% xyq = ppval(pp,tq);
+% figure, plot(xyq(1,:),xyq(2,:))
+%
+% Find more examples in the File Exchange website.
+%
+% See also CSCVN, PCHIP, MKPP, PPVAL, INTERPARC, ARCLENGTH, INTERP1.
+%
+% Author: Santiago M. Benito
+% Ruhr-Universität Bochum
+% -------------------------------------------------------------------------
+% Contact: santiago.benito@rub.de
+% -------------------------------------------------------------------------
+% Current version: 3.0
+% -------------------------------------------------------------------------
+% Last updated: 17.05.2021
+% Changes:
+% - It is now possible to compute the second moments of area of the
+% fit with this function.
+
+%% Manage input, output and catch eventual problems.
+% Check for errors in the given inputs.
+if nargin < 2
+ error('INTERPCLOSED:insufficientarguments', ...
+ 'At least x and y must be supplied.')
+end
+
+if ~isvector(x) || ~isvector(y) || (length(x) ~= length(y))
+ error('INTERPCLOSED:baddimension', ...
+ 'x and y must be vectors of the same length.')
+end
+
+% Set defaults.
+method = 'spline';
+print = false;
+geomcalc = false; tqgiven = false;
+pp = false;
+
+% Initialize output variables.
+len = 0;
+area = 0;
+c = zeros(1,2);
+Ixx = 0; Iyy = 0; Ixy = 0;
+
+% Check for other input arguments.
+if numel(varargin) > 0
+ % At least one other argument was supplied.
+ for ii = 1:numel(varargin)
+ arg = varargin{ii};
+ if ischar(arg)
+ % It can be the method or the 'pp'-flag.
+ validstrings = {'pp','linear' 'pchip' 'spline'};
+ ind = strncmp(arg,validstrings,2);
+ if isempty(ind) || (sum(ind) == 0) || (sum(ind) > 1)
+ error('INTERPCLOSED:invalidmethod', ...
+ ['Invalid method indicated. Only ''linear'',',...
+ '''pchip'',''spline'' allowed.'])
+ end
+ if ind(1) == 1
+ pp = true;
+ else
+ method = validstrings{ind>0};
+ end
+ elseif islogical(arg)
+ % It must be the print variable, set the print sampling distance.
+ if ~tqgiven, tq = 0:1/32:1; end
+ print = arg;
+ else
+ % It must be tq, defining the parametric arc-length query
+ % points
+ if ~isnumeric(arg)
+ error('INTERPCLOSED:badtq', ...
+ 'tq must be numeric.')
+ else
+ if max(arg) > 1 || min(arg) < 0
+ error('INTERPCLOSED:badtq', ...
+ 'tq elements must be bigger than 0 and smaller than 1.')
+ end
+ tqgiven = true;
+ tq = arg;
+ end
+ end
+ end
+end
+
+% If only one output variable is requested and the 'pp' flag was given, no
+% need to compute the interpolations, regardless of the definition of tq.
+% If three are given, geometry computations will be needed. If two or three
+% are given, but tq was not provided, also compute the geometry computations.
+if nargout == 1 && pp && ~print
+ tqgiven = false;
+elseif (nargout == 2 || nargout == 3 || nargout == 4) && ~tqgiven
+ geomcalc = true;
+elseif (nargout == 3 || nargout == 4 || nargout == 5)
+ geomcalc = true;
+ if ~tqgiven
+ error('INTERPCLOSED:badtq', ...
+ 'tq was not defined and is needed for interpolation.')
+ end
+end
+
+% Be sure everything is formatted correclty and group it.
+x = x(:)'; y = y(:)';
+points = [x;y];
+
+% Round to the 15th decimal position to avoid rounding errors. This is
+% necessary for the function to recongnize start and ending points
+% properly.
+points = round(points,15);
+
+% If the set of points does not describe a closed loop, close it.
+if sum(points(:,1) ~= points(:,end)) > 0
+ points = [points,points(:,1)];
+end
+
+% If less than three distinct points are given, no closed curve can be
+% formed.
+d = sum((diff(points.').^2).');
+if numel(x) - sum(d==0)-1 < 2
+ error('INTERPCLOSED:baddimension', ...
+ 'x and y must be vectors describing at least three distinct points.')
+end
+
+%% Actual program start
+% Compute the coefficients of the fit-polynomials according to the user's
+% choice.
+if strcmpi(method,'linear')
+ % Remove segments with length equal to zero, the linear interpolation
+ % has no continuous derivatives anyway.
+ points(:,d==0) = [];
+
+ % Compute the linear coefficients of the parametric versions of the
+ % lines. First compute the lengths of each segment, then the cumulative
+ % length and finally use the slope in each direction to get the coefs.
+ seglen = sqrt(sum(diff(points,[],2).^2,1));
+ cumarc = [0,cumsum(seglen)];
+ coefX = [diff(points(1,:))./diff(cumarc);points(1,1:(end-1))];
+ coefY = [diff(points(2,:))./diff(cumarc);points(2,1:(end-1))];
+
+ % Create a piecewise polynomial with the given coefficients.
+ coefs = zeros(size([coefX,coefY]));
+ coefs(:,1:2:end) = coefX;
+ coefs(:,2:2:end) = coefY;
+ curve = mkpp(cumarc,coefs',2);
+
+ % Provide the differentiation array for later use.
+ diffarray = [0 0 1;0 0 0];
+
+ % Since we already have the lenghts of the individual segments, just
+ % sum everything up and save some time.
+ len = sum(seglen);
+
+elseif strcmpi(method,'spline')
+ % MATLAB(R) already has a very useful function that makes all the work
+ % for us.
+ curve = cscvn(points);
+
+ % Provide the differentiation array for later use.
+ diffarray = [3 0 0;0 2 0;0 0 1;0 0 0];
+
+elseif strcmpi(method,'pchip')
+ % Like in the function CSCVN, if the user specified a point where the
+ % 2nd derivative is equal to zero, we have to be able to handle the
+ % situation.
+ d = sum((diff(points.').^2).');
+
+ if all(d > 0)
+ % The fit is periodic. To have the start and end slopes equal to
+ % each other, some tricks must be done. Extra points will be added
+ % right before the start and right after the end of the data set.
+ % The fit will be performed with these points, and then the extra
+ % pieces will be removed from the general fit.
+ %pointsNew = [x(end-2:end-1),x,x(2:3);y(end-2:end-1),y,y(2:3)];
+ pointsNew = [points(:,end-2:end-1),points,points(:,2:3)];
+
+ % We need the arc length of the modified dataset, therefore we will
+ % compute it here.
+ seglen = sqrt(sum(diff(pointsNew,[],2).^2,1));
+ cumarc = [0,cumsum(seglen)];
+
+ % Fit coefficients are obtained from the MATLAB(R) original pchip
+ % function.
+ temp = pchip(cumarc,pointsNew(1,:)); coefX = temp.coefs;
+ temp = pchip(cumarc,pointsNew(2,:)); coefY = temp.coefs;
+
+ % Here we remove the unnecesary pieces by removing the extra
+ % coefficients.
+ coefs = zeros(size([coefX;coefY])-[8,0]);
+ coefs(1:2:end,:) = coefX(3:end-2,:);
+ coefs(2:2:end,:) = coefY(3:end-2,:);
+
+ % Compute the actual arc length
+ seglen = sqrt(sum(diff(points,[],2).^2,1));
+ cumarc = [0,cumsum(seglen)];
+
+ else
+ % The 1st derivatives at the end points and at the specified points
+ % are not equal, while analysed from both sides. Firstly compute
+ % the arclength of the point distribution.
+ seglen = sqrt(sum(diff(points,[],2).^2,1));
+ cumarc = [0,cumsum(seglen)];
+
+ % Fit coefficients are obtained from the MATLAB(R) original pchip
+ % function, according to the desired derivative contiguity.
+ dp = find(d>0);
+ dpbig = find(diff(dp)>1);
+ dpbig = [dpbig,length(dp)];
+ idx = dp(1):(dp(dpbig(1))+1);
+ temp = pchip(cumarc(idx),points(1,idx)); coefX = temp.coefs;
+ temp = pchip(cumarc(idx),points(2,idx)); coefY = temp.coefs;
+ for j=2:length(dpbig)
+ idx = dp(dpbig(j-1)+1):(dp(dpbig(j))+1);
+ temp = pchip(cumarc(idx),points(1,idx));
+ coefX = [coefX;temp.coefs];
+ temp = pchip(cumarc(idx),points(2,idx));
+ coefY = [coefY;temp.coefs];
+ end
+
+ % Compiling the coefficients in a simple array.
+ coefs = zeros(size([coefX;coefY]));
+ coefs(1:2:end,:) = coefX(1:end,:);
+ coefs(2:2:end,:) = coefY(1:end,:);
+
+ % Update the cumulative arclength
+ cumarc(:,d==0) = [];
+ end
+
+ % Finally compute the piecewise polynomial.
+ curve = mkpp(cumarc,coefs,2);
+
+ % Provide the differentiation array for later use.
+ diffarray = [3 0 0;0 2 0;0 0 1;0 0 0];
+end
+
+% If tq is given (or a print is required), compute the interpolation using
+% the piecewise evaluation function provided in MATLAB(R) and then convert
+% the parametrization into an arc-lenght one.
+if tqgiven || print
+ step = (max(curve.breaks)-min(curve.breaks))/numel(tq)/30;
+ auxtq = min(curve.breaks):step:max(curve.breaks);
+ xyq = ppval(curve,auxtq);
+ tqp = pdearcl(auxtq,xyq,tq,0,1);
+ xyq = ppval(curve,tqp);
+end
+
+% If the geometric parameters (perimeter and area) are required, compute
+% them using some calculus.
+if geomcalc
+ for ii = 1:curve.pieces
+ % Get the coefficients of the piecewise polynomial expresions of
+ % the parametric form.
+ coefX = curve.coefs(2*ii-1,:);
+ coefY = curve.coefs(2*ii,:);
+
+ % Obtain the derivatives of the polynomials.
+ difX = coefX*diffarray;
+ difY = coefY*diffarray;
+
+ % The length in the linear case is already computed, skip this bit.
+ if ~strcmpi(method,'linear')
+ % Define the function employed in the arc length and integrate
+ % it.
+ flen = @(t) sqrt(polyval(difX,t-curve.breaks(ii)).^2 ...
+ + polyval(difY,t-curve.breaks(ii)).^2);
+ len = len + integral(flen,curve.breaks(ii),curve.breaks(ii+1));
+ end
+ % The area integral is computed here.
+ farea = @(t) polyval(conv(coefY,difX),...
+ t-curve.breaks(ii));
+ area = area + integral(farea,curve.breaks(ii),...
+ curve.breaks(ii+1));
+
+ % The centroid is computed here
+ fcx = @(t) polyval(conv(coefX,conv(coefY,difX)),...
+ t-curve.breaks(ii));
+ c(1) = c(1) + integral(fcx,curve.breaks(ii),...
+ curve.breaks(ii+1));
+ fcy = @(t) polyval(conv(coefY,conv(coefX,difY)),...
+ t-curve.breaks(ii));
+ c(2) = c(2) - integral(fcy,curve.breaks(ii),...
+ curve.breaks(ii+1));
+
+ % The area moments of inertia
+ fIxx = @(t) polyval(conv(coefY,conv(coefY,conv(coefY,difX))),...
+ t-curve.breaks(ii));
+ Ixx = Ixx + integral(fIxx,curve.breaks(ii),...
+ curve.breaks(ii+1));
+ fIyy = @(t) polyval(conv(coefX,conv(coefX,conv(coefX,difY))),...
+ t-curve.breaks(ii));
+ Iyy = Iyy - integral(fIyy,curve.breaks(ii),...
+ curve.breaks(ii+1));
+ fIxy = @(t) polyval(conv(coefY,conv(coefY,conv(coefX,difX))),...
+ t-curve.breaks(ii));
+ Ixy = Ixy + integral(fIxy,curve.breaks(ii),...
+ curve.breaks(ii+1));
+ end
+ c = c / area;
+ I = [1/3*Ixx,1/3*Iyy,1/2*Ixy]*sign(area);
+ area = abs(area);
+
+
+end
+
+%% If required, print some figures to show what the algorithm did.
+if print
+ figure
+ subplot(1,2,1)
+ plot(xyq(1,:),xyq(2,:),'*')
+ hold on
+ plot(points(1,:),points(2,:),'o')
+ plot(c(1),c(2),'x')
+ xlabel('x'), ylabel('y'), hold off, axis equal
+ title('Cartesian representation'), legend('Interpolation','Points',...
+ 'Centroid')
+ subplot(1,2,2)
+ plot(tq,xyq), hold on
+ for ii = (curve.breaks)/max(curve.breaks)
+ line([ii ii],ylim,'LineStyle','--','Color','k')
+ line([ii ii],ylim,'LineStyle','--','Color','k')
+ end
+ hold off, xlim([min(tq) max(tq)]), title('Parametric representation')
+ xlabel('t'), ylabel('x(t), y(t)'), legend('x(t)','y(t)')
+end
+
+%% Process adequately the variables to be returned.
+if nargout == 2 && ~tqgiven
+ varargout{1} = len; varargout{2} = area;
+elseif nargout == 3 && ~tqgiven
+ varargout{1} = len; varargout{2} = area; varargout{3} = c;
+elseif nargout == 3
+ varargout{1} = xyq;
+ varargout{2} = len; varargout{3} = area;
+elseif nargout == 4 && tqgiven
+ varargout{1} = xyq;
+ varargout{2} = len; varargout{3} = area; varargout{4} = c;
+elseif nargout == 4 && ~tqgiven
+ varargout{1} = len; varargout{2} = area; varargout{3} = c;
+ varargout{4} = I;
+elseif nargout == 5
+ varargout{1} = xyq;
+ varargout{2} = len; varargout{3} = area; varargout{4} = c;
+ varargout{5} = I;
+elseif nargout == 1 && pp
+ varargout{1} = curve;
+else
+ varargout{1} = xyq;
+end
\ No newline at end of file
diff --git a/lib/lofting.m b/lib/lofting.m
new file mode 100644
index 0000000..93a2cb9
--- /dev/null
+++ b/lib/lofting.m
@@ -0,0 +1,30 @@
+function model = lofting(model,data_dir)
+%Constructed Irregular 3D volumes from 2D contour curves
+% example:
+% Comsol_with_Matlab_Start;
+% import com.comsol.model.util.*
+% model = ModelUtil.create('Model1');% ModelUtil.remove('Model');
+% model.modelNode.create('mod1');
+% model.geom.create('geom1', 3);
+% model.mesh.create('mesh1', 'geom1');
+% data_dir = pwd ;
+% model = lofting(model,data_dir)
+
+CurveFiles = get_curveTxtFile(data_dir,[1:5],'curve');
+objectIPC = AndInterpolationCurve(model,CurveFiles);
+
+model.geom("geom1").create("loft1", "Loft");
+model.geom("geom1").feature("loft1").selection("profile").set(objectIPC);
+model.geom("geom1").feature("loft1").set("facepartitioning", "grid");
+
+% model.geom("geom1").create("pare1", "PartitionEdges");
+
+try
+model.component("mod1").geom("geom1").run();
+catch
+ warning('The automatic lofting failed, so the Partition Edges had to be added manually.');
+end
+ mphlaunch(model);
+
+
+end
\ No newline at end of file
diff --git a/lib/plotSlice.m b/lib/plotSlice.m
new file mode 100644
index 0000000..0e630ef
--- /dev/null
+++ b/lib/plotSlice.m
@@ -0,0 +1,27 @@
+
+
+figure;
+xyzID = ['X';'Y';'Z'];
+
+scatter(Pn(:,1),Pn( :,2),3,"filled");
+hold on;
+plot(Pn(k,1),Pn(k,2),'g--','LineWidth',2);
+hold on;
+plot(PI(1,:),PI(2,:),'k','LineWidth',2);
+
+
+title(['Silce' num2str(i)],'FontSize',12,'FontWeight','bold');
+xlabel([xyzID(planeID(1)) '(m)']);
+ylabel([xyzID(planeID(2)) '(m)']);
+l=legend('Point cloud slice','Point cloud boundary','Smooth boundary');
+set(l,'Box','off','FontSize',10);
+set(gca,'color','none','linewidth',1,'FontSize',12,'FontWeight','bold');
+set(gcf,'Position', [713.8000 224.2000 404.8000 361.6000]);
+box on;
+
+axis tight
+axis equal;
+% xlim([80,200]);
+% ylim([-120,120]);
+% set(gca,'color','none');
+% set(gcf,'color','none');
diff --git a/lib/smooth1q.m b/lib/smooth1q.m
new file mode 100644
index 0000000..660fa34
--- /dev/null
+++ b/lib/smooth1q.m
@@ -0,0 +1,241 @@
+function [z,s] = smooth1q(y,s,varargin)
+
+%SMOOTH1Q Quick & easy smoothing.
+% Z = SMOOTH1Q(Y,S) smoothes the data Y using a DCT- or FFT-based spline
+% smoothing method. Non finite data (NaN or Inf) are treated as missing
+% values.
+%
+% S is the smoothing parameter. It must be a real positive scalar. The
+% larger S is, the smoother the output will be. If S is empty (i.e. S =
+% []), it is automatically determined by minimizing the generalized
+% cross-validation (GCV) score.
+%
+% Z = SMOOTH1Q(...,'robust') carries out a robust smoothing that
+% minimizes the influence of outlying data.
+%
+% Z = SMOOTH1Q(...,'periodic') assumes that the data to be smoothed must
+% be periodic.
+%
+% [Z,S] = SMOOTH1Q(...) also returns the calculated value for the
+% smoothness parameter S so that you can fine-tune the smoothing
+% subsequently if required.
+%
+% SMOOTH1Q is a simplified and quick version of SMOOTHN for 1-D data. If
+% you want to smooth N-D arrays use SMOOTHN.
+%
+% Notes
+% -----
+% 1) SMOOTH1Q works with regularly spaced data only. Use SMOOTH1 for non
+% regularly spaced data.
+% 2) The smoothness parameter used in this algorithm is determined
+% automatically by minimizing the generalized cross-validation score.
+% See the references for more details.
+%
+% References
+% ----------
+% 1) Garcia D, Robust smoothing of gridded data in one and higher
+% dimensions with missing values. Computational Statistics & Data
+% Analysis, 2010.
+% PDF download
+% 2) Buckley MJ, Fast computation of a discretized thin-plate smoothing
+% spline for image data. Biometrika, 1994.
+% Link
+%
+% Examples:
+% --------
+% % Simple curve
+% x = linspace(0,100,200);
+% y = cos(x/10)+(x/50).^2 + randn(size(x))/10;
+% z = smooth1q(y,[]);
+% plot(x,y,'r.',x,z,'k','LineWidth',2)
+% axis tight square
+%
+% % Periodic curve with ouliers and missing data
+% x = linspace(0,2*pi,300);
+% y = cos(x)+ sin(2*x+1).^2 + randn(size(x))/5;
+% y(150:155) = rand(1,6)*5;
+% y(10:40) = NaN;
+% subplot(121)
+% z = smooth1q(y,1e3,'periodic');
+% plot(x,y,'r.',x,z,'k','LineWidth',2)
+% axis tight square
+% title('Non robust')
+% subplot(122)
+% z = smooth1q(y,1e3,'periodic','robust');
+% plot(x,y,'r.',x,z,'k','LineWidth',2)
+% axis tight square
+% title('Robust')
+%
+% % Limaon
+% t = linspace(0,2*pi,300);
+% x = cos(t).*(.5+cos(t)) + randn(size(t))*0.05;
+% y = sin(t).*(.5+cos(t)) + randn(size(t))*0.05;
+% z = smooth1q(complex(x,y),[],'periodic');
+% plot(x,y,'r.',real(z),imag(z),'k','linewidth',2)
+% axis equal tight
+%
+% See also SMOOTHN, SMOOTH1.
+%
+% -- Damien Garcia -- 2012/08, revised 2014/02/26
+% website: www.BiomeCardio.com
+
+%-- Check input arguments
+error(nargchk(2,4,nargin));
+assert(isvector(squeeze(y)),...
+ ['Y must be a 1-D array. Use SMOOTHN for non vector arrays.'])
+if isempty(s)
+ isauto = 1;
+else
+ assert(isnumeric(s),'S must be a numeric scalar')
+ assert(isscalar(s) && s>0,...
+ 'The smoothing parameter S must be a scalar >0')
+ isauto = 0;
+end
+
+%-- Order (use m>=2, m = 2 is recommended)
+m = 2; % Note: order of the smoothing process, can be modified
+
+%-- Options ('robust' and/or 'periodic')
+isrobust = 0; method = 'dct'; % default options
+%--
+if nargin>2
+ assert(all(cellfun(@ischar,varargin)),...
+ 'The options must be ''robust'' and/or ''periodic''.')
+ varargin = lower(varargin);
+ if nargin==3
+ idx = ismember({'robust','periodic'},varargin);
+ assert(any(idx),...
+ 'The options must be ''robust'' and/or ''periodic''.')
+ if idx(1), isrobust = 1; else method = 'fft'; end
+ else % nargin = 4
+ assert(all(ismember(varargin,{'robust','periodic'})),...
+ 'The options must be ''robust'' and/or ''periodic''.')
+ isrobust = 1;
+ method = 'fft';
+ end
+end
+
+n = length(y);
+siz0 = size(y);
+y = y(:).';
+
+%-- Weights
+W0 = ones(siz0);
+I = isfinite(y); % missing data (NaN or Inf values)
+if any(~I) % replace the missing data (for faster convergence)
+ X = 1:n;
+ x = X(I); xi = X(~I);
+ y(~I) = interp1(x,y(I),xi,'linear','extrap');
+end
+W0(~I) = 0; % weights for missing data are 0
+W = W0;
+
+%-- Eigenvalues
+switch method
+ case 'dct'
+ Lambda = 2-2*cos((0:n-1)*pi/n);
+ case 'fft'
+ Lambda = 2-2*cos(2*(0:n-1)*pi/n);
+end
+
+%-- Smoothing process
+nr = 3; % Number of robustness iterations
+for k = 0:nr*isrobust
+ if isrobust && k>0
+ tmp = sqrt(1+16*s);
+ h = sqrt(1+tmp)/sqrt(2)/tmp;
+ W = W0.*bisquare(y,z,I,h);
+ end
+ if ~all(W==1) % then use an iterative method
+ tol = Inf;
+ zz = y;
+ while tol>1e-3
+ switch method
+ case 'dct'
+ Y = dct(W.*(y-zz)+zz);
+ case 'fft'
+ Y = fft(W.*(y-zz)+zz);
+ end
+ if isauto
+ fminbnd(@GCVscore,-10,30,optimset('TolX',.1));
+ else
+ Gamma = 1./(1+s*Lambda.^m);
+ switch method
+ case 'dct'
+ z = idct(Gamma.*Y);
+ case 'fft'
+ if isreal(y)
+ z = ifft(Gamma.*Y,'symmetric');
+ else
+ z = ifft(Gamma.*Y);
+ end
+ end
+ end
+ tol = norm(zz-z)/norm(z);
+ zz = z;
+ end
+
+ else %---
+ % No missing values, non robust method => Direct fast method
+ %---
+ switch method
+ case 'dct'
+ Y = dct(y);
+ case 'fft'
+ Y = fft(y);
+ end
+ if isauto
+ fminbnd(@GCVscore,-10,30,optimset('TolX',.1));
+ else
+ Gamma = 1./(1+s*Lambda.^m);
+ end
+ switch method
+ case 'dct'
+ z = idct(Gamma.*Y);
+ case 'fft'
+ if isreal(y)
+ z = ifft(Gamma.*Y,'symmetric');
+ else
+ z = ifft(Gamma.*Y);
+ end
+ end
+ end
+end
+
+z = reshape(z,siz0);
+
+ function GCVs = GCVscore(p)
+ s = 10^p;
+ Gamma = 1./(1+s*Lambda.^m);
+ if any(W)
+ switch method
+ case 'dct'
+ z = idct(Gamma.*Y);
+ case 'fft'
+ if isreal(y)
+ z = ifft(Gamma.*Y,'symmetric');
+ else
+ z = ifft(Gamma.*Y);
+ end
+ end
+ RSS = norm(sqrt(W).*(y-z))^2;
+ else % No missing values, non robust method => Direct fast method
+ RSS = norm(Y.*(Gamma-1))^2;
+ end
+ TrH = sum(Gamma);
+ GCVs = RSS/(1-TrH/n)^2;
+ end
+
+end
+
+function W = bisquare(y,z,I,h)
+r = y-z; % residuals
+MAD = median(abs(r(I)-median(r(I)))); % median absolute deviation
+u = abs(r/(1.4826*MAD)/sqrt(1-h)); % studentized residuals
+W = (1-(u/4.685).^2).^2.*((u/4.685)<1); % bisquare weights
+end
\ No newline at end of file
diff --git a/lib/smooth1qExample.m b/lib/smooth1qExample.m
new file mode 100644
index 0000000..302e30d
--- /dev/null
+++ b/lib/smooth1qExample.m
@@ -0,0 +1,30 @@
+ % Simple curve
+ x = linspace(0,100,200);
+ y = cos(x/10)+(x/50).^2 + randn(size(x))/10;
+ z = smooth1q(y,[]);
+ plot(x,y,'r.',x,z,'k','LineWidth',2)
+ axis tight square
+
+ % Periodic curve with ouliers and missing data
+ x = linspace(0,2*pi,300);
+ y = cos(x)+ sin(2*x+1).^2 + randn(size(x))/5;
+ y(150:155) = rand(1,6)*5;
+ y(10:40) = NaN;
+ subplot(121)
+ z = smooth1q(y,1e3,'periodic');
+ plot(x,y,'r.',x,z,'k','LineWidth',2)
+ axis tight square
+ title('Non robust')
+ subplot(122)
+ z = smooth1q(y,1e3,'periodic','robust');
+ plot(x,y,'r.',x,z,'k','LineWidth',2)
+ axis tight square
+ title('Robust')
+
+ % Lima鏾n
+ t = linspace(0,2*pi,300);
+ x = cos(t).*(.5+cos(t)) + randn(size(t))*0.05;
+ y = sin(t).*(.5+cos(t)) + randn(size(t))*0.05;
+ z = smooth1q(complex(x,y),[],'periodic');
+ plot(x,y,'r.',real(z),imag(z),'k','linewidth',2)
+% axis equal tight
\ No newline at end of file