
function [XW,domain_struct]=cub_QMC(domain_struct,card)

% Note: "domain_struct" in input has been enforced by adding 
%       informations about bounding box and area.

% compute area and bounding box (stored in "domain_struct").
deg_area=0;
if strcmp(domain_struct.domain,'butterfly'), deg_area=1; end

[~,~,domain_struct]=define_cub_rule(domain_struct,deg_area);

% area(bounding box) and area(domain)
area_domain=domain_struct.area;
dbox=domain_struct.dbox;

xLimit=dbox(1,:); yLimit=dbox(2,:);
area_dbox=diff(xLimit)*diff(yLimit);


if area_domain > area_dbox
    fprintf(1,'\n \t area_domain: %1.5e',area_domain);
    fprintf(1,'\n \t area_dbox  : %1.5e \n',area_dbox);
    warning('domain area larger than bounding box area');
end

% ratio: area(bounding box)/area(domain)
ratio=area_dbox/area_domain;

% computing sufficiently large pointset on [0,1] x [0,1]
P=haltonset(2);
card_bbox=ceil(4*ratio*card);
t0_ref = net(P,card_bbox);

% mapping sufficiently large pointset on domain bounding box
t0=[xLimit(1)+(t0_ref(:,1))*diff(xLimit) ...
    yLimit(1)+(t0_ref(:,2))*diff(yLimit)];

if size(t0,1) < card
    warning('few points in the domain');
end

% extracting points on domain
in0=indomain_routine(domain_struct,t0);
in=find(in0 == 1);
card_dom=length(in);
t=t0(in(1:card),:);

% weights
w=area_dbox*(card_dom/card_bbox)/card*ones(size(t));

% setting output variables
XW=[t w]; 
