Changeset 25357 for issm/trunk-jpl/src/m/contrib/adhikari/yi_analytic_homogenous.m

Timestamp:

08/07/20 17:58:28 (5 years ago)

Author:

adhikari

Message:

CHG: reworded on analytic solutions for computing love_at_depth for honogeneous planet. Plus a couple of figures that show comparision between ISSM (circles) and analytic (solid lines) for both surface loading and tidal forcing

File:

• issm/trunk-jpl/src/m/contrib/adhikari/yi_analytic_homogenous.m (modified) (1 diff)

issm/trunk-jpl/src/m/contrib/adhikari/yi_analytic_homogenous.m

@@ -1 +1 @@
 
-function [yana r] = love_analytic(n,N,source); 
-        % n = SH degree 
-        % N = number of layers 
-        % source = 9 (volumentric potential) 
-        %                       = 11 (surface loads) 
+function yana = love_analytic(param); 
+        % param = required parameters 
+        % param.rho = earth density [kg/m3]
+        % param.mu = shear modulus [Pa]
+        % param.G = universal gravitational constant [m3 kg-1 s-2] 
+        % param.radius = radius of the surface and internal layers of the planet. [m] 
+        % param.g0 = acceleration due to gravity on the planetary surface. [m s-2]
+        % param.source = model forcing type [9 volumentaric potential, 11 surface loads] 
+        % param.degree = SH degree;
 
-rho=5511; %density
-mu=.75e11; %shear modulus
-G=6.67e-11;
-a=6371e3; %surface radius
-g0=9.8134357285509388; % surface gravity 
+        %dbstop if naninf
 
-r=linspace(10e3,a,N+1); 
+        rho = param.rho; 
+        mu = param.mu; 
+        G = param.G; 
+        a = max(param.radius); 
+        g0 = param.g0; 
+        r = param.radius; 
+        source = param.source; 
+        n = param.degree; 
 
-C1=0;
-C2=0;
-C6=0;
-if (source==9)
-        C3 = -1/2*n*(n+1)*rho/((2*n^2+4*n+3)*mu+n*rho*g0*a)*a^(-n);
-        C4 = 1/2/(n-1)*n^2*(n+2)*rho/((2*n^2+4*n+3)*mu+n*rho*g0*a)*a^(-n+2);
-        C5 = (1+3/2/(n-1)*n*rho*g0*a/((2*n^2+4*n+3)*mu+n*rho*g0*a))*a^(-n);
-elseif (source==11)
-        C3=(n*rho*(n^2 - 1)) / (a^n*(4*G*pi*a^2*n*rho^2 + 3*mu*(2*n^2 + 4*n + 3))); 
-        C4=-(a^2*n^2*rho*(n + 2)) / (a^n*(4*G*pi*a^2*n*rho^2 + 3*mu*(2*n^2 + 4*n + 3))); 
-        C5=(3*mu*(2*n^2 + 4*n + 3)) / (a^n*(4*G*pi*a^2*n*rho^2 + 3*mu*(2*n^2 + 4*n + 3))); 
-        %C3 = rho*(n-1)^2*(n+1)*n / (6*mu*a^n*(n-1)*(2*n^2+4*n+3) - 4*pi*G*rho^2*n*(3*a^(-n+2)-(2*n+1)*a^(n+2))); 
-        %C4 = -a^2*(n+2)*n/(n^2-1) * C3; 
-        %C5 = 1/a^n - 4*pi*G*rho*a^2/(n^2-1) * C3; 
-end
+        C1=0;
+        C2=0;
+        C6=0;
+        if (source==9)
+                C3 = -1/2*n*(n+1)*rho/((2*n^2+4*n+3)*mu+n*rho*g0*a)*a^(-n);
+                C4 = 1/2/(n-1)*n^2*(n+2)*rho/((2*n^2+4*n+3)*mu+n*rho*g0*a)*a^(-n+2);
+                C5 = (1+3/2/(n-1)*n*rho*g0*a/((2*n^2+4*n+3)*mu+n*rho*g0*a))*a^(-n);
+        elseif (source==11)
+                C3=(n*rho*(n^2 - 1)) / (a^n*(4*G*pi*a^2*n*rho^2 + 3*mu*(2*n^2 + 4*n + 3))); 
+                C4=-(a^2*n^2*rho*(n + 2)) / (a^n*(4*G*pi*a^2*n*rho^2 + 3*mu*(2*n^2 + 4*n + 3))); 
+                C5=(3*mu*(2*n^2 + 4*n + 3)) / (a^n*(4*G*pi*a^2*n*rho^2 + 3*mu*(2*n^2 + 4*n + 3))); 
+                %C3 = rho*(n-1)^2*(n+1)*n / (6*mu*a^n*(n-1)*(2*n^2+4*n+3) - 4*pi*G*rho^2*n*(3*a^(-n+2)-(2*n+1)*a^(n+2))); 
+                %C4 = -a^2*(n+2)*n/(n^2-1) * C3; 
+                %C5 = 1/a^n - 4*pi*G*rho*a^2/(n^2-1) * C3; 
+        end
 
