function [Yf, details] = Yf_from_wq(W_cell, Q_cell, alphaJ, alphaO, alphaf, params)
% YF_FROM_WQ  Compute firm output Y_f from task wages/prices via l_j, k_j, y_j, Y_{o,f}, then Y_f.
%
% Definitions used:
%   p_j = [ a*(M w)^(1-ν) + b*q^(1-ν) ]^(1/(1-ν))
%   X_{o,f} = [ Σ_j α_j p_j^(1-ψ) ]^(-1/(1-ψ)),   P_{o,f} = 1/X_{o,f}
%   Z_f = [ Σ_o α_{o,f} P_{o,f}^{1-χ} ]^(-1/(1-χ))
%   ℓ_j = a M^{-ν} w^{1-ν} [ a M^{1-ν} w^{1-ν} + b q^{1-ν} ]^{(ν-ψ)/(1-ν)}  X^{χ-ψ} Z^{θ-χ}
%   k_j = b           q^{1-ν} [ a M^{1-ν} w^{1-ν} + b q^{1-ν} ]^{(ν-ψ)/(1-ν)}  X^{χ-ψ} Z^{θ-χ}
%   y_j = [ γ_j ℓ_j^{(ν-1)/ν} + (1-γ_j) k_j^{(ν-1)/ν} ]^{ν/(ν-1)}
%   Y_{o,f} = [ Σ_j α_j^{1/ψ} y_j^{(ψ-1)/ψ} ]^{ψ/(ψ-1)}      (GM limit if ψ→1)
%   Y_f     = [ Σ_o α_{o,f}^{1/χ} Y_{o,f}^{(χ-1)/χ} ]^{χ/(χ-1)} (GM limit if χ→1)

    % ---- Params & guards ----
    a     = getd(params,'a',     []);
    b     = getd(params,'b',     []);
    nu    = getd(params,'nu',    []);
    psi   = getd(params,'psi',   []);
    chi   = getd(params,'chi',   []);
    theta = getd(params,'theta', []); %#ok<NASGU> % not used directly here
    M     = getd(params,'M',     1.0);

    if isempty(a) || isempty(b) || isempty(nu) || isempty(psi) || isempty(chi)
        error('params must include a,b,nu,psi,chi (and optionally M,theta).');
    end
    if numel(W_cell) ~= numel(Q_cell) || numel(W_cell) ~= numel(alphaJ)
        error('W_cell, Q_cell, alphaJ must have same length O.');
    end

    tol = 1e-12;
    O   = numel(W_cell);

    % Normalize α_o within firm
    alphaO = alphaO(:);
    if sum(alphaO)<=0, error('alphaO must sum positive.'); end
    alphaO = alphaO / sum(alphaO);

    % Optional γ_j per task; else infer scalar γ from a=γ^ν, b=(1-γ)^ν
    if isfield(params,'gammaJ') && ~isempty(params.gammaJ)
        gammaJ = params.gammaJ;  % cell {O}, each [J_o x 1]
        if ~iscell(gammaJ) || numel(gammaJ) ~= O
            error('params.gammaJ must be a cell {O} matching alphaJ.');
        end
    else
        gamma = a^(1/nu);  % ensures a=γ^ν, b=(1-γ)^ν
        gammaJ = cell(O,1);
        for o=1:O
            gammaJ{o} = gamma * ones(numel(alphaJ{o}),1);
        end
    end

    % ---- First pass: p_j, X_{o,f}, P_{o,f} ----
    pj_cell = cell(O,1);
    Xof     = zeros(O,1);
    Pof     = zeros(O,1);

    for o = 1:O
        w = W_cell{o}(:);
        q = Q_cell{o}(:);
        aj = alphaJ{o}(:);
        if sum(aj)<=0, aj = ones(size(aj)); end
        aj = aj / sum(aj);  % normalize α_j within occ

        % p_j
        t1 = a .* ((M.*w).^(1 - nu));
        t2 = b .* (q.^(1 - nu));
        pj = (t1 + t2).^( 1/(1 - nu) );
        pj_cell{o} = pj;

        % X_{o,f}, GM limit if ψ≈1
        if abs(psi-1) > tol
            s = sum( aj .* (pj.^(1 - psi)) );
            Xof(o) = s.^(-1/(1 - psi));
        else
            Xof(o) = exp( -sum( aj .* log(pj) ) );
        end
        Pof(o) = 1 ./ Xof(o);
    end

    % ---- Firm price aggregator: Z_f ----
    if abs(chi-1) > tol
        sZ = sum( alphaO(:) .* (Pof(:).^(1 - chi)) );
        Zf = sZ.^(-1/(1 - chi));
    else
        wts = alphaO / sum(alphaO);
        Zf  = exp( -sum( wts .* log(Pof) ) );
    end

    % ---- Second pass: l_j, k_j, y_j, then Y_{o,f} ----
    lj_cell = cell(O,1);
    kj_cell = cell(O,1);
    yj_cell = cell(O,1);
    Yof     = zeros(O,1);

    for o = 1:O
        w  = W_cell{o}(:);
        q  = Q_cell{o}(:);
        aj = alphaJ{o}(:);  aj = aj / sum(aj);
        gj = gammaJ{o}(:);
        X  = Xof(o);

        % common core
        core = a.*(M.^(1 - nu)).*(w.^(1 - nu)) + b.*(q.^(1 - nu));
        scale = (core.^((nu - psi)/(1 - nu))) .* (X.^(chi - psi)) .* (Zf.^(0 + (theta - chi))); %#ok<*NASGU>

        % l_j and k_j (note: NO M^{-nu} on k_j; correction applied)
        lj = alphaf*alphaO(o)*aj.*a .* (M.^(-nu)) .* (w.^(1 - nu)) .* scale;
        kj = alphaf*alphaO(o)*aj.*b .*              (q.^(1 - nu)) .* scale;

        % y_j CES between l and k with elasticity ν and weights γ_j
        if abs(nu-1) > tol
            yj = ( gj .* (lj.^((nu - 1)/nu)) + (1 - gj) .* (kj.^((nu - 1)/nu)) ).^( nu/(nu - 1) );
        else
            error('nu must not equal 1.');
        end

        lj_cell{o} = lj;
        kj_cell{o} = kj;
        yj_cell{o} = yj;

        % Y_{o,f} (GM limit if ψ≈1)
        if abs(psi-1) > tol
            s = sum( (aj.^(1/psi)) .* ( yj.^((psi - 1)/psi) ) );
            Yof(o) = s.^( psi/(psi - 1) );
        else
            wts = aj; % normalized already
            Yof(o) = exp( sum( wts .* log(yj) ) );
        end
    end

    % ---- Firm output Y_f from Y_{o,f} (GM limit if χ≈1) ----
    if abs(chi-1) > tol
        sY = sum( (alphaO.^(1/chi)) .* ( Yof.^((chi - 1)/chi) ) );
        Yf = sY.^( chi/(chi - 1) );
    else
        wts = alphaO;
        Yf = exp( sum( wts .* log(Yof) ) );
    end

    % ---- Pack details ----
    details = struct();
    details.pj_cell = pj_cell;
    details.Xof     = Xof;
    details.Pof     = Pof;
    details.Zf      = Zf;
    details.lj_cell = lj_cell;
    details.kj_cell = kj_cell;
    details.yj_cell = yj_cell;
    details.Yof     = Yof;
end

% --- tiny utility ---
function v = getd(s,k,d)
    if isfield(s,k) && ~isempty(s.(k)), v = s.(k); else, v = d; end
end
