Penr@}~

clc
close all
clear all
%~
%~untilable wedges for the kicks!!!~
%~
%define phi
P = (sqrt(5)+1)/2;
%define triangles below
%~A~1:P:1
%~B~P:1:P
%triangle A
TRA = [0           0;         %1~x,y
       P           0;         %2~x,y
       cos(pi/5)   sin(pi/5)];%3~x,y
%triangle B
TRB = [0           0;         %1~x,y
       P           0;         %2~x,y
       P*cos(pi/5) P*sin(pi/5)];%3~x,y
%setup fibinacci sequences
NF = 20;%not set up to ensure drawing will be completed
Fu = zeros(NF,1);
Fp = zeros(NF,1);
Fu(1:2) = [1 0];%`D~A~'
Fp(1:2) = [0 1];%,~U~T.
for q = 3:NF
    Fu(q) = Fu(q-1)+Fu(q-2);
    Fp(q) = Fp(q-1)+Fp(q-2);
end
%set up drawing scale
SC = 4;%number of blocks to draw~
%~make sure you print long enough series
%figure number
figNum = 1;
dF = 1;%0~lines,1~filled
%draw seed triangle
figNum = drawTri(TRB,0,[0,0],0,figNum,0,SC,dF);

for j = 1:SC
    scale = zeros(j,1);%scaling factor per level
    N_C = zeros(j,1);%number of bins
    type = cell(j,1);
    alph = cell(j,1);
    xOff = cell(j,1);
    yOff = cell(j,1);
    perm = cell(j,1);
    for n = 1:j
        %scale step
        q = j-n+1;
        %cell number
        Ncell = (2*n+1);
        %grab scale of the level
        scale(n) = Fu(q)+P*Fp(q);
        %set up cell count
        N_C(n) = Fu(Ncell)+Fp(Ncell);
        %setup triangle type table
        type{n} = zeros(N_C(n),1);
        %update type numbers
        num = 1;
        if n == 1
            %first triangle A
            type{n}(1,1) = 1;
            alph{n}(1,1) = 4*pi/5;
            xOff{n}(1,1) = scale(n)*P+scale(n);
            yOff{n}(1,1) = 0;
            perm{n}(1,1) = 0;
            %second triangle B
            type{n}(2,1) = 2;
            alph{n}(2,1) = 3*pi/5;
            xOff{n}(2,1) = scale(n)*P+scale(n);
            yOff{n}(2,1) = 0;
            perm{n}(2,1) = 0;
        else
            for i = 1:N_C(n-1)
                if type{n-1}(i) == 1
                    %A-A
                    ANG = alph{n-1}(i);
                    ROT = [cos(ANG) -sin(ANG);
                           sin(ANG)  cos(ANG)];
                    type{n}(num,1) = 1;
                    alph{n}(num,1) = alph{n-1}(i)+4*pi/5;
                    RAD = [P*scale(n-1);
                                      0];
                    RAD = ROT*RAD;
                    xOff{n}(num,1) = RAD(1)+xOff{n-1}(i);
                    yOff{n}(num,1) = RAD(2)+yOff{n-1}(i);
                    perm{n}(num,1) = 0;
                    num = num + 1;
                    %A-B
                    type{n}(num,1) = 2;
                    alph{n}(num,1) = alph{n-1}(i);
                    RAD = [         0;
                                    0];
                    RAD = ROT*RAD;
                    xOff{n}(num,1) = RAD(1)+xOff{n-1}(i);
                    yOff{n}(num,1) = RAD(2)+yOff{n-1}(i);
                    perm{n}(num,1) = 0;
                    num = num + 1;
                else
                    ANG = alph{n-1}(i);
                    ROT = [cos(ANG+pi/5) -sin(ANG+pi/5);
                           sin(ANG+pi/5)  cos(ANG+pi/5)];
                    %B-A
                    type{n}(num,1) = 1;
                    alph{n}(num,1) = alph{n-1}(i)-4*pi/5;
                    RAD = [scale(n-1);
                                    0];
                    RAD = ROT*RAD;
                    xOff{n}(num,1) = RAD(1)+xOff{n-1}(i);
                    yOff{n}(num,1) = RAD(2)+yOff{n-1}(i);
                    perm{n}(num,1) = 0;
                    num = num + 1;
                    %B-B1
                    ROT = [cos(ANG) -sin(ANG);
                           sin(ANG)  cos(ANG)];
                    type{n}(num,1) = 2;
                    alph{n}(num,1) = alph{n-1}(i)+3*pi/5;
                    RAD = [scale(n-1)*P;
                                      0];
                    RAD = ROT*RAD;
                    xOff{n}(num,1) = RAD(1)+xOff{n-1}(i);
                    yOff{n}(num,1) = RAD(2)+yOff{n-1}(i);
                    perm{n}(num,1) = 0;
                    num = num + 1;
                    %B-B2
                    type{n}(num,1) = 2;
                    alph{n}(num,1) = alph{n-1}(i)+4*pi/5;
                    RAD = [scale(n-1)*P;
                                      0];
                    RAD = ROT*RAD;
                    xOff{n}(num,1) = RAD(1)+xOff{n-1}(i);
                    yOff{n}(num,1) = RAD(2)+yOff{n-1}(i);
                    perm{n}(num,1) = 0;
                    num = num + 1;
                end
            end
        end
        if n == j
            %draw triangles
            for i = 1:N_C(n)
                angA = alph{n}(i);
                dOff = [xOff{n}(i) yOff{n}(i)];
                Prm = perm{n}(i);
                tp = type{n}(i);
                if tp == 1
                    figNum = drawTri(TRA,angA,dOff,Prm,figNum,n,SC,dF);
                elseif tp == 2
                    figNum = drawTri(TRB,angA,dOff,Prm,figNum,n,SC,dF);
                end
            end
        end
    end
end


function [figNum] = drawTri(drawTri,drawAngle,drawOff,Pm,figNum,n,SC,dF)
%rotate sides
R = [cos(drawAngle) -sin(drawAngle);
     sin(drawAngle)  cos(drawAngle)];
drawTri(2,:) = drawTri(2,:)*R';
drawTri(3,:) = drawTri(3,:)*R';
%apply offset
drawTri(:,1) = drawTri(:,1) + drawOff(1);
drawTri(:,2) = drawTri(:,2) + drawOff(2);
%drawing table
drawX = zeros(3,10);
drawY = zeros(3,10);
for q = 1:10
    ROT = [cos(q*pi/5) -sin(q*pi/5);
           sin(q*pi/5)  cos(q*pi/5)];
    vectA = drawTri(1,:)*ROT';
    vectB = drawTri(2,:)*ROT';
    vectC = drawTri(3,:)*ROT';
    vect = [vectA;vectB;vectC];
    drawX(:,q) = vect(:,1);
    drawY(:,q) = vect(:,2);
end
%color scheme
if SC == 1
    PF1 = 0;
    PF2 = 1;
else
    PF1 = sin(pi/2*n/SC);
    PF2 = sin(pi/2*n/SC+pi/2);
end
FCT = 1/sqrt(2);
%draw triangle
figure(figNum)
hold on
for q = 1:10
    for n = 1:3
        if n == 1
            plx = [drawX(1,q) drawX(2,q)];
            ply = [drawY(1,q) drawY(2,q)];
        elseif n == 2
            plx = [drawX(2,q) drawX(3,q)];
            ply = [drawY(2,q) drawY(3,q)];
        else
            plx = [drawX(3,q) drawX(1,q)];
            ply = [drawY(3,q) drawY(1,q)];
        end
        if dF == 0
            plot(plx,ply,'linewidth',2,'color',[(PF1+PF2) PF1 PF2])
        else
            fill(drawX(:,q),drawY(:,q),[FCT*(PF1+PF2) PF1 PF2])
        end
    end
end
axis equal
grid on

figNum = figNum;