-y1 = C1./r.^n + C2./r.^(n+2) + C3*r.^(n+1) + C4*r.^(n-1);
-y2 = 2*mu*( (-n^2-3*n+1)/(n+1)*C1*r.^(-n-1) - (n+2)*C2*r.^(-n-3) +...
-        (n^2-n-3)/n*C3*r.^n + (n-1)*C4*r.^(n-2) ) + 4/3*pi*G*rho^2*( C1*r.^(-n+1) +...
-        C2*r.^(-n-1) + C3*r.^(n+2) + C4*r.^n ) - rho*C5*r.^n - rho*C6*r.^(-n-1);
-y3 = -(n-2)/(n*(n+1))*C1*r.^(-n) - 1/(n+1)*C2*r.^(-n-2) + (n+3)/(n*(n+1))*C3*r.^(n+1) + C4/n*r.^(n-1);
-y4 = 2*mu*( (n-1)/n*C1*r.^(-n-1) + (n+2)/(n+1)*C2*r.^(-n-3) + (n+2)/(n+1)*C3*r.^n + (n-1)/n*C4*r.^(n-2) );
-y5 = C5*r.^n + C6*r.^(-n-1);
-y6 = n*C5*r.^(n-1) - (n+1)*C6*r.^(-n-2) - 4*pi*G*rho*( C1*r.^(-n) + C2*r.^(-n-2) + C3*r.^(n+1) + C4*r.^(n-1));
+        y1 = C1./r.^n + C2./r.^(n+2) + C3*r.^(n+1) + C4*r.^(n-1);
+        y2 = 2*mu*( (-n^2-3*n+1)/(n+1)*C1*r.^(-n-1) - (n+2)*C2*r.^(-n-3) +...
+                (n^2-n-3)/n*C3*r.^n + (n-1)*C4*r.^(n-2) ) + 4/3*pi*G*rho^2*( C1*r.^(-n+1) +...
+                C2*r.^(-n-1) + C3*r.^(n+2) + C4*r.^n ) - rho*C5*r.^n - rho*C6*r.^(-n-1);
+        y3 = -(n-2)/(n*(n+1))*C1*r.^(-n) - 1/(n+1)*C2*r.^(-n-2) + (n+3)/(n*(n+1))*C3*r.^(n+1) + C4/n*r.^(n-1);
+        y4 = 2*mu*( (n-1)/n*C1*r.^(-n-1) + (n+2)/(n+1)*C2*r.^(-n-3) + (n+2)/(n+1)*C3*r.^n + (n-1)/n*C4*r.^(n-2) );
+        y5 = C5*r.^n + C6*r.^(-n-1);
+        y6 = n*C5*r.^(n-1) - (n+1)*C6*r.^(-n-2) - 4*pi*G*rho*( C1*r.^(-n) + C2*r.^(-n-2) + C3*r.^(n+1) + C4*r.^(n-1));
 
-yana = [y1; y2; y3; y4; y5; y6]'; 
+        yana = [y1 y2 y3 y4 y5 y6];