How can I feed MATLAB's integral2 a vector of parameters via nested functions?

But technically, there are no matrices involved

There are. The documentation of integral2 specifies that the function being integrated

must accept two arrays of the same size and return an array of corresponding values. It must perform element-wise operations.

That is, integral2 passes some arrays of equal size as alpha and phi into your function lattice_c. This immediately breaks sin(alpha)*cos(phi) - you need .* there.

But that's not the end of story because you also want to multiply that by x. Your idea of introducing position argument and then passing : into it doesn't seem to help (at least in my R2013a, the effect of : as function argument is neither documented nor clear to me).

Without integration, you can check that lattice_b(-1,0:1,0:1) throws an error: the function cannot handle array arguments.

The underlying reason for most of these troubles is that you are stuffing too much logic into anonymous functions. Here is my rewrite: I split off lattice_b, folding lattice_a into it. Due to the presence of a non-anonymous function, I wrapped the top level code into function main.

Oh, and the entries of q should be dealt with one at a time, for which a for loop is perfectly appropriate. There is no sense in trying to vectorize sophisticated operations such as evaluation of double integrals (I explained why elsewhere).

function main()
clear all
L=1000;
R=80;
formfactor = @(q,alpha) 2*sin(q.*cos(alpha)*L/2)./(q.*cos(alpha)*L/2).*besselj(1,q.*sin(alpha)*R)./(q.*sin(alpha)*R);
result = @(q) integral(@(alpha) formfactor(q,alpha).^2.*sin(alpha),0,pi/2,'ArrayValued',true);
lattice_c = @(q,alpha,phi) formfactor(q,alpha).*lattice_b(q,alpha,phi);
lattice_d = @(q) integral2(@(alpha,phi) lattice_c(q,alpha,phi).^2.*sin(alpha),0,pi/2,0,pi/3);

qbins = 100;
q = logspace(-2,0,qbins);
Inew = zeros(size(q));
for i=1:numel(q)
    Inew(i) = lattice_d(q(i));
end
disp(Inew)
end

function z = lattice_b(q,alpha,phi)
a = 24;
b = a*sqrt(3)/2;
x=[-2.5*a, -2*a, -a, -a/2, a/2, a, 2*a, 2.5*a, 2*a, 2.5*a, 2*a, a, a/2, -a/2, -a, -2*a, -2.5*a, -2*a, -a, -a/2, a/2, a, a/2, -a/2 ];
y=[b ,2*b ,2*b ,3*b ,3*b ,2*b ,2*b ,b ,0 ,-b ,-2*b ,-2*b ,-3*b ,-3*b ,-2*b ,-2*b ,-b ,0 ,0 ,b ,b ,0 ,-b ,-b ];
z = 0;
for i=1:numel(x)
    z = z + exp(sqrt(-1)*q*(x(i)*sin(alpha).*cos(phi) + y(i)*sin(alpha).*sin(phi)));
end
end