Каков наименьший по времени способ решения сотен нелинейных уравнений?

Question

Я хочу решить нелинейное уравнение следующим образом: EQ

параметр "ta" является неизвестным, а параметры "Ps" и "PTOA" являются матрицами m-by-n, а другие параметры известны. На самом деле, я должен решить это уравнение для каждого элемента в матрицах "Ps" и "PTOA", чтобы в итоге получить матрицу "ta", имеющую размерность m-на-n. Мой код занимает много времени для расчета "та" матрицы. Кто-нибудь может предложить лучший способ сделать это? вот мой код:

clear; clc; close all;

% % % % sza: sun zenith angle
% % % % u_s: cosine of sun zenith angle
% % % % u_v: cosine of view/sensor zenith angle
% % % % w0: aerosol single scattering albedo(ssa)
% % % % wR: rayleigh single scattering albedo
% % % % theta: scattering angle
% % % % Pa: aerosol scattering phase function
% % % % PR: rayleigh scattering phase function
% % % % R_TOA: top of atmosphere reflectance at specified wavelength
% % % % R_RAY: rayleigh reflectance
% % % % tR: rayleigh optical depth
% % % % ta: aerosol optical thickness
% % % % R_sur: surface reflectance
% % % % g: asymmetry factor
% % % % lambda: wavelength

wave = xlsread('E:\PayanNameh\SSA_AsymmetryFactor\aerosolModels.xlsx','6S urban', 'A3:A13');
ssa = xlsread('E:\PayanNameh\SSA_AsymmetryFactor\aerosolModels.xlsx','6S urban', 'B3:B13');
g = xlsread('E:\PayanNameh\SSA_AsymmetryFactor\aerosolModels.xlsx','6S urban', 'C3:C13');
% % % % g_660 = interp1(wave, g, 0.660);

lambda = [0.480];
sza = input('Enter sun zenith angle(degree): ');
u_s = cosd(sza);
u_v = 1;
w0 = interp1(wave, ssa, lambda);
g = interp1(wave, g, lambda);
Pa = (1-g^2)/(1+g^2-2*g*(u_s*u_v))^1.5;
PR = 0.75*(1+(u_s*u_v)^2);
tR = 0.00864/lambda^(3.916+0.074*lambda+0.05/lambda);
R_RAY = tR*PR/(4*u_s*u_v);
R_TOA = imread('E:\PayanNameh\Data\landsat8_toaReflectance_Blue_20140122_60.tif'); R_TOA = double(R_TOA);
R_sur = imread('E:\PayanNameh\Data\landsat8_srReflectance_Blue_20140122_60.tif'); R_sur = double(double(R_sur)/10000);
R_sur(R_sur<0) = 0;
R_TOA(isnan(R_TOA)) = 0;
% R_TOA1 = imread('E:\PayanNameh\Data\landsat8_toaReflectance_Blue_20140122_60.tif'); R_TOA1 = double(R_TOA1);
% R_sur1 = imread('E:\PayanNameh\Data\landsat8_srReflectance_Blue_20140122_60.tif'); R_sur1 = double(R_sur1)/10000;
% R_sur1(R_sur1<0) = 0;
% R_TOA = R_TOA1(151:162,599:610); R_sur = R_sur1(151:162,599:610); R_TOA(isnan(R_TOA)) = 0;


% aot = double(vpasolve(ta == 4*u_s*u_v*(R_TOA-R_sur-(exp(-(tR+ta)/u_s)*exp(-(tR+ta)))*R_sur/(1-R_sur*(0.92*tR+(1-g)*ta)*exp(-(tR+ta))))/(w0*Pa), ta));


syms ta
aot = zeros(size(R_TOA,1),size(R_TOA,2));
tic
for i = 1:size(R_TOA,1)
    for j = 1:size(R_TOA,2)
        if R_TOA(i,j) ~= 0
            bb = double(vpasolve(ta == 4*u_s*u_v*(R_TOA(i,j)-R_RAY-(exp(-(tR+ta)/u_s)*exp(-(tR+ta)))*R_sur(i,j)/(1-R_sur(i,j)*(0.92*tR+(1-g)*ta)*exp(-(tR+ta))))/(w0*Pa), ta));
            if size(bb,1) == 0
                bb = NaN;
                aot(i,j) = bb;
            else
                aot(i,j) = bb;
            end
        end
    end
end
toc
aot(aot == 0) = NaN;

Есть ли кто-нибудь, кто мог бы написать версию этого кода на javascript?!