% "all_ensembles.m" - Creats all 40x120 Sverdlik/Levanon arrays
clear all
hh_4=[0 0 0 0 ; 1 0 0 0 ; 0 1 0 0  ; 0 0 1 0  ; 0 0 0 1  ]; % first 4 columns
h_last_col=[2 1 0 0 0 ;2 2 1 0 0 ; 2 2 0 2 0 ; 2 0 1 2 0; 2 1 1 2 0; 2 0 2 1 0; 2 2 1 1 0; ...
            2 1 0 0 2; 2 0 2 0 2; 2 2 2 0 2; 2 2 0 2 2; 2 1 2 2 2; 2 2 1 2 2; 2 0 0 1 2; ...
            2 2 2 1 2; 2 0 1 1 2; 2 0 0 0 1; 2 2 0 0 1; 2 1 1 0 1; 2 2 2 2 1; 2 1 0 1 1; 2 0 2 1 1]';
% mp are all the possible multipliers to check linear dependence        
mp=[1 1 1 1; 1 1 1 2; 1 1 2 1; 1 1 2 2; 1 2 1 1; 1 2 1 2; 1 2 2 1; 1 2 2 2; ...
          2 1 1 1; 2 1 1 2; 2 1 2 1; 2 1 2 2; 2 2 1 1; 2 2 1 2; 2 2 2 1]';
[ss,tt]=size(mp);
for pol=1:22
            hh(:,:,pol)=[hh_4 h_last_col(:,pol)];
            aa(:,1,pol)=[1 0 0 0 0]';
            bb(:,1,pol)=[1 0 0 0 0]';
            for k=1:120
                    bb(:,k+1,pol)=hh(:,:,pol)*aa(:,k,pol);
                    aa(:,k+1,pol)=mod(bb(:,k+1,pol),3);
            end
end % of for pol=1:22
% aa are 22 arrays each 5x121.  The columns are all the powers of x in the corresponding polynom
pa=1; pb=2; pc=2;
flag_result=0;
results=input('How many ensembles to search for (1 to 22) = ?  ');
while flag_result<results
    if pc<22
            pc=pc+1;
    else
            pb=pb+1;
            pc=pb+1;
    end
    poly3=[pa  pb  pc];
    disp('   ')
    disp(' Polynomials used   '  ), disp(poly3)
    enssig_temp=[ ];
    enssig=[ ];
    r=[1:120];
    lr=length(r);
    aap=aa(:,:,poly3);
    while lr>11
            m=1+min(r);
            for pol=1:3
                    d=1; % the row number
                    for n=(m+1):120
                            for p=(n+1):121
                                    qq(:,:,pol)=[aap(:,1,pol) aap(:,m,pol) aap(:,n,pol) aap(:,p,pol)]'  ;
                                    t=1;
                                    while t<(tt+1)  % tt is the number of columns of mp (The multiplying options)
                                            qqq(:,:,pol)=diag(mp(:,t))*qq(:,:,pol);
                                            sq=mod(sum(qqq(:,:,pol)),3);
                                            if sq==zeros(1,5)  % implying linear dependence
                                                   t=tt+1; % a multiplier was found. No need to try other multipliers
                                                   dd(d,:,pol)=sort([0  m-1  n-1 p-1]);
                                                   d=d+1;
                                            else
                                                    t=t+1;
                                            end
                                   end  % while  t
                         end % the loop for p 
                end % the loop for n
                ddd(:,:,pol)=dd(:,2:4,pol);
            end % the loop for pol = 
            dmax=d;
            for p1=1:dmax-1
                    for p2=1:dmax-1
                            same_ens(p1,p2)=sum(ddd(p1,:,1)==ddd(p2,:,2));
                            if same_ens(p1,p2)==3 % all three elements identical
                                    for p3=1:dmax-1
                                            samep1p3(p3)=sum(ddd(p1,:,1)==ddd(p3,:,3));
                                            if samep1p3(p3)==3
                                                   enssig_temp=ddd(p1,:,1);
                                                   % creating the automorphisms
                                                   for k=1:3
                                                         stemp(k,:)=mod( [dd(p1,:,1)-dd(p1,k+1,1)] , 121);
                                                         stemp(k,:)=sort(stemp(k,:));
                                                         enssig_temp=[enssig_temp;stemp(k,2:4)];
                                                   end % of creating automorphism
                                            end
                                    end
                            end
                    end
            end
            if length(enssig_temp)<4
                    disp(' No results ')
                    enssig=[ ];
                    break
            end
            enssig=[enssig;enssig_temp];
            dline=reshape(enssig_temp,1,12);
            enssig_temp=[ ]; 
            flag=0;
            for q=1:12
                    el=find(r==dline(q));
                    [s1 s2]=size(el);
                    if s2==1 % implying not an empty matrix
                            elim(q)=el;
                    else
                        elim(q)=0;
                        flag=1;  % implying that at least one element was already used
                    end
            end
            if flag==0
                     r(elim)=[ ]; % removing the elements used now
                     lr=length(r); % updating the length of r
            end  % the elements used were taken out of the vector r
            if flag==1
                disp(' No results ')
                        break
            end
    end % of while lr>11
    if length(enssig)==40
         disp(' ')
        disp(' Resulted array ')
        disp(' ')
        disp(enssig) % displays the entire signal
        flag_result=flag_result+1;
    end
end  % of while flag_result==results-1