Create computeK.m

Pre-compute the integrand with K kernel

Solver/computeK.m

@@ -0,0 +1,42 @@
+function Kernel = computeK(cxi, CHI, RotationMatrix, ...
+                           eta, wg, mask_a, mask_b, ...
+                           nlat, nlon, N, NGSphere, NGthet, NGphi)
+%% Pre-compute the integrand with K kernel
+Kernel = cell(nlat*nlon,1);
+
+% parfor m = 1:nlat*nlon 
+% parfor works and it is faster with 100% cpu usage
+for m = 1:nlat*nlon
+    KK = zeros(NGthet, NGphi,6);
+    Normal = zeros(NGthet, NGphi,3);
+    chi = reshape(CHI(:,m),1,1,3);
+    Rot = RotationMatrix{m};
+
+    %% Sample points 
+    [Xihat, ahatXi_, bhatXi_] = ...
+                      InvTransformUpSampleRotateSHCoeff(cxi, Rot, ...
+                                                        mask_a, mask_b, ...
+                                                        N, NGSphere);
+    %%
+    [gXihat_thet,gXihat_phi] = gradgcm(ahatXi_,bhatXi_);
+    rr = chi - Xihat;
+    normrr = sqrt(sum(rr.^2,3));
+    rrhat = rr./normrr;
+
+    Normal(:,:,1) = gXihat_thet(:,:,2).*gXihat_phi(:,:,3) - ...
+                    gXihat_thet(:,:,3).*gXihat_phi(:,:,2);
+    Normal(:,:,2) = gXihat_thet(:,:,3).*gXihat_phi(:,:,1) - ...
+                    gXihat_thet(:,:,1).*gXihat_phi(:,:,3);
+    Normal(:,:,3) = gXihat_thet(:,:,1).*gXihat_phi(:,:,2) - ...
+                    gXihat_thet(:,:,2).*gXihat_phi(:,:,1);
+
+    KK(:,:,1) = rrhat(:,:,1).*rrhat(:,:,1);
+    KK(:,:,2) = rrhat(:,:,2).*rrhat(:,:,2);
+    KK(:,:,3) = rrhat(:,:,3).*rrhat(:,:,3);
+    KK(:,:,4) = rrhat(:,:,1).*rrhat(:,:,2);
+    KK(:,:,5) = rrhat(:,:,2).*rrhat(:,:,3);
+    KK(:,:,6) = rrhat(:,:,1).*rrhat(:,:,3);
+    Kernel{m} = 6*KK.*(sum(rrhat.*Normal,3)./normrr) ...
+                    .*(eta./normrr).*wg*(2*pi/NGphi);
+end
+end