function D = match_Ds(seq, genD, V_end, J_start, scorematrix, D_maximum_deletion, negative_excess_deletions_max, max_aligns)
% Given end of V alignment and beginning of J alignment, this function
% matches genomic D sequences to the segment in the middle. Alignments use
% scorematrix.

lseq = length(seq);
flank_length = 10;
seqNDN=seq( V_end +1 - flank_length : min( length(seq), J_start - 1 + flank_length) ); % segment in the middle with flanks.

flank_offset = V_end - flank_length;
n_D_alleles = length(genD);
D.align_length = zeros(max_aligns,n_D_alleles);
D.score = zeros(max_aligns,n_D_alleles);
D.n_errors = zeros(max_aligns,n_D_alleles);
%maximum no. of possible errors with this scoring matrix. assuming maximum genD length is 16.
D_max_errors = max(1,floor(16*scorematrix(1,1)/(scorematrix(1,1) -scorematrix(1,2) )));

D.error_positions = zeros(max_aligns, D_max_errors , n_D_alleles); % positions of mismatches
D.excess_error_positions_left = zeros(max_aligns,negative_excess_deletions_max, n_D_alleles);
D.excess_error_positions_right = zeros(max_aligns,negative_excess_deletions_max, n_D_alleles);

D.align_position_left = zeros(max_aligns,n_D_alleles);
D.align_position_right = zeros(max_aligns,n_D_alleles);
D.deletions_left = zeros(max_aligns,n_D_alleles);
D.deletions_right = zeros(max_aligns, n_D_alleles);
D.p_region_max_length_left = zeros(max_aligns, D_maximum_deletion + 1, n_D_alleles);
D.p_region_max_length_right = zeros(max_aligns, D_maximum_deletion + 1, n_D_alleles);
D.n_alignments = zeros(n_D_alleles,1);

% Get ALL local alignments of length > 0.
if ~isempty(seqNDN)
    for d=1:length(genD)
        A = localalign(seqNDN,genD(d).seq,'Alphabet','nt','ScoringMatrix',scorematrix,'GapOpen',1000000,'NUMALN',max_aligns,'DOALIGNMENT',false);
        n_aligns_plus_one = find(A.Stop(:,1)==0,1);
        if isempty(n_aligns_plus_one)
            n_aligns = max_aligns;
        else
            n_aligns = n_aligns_plus_one -1;
        end
        
        D.align_length(1:n_aligns,d) = A.Stop(1:n_aligns,1) - A.Start(1:n_aligns,1) + 1;
        D.score(1:n_aligns,d) = A.Score(1:n_aligns);
        D.n_errors(1:n_aligns,d) =  round((D.score(1:n_aligns,d) -D.align_length(1:n_aligns,d)*scorematrix(1,1) )/(scorematrix(1,2) - scorematrix(1,1)));
        
        D.align_position_left(1:n_aligns,d) = A.Start(1:n_aligns,1) + flank_offset;
        D.align_position_right(1:n_aligns,d) = A.Stop(1:n_aligns,1) + flank_offset;
        D.deletions_left(1:n_aligns,d) = A.Start(1:n_aligns,2)-1;
        D.deletions_right(1:n_aligns,d) = length(genD(d).seq) - A.Stop(1:n_aligns,2);
        D.n_alignments(d) = n_aligns;
    end
end

[~,D.allele_order] = sort(D.score(1,:),'descend');
%Now go through and find palindromic nucleotides for each alignment
for d=1:length(genD)
    if D.align_length(1,d)>0
        % First on the left
        
        % Loop over number of deletions
        %for nd = max(0,D.min_deletions_left(d) - negative_excess_deletions_max):min(D_maximum_deletion, D.align_length(d) + D.min_deletions_left(d))
        n_aligns = D.n_alignments(d);
        for na=1:n_aligns
            
            if D.n_errors(na,d)>0
                D.error_positions(na,1:D.n_errors(na,d),d) = D.align_position_left(na,d) - 1 +  find(seq(D.align_position_left(na,d):D.align_position_right(na,d))~=genD(d).seq(D.deletions_left(na,d)+1: (end-D.deletions_right(na,d)) ));
            end
            
            % For each value of deletions, find longest half-palindrome
            % from the implied end of the gene sequence            
            for nd=(max(0,D.deletions_left(na,d) - negative_excess_deletions_max)):min(D_maximum_deletion, D.align_length(na,d) + D.deletions_left(na,d))
                p=0;
                still_palindrome=true;
                while still_palindrome==1 && p < min((D.align_length(na,d)-(nd-D.deletions_left(na,d))), D.align_position_left(na,d)-1 + nd - D.deletions_left(na,d)  )
                    still_palindrome = genD(d).seq( nd + 1 + p) == dnacomplement(seq(D.align_position_left(na,d) - p + nd - D.deletions_left(na,d) - 1));
                    if still_palindrome
                        p = p + 1;
                    end
                end
                D.p_region_max_length_left(na, 1 + nd, d) = p;
                
            end
            
            
            % Now for the 'negative' deletions, store where the mismatches are.
            % This is so we can count number of errors for these possibilities.
            
            % Number of 'negative' deletions to consider:
            n_excess_left = min([negative_excess_deletions_max, D.deletions_left(na,d), D.align_position_left(na,d)-1 ]);
            
            % Find mismatches between sequence and genomic J at those positions.
            excess_err_pos_left = D.align_position_left(na,d) - ( n_excess_left + 1 - find((genD(d).seq((D.deletions_left(na,d)-n_excess_left+1): D.deletions_left(na,d) )~=seq((D.align_position_left(na,d) -n_excess_left): (D.align_position_left(na,d)-1) )) ) );
            if ~isempty(excess_err_pos_left)
                D.excess_error_positions_left(na, 1:length(excess_err_pos_left),d) = excess_err_pos_left;
            end
            
            for nd= (max(0,D.deletions_right(na,d) - negative_excess_deletions_max)):min(D_maximum_deletion, D.align_length(na,d) + D.deletions_right(na,d)  )
                
                % right palindrome
                p=0;
                still_palindrome=true;
                while still_palindrome==1 && p < min(D.align_length(na,d) - (nd - D.deletions_right(na,d)), lseq - D.align_position_right(na,d) + nd - D.deletions_right(na,d) )
                    still_palindrome = genD(d).seq( end - nd - p) == dnacomplement(seq(D.align_position_right(na,d)+p -(nd - D.deletions_right(na,d)) + 1));
                    if still_palindrome
                        p = p + 1;
                    end
                end
                D.p_region_max_length_right(na,1 + nd,d) = p;
                
            end
            
            % Now for the 'negative' deletions, store where the mismatches are.
            % This is so we can count number of errors for these possibilities.
            
            % Number of 'negative' deletions to consider:
            n_excess_right = min([negative_excess_deletions_max, D.deletions_right(na,d), length(seq) - D.align_position_right(na,d) ]);
            % Find mismatches between sequence and genomic V at those positions.
            excess_err_pos_right = D.align_position_right(na,d) + find((genD(d).seq((end-D.deletions_right(na,d)+1): (end - D.deletions_right(na,d) + n_excess_right)  )~=seq( (D.align_position_right(na,d)+1): (D.align_position_right(na,d) + n_excess_right) )));
            if ~isempty(excess_err_pos_right)
                D.excess_error_positions_right(na,1:length(excess_err_pos_right),d) = excess_err_pos_right;
            end
            
            
            
        end
    end
end

end
