将上述代码向量化

以下是对上述代码的向量化改进：

% 预处理一些固定的变量
z = (0:nj-2)' * dz;
a = (0:2:(2*nj-4))' * dz * 0.5;
b = (2:2:(2*nj-2))' * dz * 0.5;
rydis = (a + b) * 0.5;
jsdis = repelem(z, nj-2);

% 计算 Tj
Tj = zeros(LL, nj, nj, zx*zy);
JSX = 0; JSY = 0;
Xt = abs(JSX - (1:zx)' * gj);
rbs = ((1:zx)' - a1).^2 * gj^(-2) + 1;
parfor r = 1:zx*zy
    ZKX = (ceil(r/zy) - 1) * gj;
    ZKY = mod(r-1, zy) * gj;
    rr = sqrt((ZKX-JSX).^2 + (ZKY-JSY).^2);
    if rr == 0 % 避免除以 0 出错
        Tj(:, :, :, r) = 0;
        continue;
    end
    for k = 1:nj
        x = z(k);
        ax = repmat(a(k), nj-2, 1);
        bx = repmat(b(k), nj-2, 1);
        rydisx = repmat(rydis(k), nj-2, 1);
        jsdisx = jsdis;
        [v, Rap, Iap, Rlamd, Ilamd] = untitled55(rydis(k), jsdis(k));
        aa = trapz(bx - ax, ...
            0.25 * exp(v*Xt*0.5/Rap) .* exp(-v*sqrt(rr^2 + (z(k) - jsdisx).^2)*0.5/Rap) .* ...
            erfc((sqrt(rr^2 + (z(k) - jsdisx).^2) - v*(0:LL-1)'*dt)*0.5/sqrt(Rap*(0:LL-1)'*dt)) ./ ...
            sqrt(rr^2 + (z(k) - jsdisx).^2) ./ (2*pi*Rlamd), 1);
        ab = trapz(bx - ax, ...
            0.25 * exp(v*Xt*0.5/Rap) .* exp(v*sqrt(rr^2 + (z(k) - jsdisx).^2)*0.5/Rap) .* ...
            erfc((sqrt(rr^2 + (z(k) - jsdisx).^2) + v*(0:LL-1)'*dt)*0.5/sqrt(Rap*(0:LL-1)'*dt)) ./ ...
            sqrt(rr^2 + (z(k) - jsdisx).^2) ./ (2*pi*Rlamd), 1);
        ac = trapz(bx - ax, ...
            0.25 * exp(v*Xt*0.5/Iap) .* exp(-v*sqrt(rr^2 + (z(k) + jsdisx).^2)*0.5/Iap) .* ...
            erfc((sqrt(rr^2 + (z(k) + jsdisx).^2) - v*(0:LL-1)'*dt)*0.5/sqrt(Iap*(0:LL-1)'*dt)) ./ ...
            sqrt(rr^2 + (z(k) + jsdisx).^2) ./ (2*pi*Ilamd), 1);
        ad = trapz(bx - ax, ...
            0.25 * exp(v*Xt*0.5/Iap) .* exp(v*sqrt(rr^2 + (z(k) + jsdisx).^2)*0.5/Iap) .* ...
            erfc((sqrt(rr^2 + (z(k) + jsdisx).^2) + v*(0:LL-1)'*dt)*0.5/sqrt(Iap*(0:LL-1)'*dt)) ./ ...
            sqrt(rr^2 + (z(k) + jsdisx).^2) ./ (2*pi*Ilamd), 1);
        if (a1==2 && a2==1)
            aa(:) = 0;
            ab(:) = 0;
            ac(:) = 0;
            ad(:) = 0;
        end
        Tj(:, 1, k, r) = aa + ab - ac - ad;
    end
end

这里我们使用了 `repmat` 函数和 `trapz` 函数，将循环中的一些计算向量化。此外，我们在循环外部预处理了一些固定的变量，避免在循环内部重复计算，提高了计算效率。