Create computeGf.m

aligurbu · Dec 21, 2022 · 0d78898 · 0d78898
1 parent bf0d3b9
commit 0d78898
Showing 1 changed file with 32 additions and 0 deletions.
diff --git a/Solver/computeGf.m b/Solver/computeGf.m
@@ -0,0 +1,32 @@
+function cGf = computeGf(cf, cxi, CHI, RotationMatrix, ...
+                         eta, wg, mask_a, mask_b, mu, ...
+                         nlat, nlon, N, NGSphere, NGphi)
+%% Computation of integral with Gf
+%%
+BI = zeros(nlat*nlon,3);
+% parfor m = 1:nlat*nlon
+% parfor works but just little bit faster with a lot of cpu usage
+for m = 1:nlat*nlon
+    chi = reshape(CHI(:,m),1,1,3);
+    Rot = RotationMatrix{m};
+    %% Sample points
+    xihat = InvTransformUpSampleRotateSHCoeff(cxi, Rot, ...
+                                              mask_a, mask_b, ...
+                                              N, NGSphere);
+    %% Density
+    fhat = InvTransformUpSampleRotateSHCoeff(cf, Rot, ...
+                                             mask_a, mask_b, ...
+                                             N, NGSphere);
+    %%
+    rr = chi - xihat;
+    normrr = sqrt(sum(rr.^2,3));
+    rrhat = rr./normrr;
+    BI_ = sum(sum((fhat + rrhat.*sum((rrhat.*fhat),3)).* ...
+                  (eta./normrr).*wg*(2*pi/NGphi)));
+    BI(m,:) = reshape(BI_,1,3);
+end
+%%
+BI = reshape(BI,nlat,nlon,3);
+[aBI,bBI] = shagcm(BI);
+cGf = [aBI(mask_a); bBI(mask_b)]/(8*pi*mu);
+end