function c_1minusalpha = find_crit_RSW_20Dec2020(n,mbar,Ztilde1,Ztilde2,D,W_order_B,B,Omega_hat,delta_1,gamma,p_eta,lambda,alpha,alpha1)

%Create vectors for storing
mbar_star = zeros(size(mbar,2),B); %these would be m x j x n if Gamma_LRR(beta) was non-singleton. (See Sep 25 code)
Qhat_star = zeros(B,1);
T_LRR_star = zeros(B,1);

for b = 1:B %for each bootstrap
    
%Redrawn covariates
Db = D(W_order_B(:,b));

%Indicator functions for D.
D_eq_1 = (Db==1);
D_eq_0 = (Db==0);

%Redrawn Ztilde1, Ztilde2
Z1t = Ztilde1(W_order_B(:,b));
Z2t = Ztilde2(W_order_B(:,b));
  
 %Moment conditions for those parameters
m1 = (Z1t.*D_eq_1 - (p_eta*(delta_1+gamma)+ (1-p_eta)*(log(delta_1 + gamma))).*D_eq_1);
m2 = (Z1t.*D_eq_0 - (p_eta*delta_1 + (1-p_eta)*log(delta_1)).*D_eq_0);
m3 = (p_eta*(delta_1+gamma) + (1-p_eta)*log(delta_1+gamma)).*D_eq_1-Z2t.*D_eq_1;
m4 = (p_eta*delta_1 + (1-p_eta)*log(delta_1)).*D_eq_0-Z2t.*D_eq_0;
  
%mbar_star average bootstrap moments
mbar_star(:,b) = [mean(m1),mean(m2),mean(m3),mean(m4)];

%Test Statistic
Qhat_star(b) = sum(max(((mbar_star(:,b)'-mbar+lambda)./sqrt(Omega_hat)),0));
T_LRR_star(b) = sqrt(n)*Qhat_star(b); %trivial for this case, unique minimizer.

end

% Critical Value is 1-alpha percentile of statistic
T_sorted = sort(T_LRR_star);
c_1minusalpha = T_sorted(ceil((1-alpha+alpha1)*B)); %function of beta value.
end
