clear
zoom = 0.01;
sim.domain = [0 2 0 2]; %[xmin xmax ymin ymax] to set the x-axis and y-axis limits
%!!!reading from .mat file, all needed is wall.corner & wall.corner_num
%and wall.num(real number)
wall.num = 75;
wall.corner_setting = [10,12];  %min and max corner number which polygon has
long = 0.14;
short = 0.03;
wall.size_setting = [long,short; 0.6,0.2];  %1-shrinking rate*rand
%[long side, short side; long shrinking rate, short shrinking rate]
% rho_setting = 1000;%2D kg/m2
%
wall.rho = zeros(1,wall.num);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wall = polygon_generation(wall);%%%%%%%%%%%%%

wall_pos = zeros(2,wall.num);
center_zone = [0 2 0.4 1.4]  ;%[xmin xmax ymin ymax] to set the x-axis and y-axis limits

%rotmat = [cos(theta_rot),-sin(theta_rot);sin(theta_rot),cos(theta_rot)]
for i = 1:wall.num
    theta_rot = 0.25*pi+rand*0.5*pi;
    %   theta_rot = rand*2*pi;
    wall.corner(:,1:wall.corner_num(i)+1,i) = [cos(theta_rot),-sin(theta_rot);sin(theta_rot),cos(theta_rot)]*wall.corner(:,1:wall.corner_num(i)+1,i);
    
    wall.corner(1,1:wall.corner_num(i)+1,i) = wall.corner(1,1:wall.corner_num(i)+1,i)+wall_pos(1,i)+center_zone(1)+rand*(center_zone(2)-center_zone(1));
    wall.corner(2,1:wall.corner_num(i)+1,i) = wall.corner(2,1:wall.corner_num(i)+1,i)+wall_pos(2,i)+center_zone(3)+rand*(center_zone(4)-center_zone(3));
end
%%
i=1;
while i<=wall.num
    
    add = 1;
    
    corner1 = wall.corner(1,1:wall.corner_num(i),i);
    corner2 = wall.corner(2,1:wall.corner_num(i),i);
    
    mat = zeros(4,length(corner1));
    mat(1,:)=(corner1<center_zone(1));
    mat(2,:)=(corner1>center_zone(2));
    mat(3,:)=(corner2<center_zone(3));
    mat(4,:)=(corner2>center_zone(4));
    corner1(mat(1,:)==1)=center_zone(1);
    corner1(mat(2,:)==1)=center_zone(2);
    corner2(mat(3,:)==1)=center_zone(3);
    corner2(mat(4,:)==1)=center_zone(4);
    
    if sum(sum(mat))>0
        for ii = 1:4
            if sum(mat(ii,:))>2
                while mat(ii,1)==1 || mat(ii,end)==0
                    mat = circshift(mat,[0,-1]);
                    corner1 = circshift(corner1,[0,-1]);
                    corner2 = circshift(corner2,[0,-1]);
                end
                nnum = sum(mat(ii,:)==0)+1;
                corner1_new = [corner1(1:nnum) corner1(end)];
                corner2_new = [corner2(1:nnum) corner2(end)];
                nnum = nnum+1;
                
                corner1 = corner1_new;
                corner2 = corner2_new;
                wall.corner_num(i) = nnum;
                add=0;
                break
            end
        end
    end
    wall.corner(:,:,i) = 0;
    wall.corner(1,1:wall.corner_num(i)+1,i)=[corner1 corner1(1)];
    wall.corner(2,1:wall.corner_num(i)+1,i)=[corner2 corner2(1)];
    i=i+add;
end

wall.num = wall.num;
wall.corner = wall.corner.*zoom;
wall.corner_num = wall.corner_num;
%% figure output
hold on;
axis equal;
for i=1:wall.num
    corner_num = wall.corner_num(i)+1;
    fill(wall.corner(1,1:corner_num,i),...
        wall.corner(2,1:corner_num,i),[.6 .6 .6]);%,'EdgeColor','none');
    %       plot(walll.corner(1,1:corner_num,i),...
    %           walll.corner(2,1:corner_num,i),'black');
end
axis(sim.domain.*zoom);      %[xmin xmax ymin ymax] to set the x-axis and y-axis limits

%save('wall.mat','wall');

function [ part ] = polygon_generation( part )

%%%%%%% parameter input %%%%%%%%%%%
% part.num = 10;
% part.corner_setting = [6,10];  %min and max corner number which polygon has
% part.corner = zeros(2,part.corner_setting(2),part.num);
% part.corner_num = zeros(1,part.num);
% part.area = zeros(1,part.num);
% part.pos = zeros(2,part.num);

