initialize_basis_tables Subroutine

  subroutine initialize_basis_tables(plan, order)
    type(fmm_plan_type), intent(inout) :: plan
    integer(i32), intent(in) :: order
    integer(i32) :: idx, alpha_idx, beta_idx, ax, term_idx
    integer(i32) :: sum_exp(3)

    plan%ncoef = count_multi_indices(order)
    plan%nderiv = count_multi_indices(2_i32*order)
    call build_multi_index_table(order, plan%alpha, plan%alpha_degree, plan%alpha_factorial, plan%alpha_map)
    call build_multi_index_table(2_i32*order, plan%deriv_alpha, plan%deriv_degree, plan%deriv_factorial, plan%deriv_map)
    allocate (plan%alpha_sign(plan%ncoef), plan%alpha_plus_axis(3, plan%ncoef), &
              plan%alpha_beta_deriv_idx(plan%ncoef, plan%ncoef))
    do idx = 1_i32, plan%ncoef
      if (mod(plan%alpha_degree(idx), 2_i32) == 0_i32) then
        plan%alpha_sign(idx) = 1.0d0
      else
        plan%alpha_sign(idx) = -1.0d0
      end if
      do ax = 1_i32, 3_i32
        sum_exp = plan%alpha(:, idx)
        sum_exp(ax) = sum_exp(ax) + 1_i32
        if (sum(sum_exp) <= order) then
          plan%alpha_plus_axis(ax, idx) = plan%alpha_map(sum_exp(1), sum_exp(2), sum_exp(3))
        else
          plan%alpha_plus_axis(ax, idx) = 0_i32
        end if
      end do
    end do

    do alpha_idx = 1_i32, plan%ncoef
      do beta_idx = 1_i32, plan%ncoef
        sum_exp = plan%alpha(:, alpha_idx) + plan%alpha(:, beta_idx)
        plan%alpha_beta_deriv_idx(alpha_idx, beta_idx) = plan%deriv_map(sum_exp(1), sum_exp(2), sum_exp(3))
      end do
    end do

    plan%eval_term_count = count(plan%alpha_degree < order)
    allocate (plan%eval_exp(3, max(1_i32, plan%eval_term_count)))
    allocate (plan%eval_deriv_idx(3, max(1_i32, plan%eval_term_count)))
    allocate (plan%eval_inv_factorial(max(1_i32, plan%eval_term_count)))
    plan%eval_exp = 0_i32
    plan%eval_deriv_idx = 0_i32
    plan%eval_inv_factorial = 0.0d0

    term_idx = 0_i32
    do alpha_idx = 1_i32, plan%ncoef
      if (plan%alpha_degree(alpha_idx) >= order) cycle
      term_idx = term_idx + 1_i32
      plan%eval_exp(:, term_idx) = plan%alpha(:, alpha_idx)
      plan%eval_deriv_idx(:, term_idx) = plan%alpha_plus_axis(:, alpha_idx)
      plan%eval_inv_factorial(term_idx) = 1.0d0/plan%alpha_factorial(alpha_idx)
    end do
  end subroutine initialize_basis_tables