%This function calculates an affine Reccurent Fractal Interpolation Function.
%Input:IP, AD, S, N, steps
%IP is a (N+1)x2 matrix representing the interpolation points
%N is the number of mappings and N+1 is the number of points
%S is a Nx1 matrix representing the contraction factors
%AD: is a Nx4 matrix representing the interval that is mapped to each section. 
%i.e. if [...;x_0,x_1,y_0,y_1;...] then the [x_0,x_1] is the interval.
%steps: is the number of iterations
%Output: P
%The algorithm computes and returns N^(steps+1)+1 points
%representing the attractor of the RFIF.

function P=myRFIF(IP,AD,S,N,steps)
%Here are the map parameters.
for i=1:N
    a(i)=(IP(i+1,1)-IP(i,1))/(AD(i,2)-AD(i,1));
	e(i)=(AD(i,2)*IP(i,1)-AD(i,1)*IP(i+1,1))/(AD(i,2)-AD(i,1));
	c(i)=(IP(i+1,2)-IP(i,2))/(AD(i,2)-AD(i,1))-S(i)*(AD(i,4)-AD(i,3))/(AD(i,2)-AD(i,1));
	f(i)=(AD(i,2)*IP(i,2)-AD(i,1)*IP(i+1,2))/(AD(i,2)-AD(i,1))-S(i)*(AD(i,2)*AD(i,3)-AD(i,1)*AD(i,4))/(AD(i,2)-AD(i,1));
end;

A_0=IP;
points=N+1;
E=10^-14;
for i=1:steps
    %for each section
    l=1; 
    for j=1:N
        %find the start of the respective interval
        for start=1:points
            if abs(A_0(start,1)-AD(j,1))<E
                break;
            end;
        end;
        %find the end of the respective interval
        for end_=start+1:points
            if abs(A_0(end_,1)-AD(j,2))<E
                break;
            end;
        end;
        %map all the points to the section j
        for k=start:end_-1
            P(l,1)=a(j)*A_0(k,1)+e(j);
            P(l,2)=c(j)*A_0(k,1)+S(j)*A_0(k,2)+f(j);
            l=l+1;
        end;
    end;
    P(l,1)=A_0(points,1);
    P(l,2)=A_0(points,2);
    points=l;
    A_0=P;
end;
