123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990 |
- function [sph_pnt,energy] = CPS_sampling(varargin)
- p = inputParser;
- default.pnt_num = 240; % Number of particles; for this algorithm, more points the better
- default.thre = 'auto'; % Convergence threshold
- default.RepulF = 0.1; % Repulsive force factor, the higher the stronger force
- default.RstF = 1.5; % Resistance force, higher stronger
- default.MaxIter = 500; % Maximum number of iteration if not converged
- default.VisIt = false; % If make visulization
- addOptional(p,'pnt_num',default.pnt_num);
- addOptional(p,'thre',default.thre);
- addOptional(p,'RepulF',default.RepulF,@(n) n>0&&n<=1);
- addOptional(p,'RstF',default.RstF,@(n) n>0&&n<=5);
- addOptional(p,'MaxIter',default.MaxIter,@(n) n>0&&(floor(n)==n));
- addOptional(p,'VisIt',default.VisIt,@(n) islogical(n));
- parse(p,varargin{:});
- pnt_num = p.Results.pnt_num;
- thre = p.Results.thre;
- RepulF = p.Results.RepulF;
- RstF = p.Results.RstF;
- MaxIter = p.Results.MaxIter;
- VisIt = p.Results.VisIt;
- if strcmp(thre,"auto")
- thre = max(1/pnt_num,0.0065);
- end
- sph_pnt = rand(pnt_num,3)-.5; % Randomly assign particle position on 3D space
- sph_pnt = sph_pnt./sqrt(sum(sph_pnt.^2,2)); % Projected particles onto sphere
- sph_pntQN = quaternion([zeros(pnt_num,1),sph_pnt]); % Convert 3-element vector to quaternion
- currentSpeed = 0; % Initial moving speed = 0
- energy = 1; % The energy will be lower when the particles are more evenly distributed
- iterN = 0; % Iteration number
- while energy > thre && iterN <= MaxIter
- % Core for charged particle system
- diff_vec = sph_pntQN - (sph_pntQN.conj)'; % Calculate the differnece vector between each pair of particle
- vec_dist = asin(diff_vec.norm./2).*2; % Calculating the angle between each pair of particles
- proj_diff_vec = (diff_vec+sph_pntQN.*diff_vec.*sph_pntQN)./2; % Project the difference vectors onto the tangent plane for each particle
- vecForce = 1./vec_dist.^2.*proj_diff_vec./proj_diff_vec.norm; % Calculating the repulsive force vector from other particles weighted by 1/(the distance)^2
- vecForce = quaternion(squeeze(nansum(reshape(vecForce.compact,[size(proj_diff_vec),4]),2))); % Calculate the sum of all force vectors as the overall force vector
- vecForce = vecForce./pnt_num.*RepulF; % Adjust the force to an appriate level so the particles won't bounce around or move very slowly
- currentSpeed = currentSpeed+vecForce; % Updating the moving velocity for each particle
- rotAxis = sph_pntQN.*currentSpeed; % Compute the rotation axis for particle movement on sphere
- rotAxis = (rotAxis.compact.*[0 1 1 1]); % Extract only the cross product
- rotAxis = quaternion((rotAxis./sqrt(sum(rotAxis.^2,2))+[1 0 0 0])... % Compute the rotation quaternion
- .*[cos(vecForce.norm./2),sin(vecForce.norm./2).*[1,1,1]]); % vecForce.norm is the rotation angle for the movement
- sph_pntQN = rotAxis.*sph_pntQN.*rotAxis.conj; % Rotate the particles by the rotation quaternion
- currentSpeed = currentSpeed./RstF; % "Resistance force": prevent particles move too much
-
- energy_all = sort(vec_dist); % These two lines are for computing the energy after movement
- energy = std(energy_all(2,:));
- iterN = iterN+1;
- if rem(iterN,50) ==0
- energy
- end
- if iterN>3000&&energy>.01
- break
- end
- % Visualization part below
- if VisIt
- sph_pnt = sph_pntQN.compact;
- sph_pnt = sph_pnt(:,2:4);
- k = delaunayTriangulation(sph_pnt);
- [K,~] = convexHull(k);
- if iterN == 1
- clf
- sct = patch('vertices',sph_pnt,'faces',K,'Facecolor',[0,.25,.5],...
- 'FaceAlpha',0.7,'edgecolor',[0,.1,.3],'linewidth',2); % Visualize the solid formed with these particles
- axis equal vis3d off
- view(120,30)
- elseif rem(iterN,3) == 0
- sct.Faces = K;
- sct.Vertices = sph_pnt;
- end
- title(iterN)
- pause(0.001)
- end
- end
- if energy > thre
- warning("MaxIter reached; Not Converged")
- end
- sph_pnt = sph_pntQN.compact;
- sph_pnt = sph_pnt(:,2:4);
- end
|