SolveSpectralBIE.m

function [u, au, bu, cu, FLAG, RELRES, ITER] = ...
                           SolveSpectralBIE(axi, bxi, xi, ...
                                            RotationMatrix, ...
                                            NGSphere, NGthet, NGphi, ...
                                            eta, wg, ...
                                            af, bf, auinf, buinf, ...
                                            nlat, nlon, mask_a, mask_b, ...
                                            mu, lam, N, ToleranceGMRES, ...
                                            cuprev)
cxi = [axi(mask_a); bxi(mask_b)];
CHI = reshape(xi,nlat*nlon,3)';
%% Compute uinf - Gf
cf = [af(mask_a); bf(mask_b)];
cGf = computeGf(cf, cxi, CHI, RotationMatrix, ...
                eta, wg, mask_a, mask_b, mu, ...
                nlat, nlon, N, NGSphere, NGphi);
cuinf = [auinf(mask_a); buinf(mask_b)];
rhs = cuinf - cGf;

%% Compute u + (lam-1)/8/pi*Ku
Kernel = computeK(cxi, CHI, RotationMatrix, ...
                  eta, wg, mask_a, mask_b, ...
                  nlat, nlon, N, NGSphere, NGthet, NGphi);

%% Solve spectal Galerkin BIE 
[cu, FLAG, RELRES, ITER] = gmres(@computeSpectralBIE,rhs,[], ...
                                 ToleranceGMRES,3*(N+1)^2,[],[],cuprev);

au = zeros(size(axi));  bu = zeros(size(bxi));
au(mask_a) = cu(1:3*(N+1)*(N+2)/2);
bu(mask_b) = cu(3*(N+1)*(N+2)/2+1:3*(N+1)^2);
u = shsgcm(au, bu);
    %%
    function c = computeSpectralBIE(cu)
        %%
        au = zeros(size(axi));  bu = zeros(size(bxi));
        au(mask_a) = cu(1:3*(N+1)*(N+2)/2);
        bu(mask_b) = cu(3*(N+1)*(N+2)/2+1:3*(N+1)^2);
        u = shsgcm(au, bu);
        c = computeKu(Kernel, cu, u, RotationMatrix, ...
                      mask_a, mask_b, lam, N, nlat, nlon, NGSphere);
    end
end