%%%%%%% parameter initialzation %%%%%%%%%%%
part.corner = zeros(2,part.corner_setting(2),part.num);
part.corner_num = zeros(1,part.num);
part.area = zeros(1,part.num);
part.pos = zeros(2,part.num);
part.I = zeros(1,part.num);      %Moment of inertia
part.bound_box = zeros(4,part.num);

for i=1:part.num
    
    n = randi([part.corner_setting(1) part.corner_setting(2)],1,1);
    theta = linspace(0,2*pi,n+1);
    
    %random theta%
    theta = theta+pi/n*(rand(size(theta))-0.5);
    theta = sort(theta);
    theta(end) = theta(1);
    %end% %need judgement!!!!!%
    
    x = part.size_setting(1,1)*(1-part.size_setting(2,1)*rand)*cos(theta);
    y = part.size_setting(1,2)*(1-part.size_setting(2,2)*rand)*sin(theta);
    
    % %randomly relocate--unavailable now!!!%
    % deltar=0.55;
    % x=a*cos(theta);
    % y=b*sin(theta);
    % x = x+deltar*a*(rand(size(x))-.5);
    % y = y+deltar*b*(rand(size(y))-.5);
    % x(end)=x(1);
    % y(end)=y(1);
    % %end%
    
    corner = zeros(n+1,2);
    corner(:,1) = x';
    corner(:,2) = y';
    corner = corner';
    
    % %rotation%
    % theta_rot = rand(1)*2*pi;
    % rotmat = [cos(theta_rot)  -sin(theta_rot)
    %         sin(theta_rot)   cos(theta_rot)];
    % corner = corner*rotmat;
    % corner = corner';
    % %end%
    
    %get the center of mass
    to_inertia = 1;
    [center,area,inertia] = polygon_center_area(corner,to_inertia); %mass_center function
    corner(1,:) = corner(1,:)-center(1,1);
    corner(2,:) = corner(2,:)-center(2,1);
    %end%
    
    aa = corner.*corner;
    aa2 = (max(aa(1,:)+aa(2,:)))^0.5;
    
    part.bound_box(:,i) = aa2;  %initialize the bounding box!!!!
    
    part.corner(:,1:n+1,i) = corner;
    part.corner_num(1,i) = n;
    part.area(1,i) = area;
    part.I(1,i) = inertia*part.rho(i);
    
    % plot([corner(1,:),0],[corner(2,:),0]);
    % hold on
    % axis equal
    
end
end

function [ center,area,inertia ] = polygon_center_area( corner,to_inertia )

n = length(corner(1,:))-1;
%aver_point = [sum(corner(1,1:n))/n sum(corner(2,1:n))/n];
aver_point = sum(corner(:,1:n),2)'./n;
cm = zeros(2,n);
tri_center = zeros(2,n);
tri_area = zeros(1,n);
tri_J = zeros(1,n);   %moment of inertia J

for i=1:n 
    %%%mass center%
%     tri_area(i) = polyarea([corner(1,i),corner(1,i+1),aver_point(1,1)],...
%         [corner(2,i),corner(2,i+1),aver_point(1,2)]);
    x1=corner(1,i);
    x2=corner(1,i+1);
    x3=aver_point(1,1);
    y1=corner(2,i);
    y2=corner(2,i+1);
    y3=aver_point(1,2);
    tri_area(i) = 0.5*abs(x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3);
    x_e=(x1+x2+x3)/3;
    y_e=(y1+y2+y3)/3;
    tri_center(:,i) = [x_e;y_e];
    cm(:,i) = tri_area(i).*tri_center(:,i);
end

area = sum(tri_area);
% center(1,1) = sum(cm(1,:))./area;
% center(2,1) = sum(cm(2,:))./area;

if area==0
    center = aver_point';
else
    center = sum(cm,2)./area;
end

inertia = [];

if to_inertia==1
 for i=1:n
%     tri_center(:,i) = [(corner(1,i)+corner(1,i+1)+center(1,1))/3
%                      (corner(2,i)+corner(2,i+1)+center(2,1))/3];
%     tri_area(i) = polyarea([corner(1,i),corner(1,i+1),center(1,1)],...
%                      [corner(2,i),corner(2,i+1),center(2,1)]);
    r_ = corner(:,i)-center;
    s_ = corner(:,i+1)-corner(:,i);
    b_= r_+s_;
    b_proj = dot(b_,r_).*r_./(norm(r_))^2;
    b_rej = b_-b_proj;
    a = sign(dot(b_proj,r_))*norm(b_proj);  %maybe costy
    h = norm(b_rej);
    b = norm(r_);
    tri_J(i) = (h*b^3+h*a*b^2+h*(a^2)*b+b*h^3)/12;    %moment of inertia J of tri
 end
inertia = sum(tri_J);
end

end
