cgcoeff_stage2( J, init_rhs_i, init_rhs_c )

function [ lhs, cell_rhs_i, cell_rhs_c ] = cgcoeff_stage2( J, init_rhs_i, init_rhs_c )
%CGCOEFF_STAGE2 for a given J, it computes RHS of |J, M> for M=J:-1:-J
%in
%   J: J of lhs
%   init_rhs_i: rhs index matix, nx4 (j1, m1, j2, m2)
%   init_rhs_c: rhs coeff matrix
%out
%   cell_rhs_i: all the RHS indices (|j1, m1; j2, m2>)
%   cell_rhs_c: all the RHS coefficients

cell_rhs_i = {init_rhs_i}; %set of 4 (j1, m1, j2, m2), %ij1 = 1, im1 = 2, ij2 = 3, im2 = 4
cell_rhs_c = {sym(init_rhs_c)}; %const matrix

Ms = (J:-1:-J);
for r = 2:numel(Ms)
    prevM = Ms(r-1);
    prev_rhs_i = cell_rhs_i{r-1};
    prev_rhs_c = cell_rhs_c{r-1};
    [numprevtouble, temp] = size(prev_rhs_i); %temp will be 4
    curr_rhs_i = [];
    curr_rhs_c = [];
    for c=1:numprevtouble
        [children_i children_c] = op_jminus(prev_rhs_i(c,:), prev_rhs_c(c));
        for ichild=1:numel(children_c)
            achild_i = children_i(ichild,:);
            achild_c = children_c(ichild);
            iloc = locateindex(curr_rhs_i, achild_i);
            if (iloc == 0)
                curr_rhs_i(end+1,:) = achild_i;
                curr_rhs_c(end+1) = achild_c;
            else
                curr_rhs_c(iloc) = curr_rhs_c(iloc) + achild_c;
            end
        end
    end
    curr_rhs_c = curr_rhs_c / coef_jminus(J, prevM);
    curr_rhs_c = simple(curr_rhs_c);
    cell_rhs_i = {cell_rhs_i{:} curr_rhs_i};
    cell_rhs_c = {cell_rhs_c{:} curr_rhs_c};
    r = r + 1;
end

%build lhs
lhs =[];
for i = 1:numel(Ms)
    lhs = [lhs ; [J, Ms(i)]];
end

%display
% for i = 1:numel(Ms)
%     disp(['J=' num2str(J), ',M=' num2str(Ms(i))]);
%     [cell_rhs_i{i}, cell_rhs_c{i}']
% end