diff --git a/README.md b/README.md index 6ffb4e3..0461a44 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Global Spherical Harmonic package (GSH) -GSH is a MATLAB package to do **Global Spherical Harmonic Analyses** (GSHA) and **Synthesis** (GSHS) for Crust1.0. +GSH is a MATLAB package to do **Global Spherical Harmonic Analyses** (GSHA) and **Synthesis** (GSHS) for Crust1.0. This version supports running the scripts with Octave. ## Requirements diff --git a/Tools/GSHA.m b/Tools/GSHA.m index 00a06fb..59b277e 100644 --- a/Tools/GSHA.m +++ b/Tools/GSHA.m @@ -1,92 +1,91 @@ -function cs = GSHA(f,lmax) - -% GSHA global spherical harmonic analysis (wls) -% -% HOW cs = gsha(f,method,grid) -% -% IN f - global field of size N*2N -% lmax - maximum degree of development -% OUT cs - Clm & Slm in |C\S| format -% -%-------------------------------------------------------------------------- -% uses Legendre_functions, SC2CS -%-------------------------------------------------------------------------- - -%-------------------------------------------------------------------------- -% Grid definition -%-------------------------------------------------------------------------- -[rows,cols] = size(f); - -n = rows; -dt = 180 / n; -theta = (dt/2:dt:180)'; -lam = (dt/2:dt:360); % dt = dlam - -%-------------------------------------------------------------------------- -% further diagnostics -%-------------------------------------------------------------------------- -if lmax > n - error('Maximum degree of development is higher than number of rows of input.') -end - -if lmax == n - error('max. degree 180 is not possible for block grid, problem for m = 0') -end -%-------------------------------------------------------------------------- -% Init. -%-------------------------------------------------------------------------- -% L = n; -L = lmax; -clm = zeros(L+1,L+1); -slm = zeros(L+1,L+1); - -%-------------------------------------------------------------------------- -% 1st step analysis: Am(theta) & Bm(theta) -%-------------------------------------------------------------------------- -m = 0:L; -c = cos(lam'*m*pi/180); -s = sin(lam'*m*pi/180); - -% preserving the orthogonality - -c = c/rows; -s = s/rows; - -c(:,1) = c(:,1)/2; -s(:,L+1) = s(:,L+1)/2; % not sure if this is needed? Need to take a look at the equations -%c(:,L+1) = zeros(2*n,1); % this line is not needed otherwise Clamx,lamx is not computed. -s(:,1) = zeros(2*n,1); - -a = f * c; -b = f * s; - -%-------------------------------------------------------------------------- -% 2nd step analysis: Clm & Slm -%-------------------------------------------------------------------------- - -% predefine weights for the weighted least squares analysis -si = sin(theta*pi/180); -si = 2*si/sum(si); - -% loop over the order of the Spherical Harmonics -for m = 0:L - % construct the legendre polynomials - p = Legendre_functions(m:L,m,theta); - % Select particular colum of the signal - ai = a(:,m+1); - bi = b(:,m+1); - d = 1:length(theta); - pts = p' * sparse(d,d,si); - - % Estimate the coefficients - clm(m+1:L+1,m+1) = (pts * p) \ pts * ai; - slm(m+1:L+1,m+1) = (pts * p) \ pts * bi; -end - -%-------------------------------------------------------------------------- -% Write the coefficients Clm & Slm in |C\S| format -%-------------------------------------------------------------------------- -slm = fliplr(slm); -%sc = [slm(:,1:L) clm]; -cs = sc2cs([slm(:,1:L) clm]); -cs = cs(1:lmax+1,1:lmax+1); \ No newline at end of file +function cs = GSHA(f,lmax) + +% GSHA global spherical harmonic analysis (wls) +% +% HOW cs = gsha(f,method,grid) +% +% IN f - global field of size N*2N +% lmax - maximum degree of development +% OUT cs - Clm & Slm in |C\S| format +% +%-------------------------------------------------------------------------- +% uses Legendre_functions, SC2CS +%-------------------------------------------------------------------------- + +%-------------------------------------------------------------------------- +% Grid definition +%-------------------------------------------------------------------------- +[rows,cols] = size(f); +n = rows; +dt = double(180 / n); +theta = double((dt/2:dt:180))'; +lam = double((dt/2:dt:360)); % dt = dlam +%-------------------------------------------------------------------------- +% further diagnostics +%-------------------------------------------------------------------------- +if lmax > n + error('Maximum degree of development is higher than number of rows of input.') +end + +if lmax == n + error('max. degree 180 is not possible for block grid, problem for m = 0') +end +%-------------------------------------------------------------------------- +% Init. +%-------------------------------------------------------------------------- +% L = n; +L = lmax; +clm = zeros(L+1,L+1); +slm = zeros(L+1,L+1); + +%-------------------------------------------------------------------------- +% 1st step analysis: Am(theta) & Bm(theta) +%-------------------------------------------------------------------------- +m = double(0:L); + +c = cos(lam'*m*pi/180); +s = sin(lam'*m*pi/180); + +% preserving the orthogonality + +c = c/rows; +s = s/rows; + +c(:,1) = c(:,1)/2; +s(:,L+1) = s(:,L+1)/2; % not sure if this is needed? Need to take a look at the equations +%c(:,L+1) = zeros(2*n,1); % this line is not needed otherwise Clamx,lamx is not computed. +s(:,1) = zeros(2*n,1); + +a = f * c; +b = f * s; + +%-------------------------------------------------------------------------- +% 2nd step analysis: Clm & Slm +%-------------------------------------------------------------------------- + +% predefine weights for the weighted least squares analysis +si = sin(theta*pi/180); +si = 2*si/sum(si); + +%% loop over the order of the Spherical Harmonics +for m = double(0:L) + % construct the legendre polynomials + p = Legendre_functions(m:L,m,theta); + % Select particular colum of the signal + ai = a(:,m+1); + bi = b(:,m+1); + d = 1:length(theta); + pts = p' * sparse(d,d,si); + + % Estimate the coefficients + clm(m+1:L+1,m+1) = (pts * p) \ (pts * ai); + slm(m+1:L+1,m+1) = (pts * p) \ (pts * bi); +end + +%-------------------------------------------------------------------------- +% Write the coefficients Clm & Slm in |C\S| format +%-------------------------------------------------------------------------- +slm = fliplr(slm); +%sc = [slm(:,1:L) clm]; +cs = sc2cs([slm(:,1:L) clm]); +cs = cs(1:lmax+1,1:lmax+1); diff --git a/Tools/GSHS.m b/Tools/GSHS.m index dc6e8ce..747df1a 100644 --- a/Tools/GSHS.m +++ b/Tools/GSHS.m @@ -53,9 +53,8 @@ end % prepare coordinates -th = sort(th(:)); -lam = lam(:).*pi/180; - +th = double(sort(th(:))); +lam = double(lam(:).*pi/180); % use the desired degree if nargin < 5 || isempty(ldesired), ldesired = lmax; end if ldesired < lmax @@ -68,7 +67,7 @@ [row,~] = size(field); % prepare cosine and sine --> cos(m*lam) and sin(m*lam) -m = 0:lmax; +m = double(0:lmax); l = m'; mlam = (lam*m)'; cosmlam = cos(mlam); @@ -81,7 +80,7 @@ %---------------------------------------------------------------------------- % CALCULATION %---------------------------------------------------------------------------- -for m = 0:row-1 +for m = double(0:row-1) if m==0 Cnm = field(:,row+m); % get column with order 0 Snm = zeros(row,1); % there are no Sn0 coefficients diff --git a/Tools/Legendre_functions.m b/Tools/Legendre_functions.m index e5aba22..8cab148 100644 --- a/Tools/Legendre_functions.m +++ b/Tools/Legendre_functions.m @@ -1,140 +1,141 @@ -function [p, dp, ddp] = Legendre_functions(l,m,th) - -% PLM Fully normalized associated Legendre functions for a selected order M -% -% HOW p = plm(l,th) - assumes M=0 -% p = plm(l,m,th) -% [p,dp] = plm(l,m,th) -% [p,ddp] = plm(l,m,th) -% -% -% IN l - degree (vector). Integer, but not necessarily monotonic. -% For l < m a vector of zeros will be returned. -% m - order (scalar). If absent, m=0 is assumed. -% th - co-latitude [deg] (vector) -% OUT p - Matrix with Legendre functions. The matrix has length(TH) rows -% and length(L) columns, unless L or TH is scalar. Then the output -% vector follows the shape of respectively L or TH. -% dp - Matrix with first derivative of Legendre functions. The matrix -% has length(TH) rows and length(L) columns, unless L or TH is -% scalar. Then the output vector follows the shape of respectively -% L or TH. -% ddp - Matrix with second (colatitude) derivative of Legendre functions. The matrix -% has length(TH) rows and length(L) columns, unless L or TH is -% scalar. Then the output vector follows the shape of respectively -% L or TH. -% -% Derivatives are based on GRafarend and Novak (2006) paper. -%----------------------------------------------------------------------------- -% Uses none -%----------------------------------------------------------------------------- -% Revision history: -% - code based on visu2plm_ww from Nico Sneeuw and Wouter van der Wal -% (23/11/2020) by Bart Root -%----------------------------------------------------------------------------- - - -% Input check for validatity of the program -if min(size(l)) ~= 1; error('Degree l must be vector (or scalar)'); end -if any(rem(l,1) ~= 0); error('Vector l contains non-integers.'); end -if max(size(m)) ~= 1; error('Order m must be scalar.'); end -if rem(m,1) ~=0; error('Order m must be integer.'); end - - -% Preliminaries. -[lrow,lcol] = size(l); -[trow,tcol] = size(th); -lmax = max(l); -if lmax < m; error('Largest degree still smaller than order m.'); end -x = cos(deg2rad(th(:))); -y = sin(deg2rad(th(:))); -lvec = l(:)'; % l can be used now as running index. - -% Recursive computation of the temporary matrix ptmp, containing the Legendre -% functions in its columns, with progressing degree l. The last column of -% ptmp will contain zeros, which is useful for assignments when l < m. -ptmp = zeros(length(th),lmax-m+2); -dptmp = zeros(length(th),lmax-m+2); -ddptmp = zeros(length(th),lmax-m+2); - -%-------------------------------------------------------------------- -% sectorial recursion: PM (non-recursive, though) -%-------------------------------------------------------------------- -% WW: produces sqrt( (2n+1)/2n ) -% Novak and Grafarend (2006), eq 64 -% recursion is beta_n,n * beta_n-1,n-1 * ... -% sin(theta) * sin(theta) * ... -% * cos(theta) * P_0,0 (which is 1, dP_0,0 is 0) -% note that the term beta_n,n*cos(theta)*P_n-1,n-1 of Novak and -% Grafarend(2006) equation 72 is taken care of by the m. - -% Calculate extra factor for derivatives of the Legendre functions. -if m == 0 - fac = 1; -else - mm = 2*(1:m); - fac = sqrt(2*prod((mm+1)./mm)); % extra sqrt(2) because summation not over negative orders -end - -% Calculate first column -ptmp(:,1) = fac*y.^m; % The 1st column of ptmp. -dptmp(:,1) = m*fac*(y.^(m-1).*x); -ddptmp(:,1) = -m*fac*(y.^m) + m*(m-1)*fac*(y.^(m-2).*x.^2); - - -%-------------------------------------------------------------------- -% l-recursion: P -%-------------------------------------------------------------------- -for l = m+1:lmax - col = l - m + 1; % points to the next column of ptmp - root1 = sqrt( (2*l+1)*(2*l-1)/((l-m)*(l+m)) ) ; % beta_n,m (65) - root2 = sqrt( (2*l+1)*(l+m-1)*(l-m-1) / ( (2*l-3)*(l-m)*(l+m) ) ); % beta_n,m (65) * gamma_n,m (66) - - % recursion P - if l == m+1 - ptmp(:,col) = root1 *x.*ptmp(:,col-1); - else - ptmp(:,col) = root1 *x.*ptmp(:,col-1) - root2 *ptmp(:,col-2); - end - - % recursion DP - if l == m+1 - dptmp(:,col) = root1 *(x.*dptmp(:,col-1)-y.*ptmp(:,col-1)); - else - dptmp(:,col) = root1 *(x.*dptmp(:,col-1)-y.*ptmp(:,col-1)) - root2 *dptmp(:,col-2); - end - - % recursion DDP - if l == m+1 - ddptmp(:,col) = root1 *(-x.*ptmp(:,col-1) -2*y.*dptmp(:,col-1) + x.*ddptmp(:,col-1) ); - else - ddptmp(:,col) = root1 *(-x.*ptmp(:,col-1) -2*y.*dptmp(:,col-1) + x.*ddptmp(:,col-1) ) - root2*ddptmp(:,col-2); - end - -end - - -% The Legendre functions have been computed. What remains to be done, is to -% extract the proper columns from ptmp, corresponding to the vector lvec. -% If l or theta is scalar the output matrix p reduces to a vector. It should -% have the shape of respectively theta or l in that case. - -lind = find(lvec < m); % index into l < m -pcol = lvec - m + 1; % index into columns of ptmp -pcol(lind) = (lmax-m+2)*ones(size(lind)); % Now l < m points to last col. - -p = ptmp(:,pcol); % proper column extraction -dp = dptmp(:,pcol); % proper column extraction -ddp = ddptmp(:,pcol); % proper column extraction - -% Change order of vectors -if max(size(lvec))==1 && min(size(th))==1 && (trow == 1) - p = p'; - dp = dp'; - ddp = ddp'; -end -if max(size(th))==1 && min(size(lvec))==1 && (lcol == 1) - p = p'; - dp = dp'; - ddp = ddp'; -end +function [p, dp, ddp] = Legendre_functions(l,m,th) + +% PLM Fully normalized associated Legendre functions for a selected order M +% +% HOW p = plm(l,th) - assumes M=0 +% p = plm(l,m,th) +% [p,dp] = plm(l,m,th) +% [p,ddp] = plm(l,m,th) +% +% +% IN l - degree (vector). Integer, but not necessarily monotonic. +% For l < m a vector of zeros will be returned. +% m - order (scalar). If absent, m=0 is assumed. +% th - co-latitude [deg] (vector) +% OUT p - Matrix with Legendre functions. The matrix has length(TH) rows +% and length(L) columns, unless L or TH is scalar. Then the output +% vector follows the shape of respectively L or TH. +% dp - Matrix with first derivative of Legendre functions. The matrix +% has length(TH) rows and length(L) columns, unless L or TH is +% scalar. Then the output vector follows the shape of respectively +% L or TH. +% ddp - Matrix with second (colatitude) derivative of Legendre functions. The matrix +% has length(TH) rows and length(L) columns, unless L or TH is +% scalar. Then the output vector follows the shape of respectively +% L or TH. +% +% Derivatives are based on GRafarend and Novak (2006) paper. +%----------------------------------------------------------------------------- +% Uses none +%----------------------------------------------------------------------------- +% Revision history: +% - code based on visu2plm_ww from Nico Sneeuw and Wouter van der Wal +% (23/11/2020) by Bart Root +%----------------------------------------------------------------------------- + +% Input check for validatity of the program +if min(size(l)) ~= 1; error('Degree l must be vector (or scalar)'); end +if any(rem(l,1) ~= 0); error('Vector l contains non-integers.'); end +if max(size(m)) ~= 1; error('Order m must be scalar.'); end +if rem(m,1) ~=0; error('Order m must be integer.'); end + + +% Preliminaries. +[lrow,lcol] = size(l); +[trow,tcol] = size(th); +lmax = max(l); +if lmax < m; error('Largest degree still smaller than order m.'); end +x = cos(deg2rad(th(:))); +y = sin(deg2rad(th(:))); +lvec = l(:)'; % l can be used now as running index. + +% Recursive computation of the temporary matrix ptmp, containing the Legendre +% functions in its columns, with progressing degree l. The last column of +% ptmp will contain zeros, which is useful for assignments when l < m. +ptmp = zeros(length(th),lmax-m+2); +dptmp = zeros(length(th),lmax-m+2); +ddptmp = zeros(length(th),lmax-m+2); + +%-------------------------------------------------------------------- +% sectorial recursion: PM (non-recursive, though) +%-------------------------------------------------------------------- +% WW: produces sqrt( (2n+1)/2n ) +% Novak and Grafarend (2006), eq 64 +% recursion is beta_n,n * beta_n-1,n-1 * ... +% sin(theta) * sin(theta) * ... +% * cos(theta) * P_0,0 (which is 1, dP_0,0 is 0) +% note that the term beta_n,n*cos(theta)*P_n-1,n-1 of Novak and +% Grafarend(2006) equation 72 is taken care of by the m. + +% Calculate extra factor for derivatives of the Legendre functions. +if m == 0 + fac = 1; +else + mm = 2*(1:m); + fac = sqrt(2*prod((mm+1)./mm)); % extra sqrt(2) because summation not over negative orders +end + +% Calculate first column +ptmp(:,1) = fac*y.^m; % The 1st column of ptmp. +dptmp(:,1) = m*fac*(y.^(m-1).*x); +ddptmp(:,1) = -m*fac*(y.^m) + m*(m-1)*fac*(y.^(m-2).*x.^2); + + +%-------------------------------------------------------------------- +% l-recursion: P +%-------------------------------------------------------------------- +for l = double(m+1:lmax) + col = l - m + 1; % points to the next column of ptmp + root1 = sqrt( (2*l+1)*(2*l-1)/((l-m)*(l+m)) ) ; % beta_n,m (65) + root2 = sqrt( (2*l+1)*(l+m-1)*(l-m-1) / ( (2*l-3)*(l-m)*(l+m) ) ); % beta_n,m (65) * gamma_n,m (66) + + % recursion P + if l == m+1 + ptmp(:,col) = root1 *x.*ptmp(:,col-1); + else + ptmp(:,col) = root1 *x.*ptmp(:,col-1) - root2 *ptmp(:,col-2); + end + + % recursion DP + if l == m+1 + dptmp(:,col) = root1 *(x.*dptmp(:,col-1)-y.*ptmp(:,col-1)); + else + dptmp(:,col) = root1 *(x.*dptmp(:,col-1)-y.*ptmp(:,col-1)) - root2 *dptmp(:,col-2); + end + + % recursion DDP + if l == m+1 + ddptmp(:,col) = root1 *(-x.*ptmp(:,col-1) -2*y.*dptmp(:,col-1) + x.*ddptmp(:,col-1) ); + else + ddptmp(:,col) = root1 *(-x.*ptmp(:,col-1) -2*y.*dptmp(:,col-1) + x.*ddptmp(:,col-1) ) - root2*ddptmp(:,col-2); + end + +end + + +% The Legendre functions have been computed. What remains to be done, is to +% extract the proper columns from ptmp, corresponding to the vector lvec. +% If l or theta is scalar the output matrix p reduces to a vector. It should +% have the shape of respectively theta or l in that case. + +lind = find(lvec < m); % index into l < m +pcol = lvec - m + 1; % index into columns of ptmp +pcol(lind) = (lmax-m+2)*ones(size(lind)); % Now l < m points to last col. + +p = ptmp(:,pcol); % proper column extraction +dp = dptmp(:,pcol); % proper column extraction +ddp = ddptmp(:,pcol); % proper column extraction + +% Change order of vectors +if max(size(lvec))==1 && min(size(th))==1 && (trow == 1) + p = p'; + dp = dp'; + ddp = ddp'; +end +if max(size(th))==1 && min(size(lvec))==1 && (lcol == 1) + p = p'; + dp = dp'; + ddp = ddp'; +end + +return diff --git a/Tools/gravityModule.m b/Tools/gravityModule.m index 1d2a4d5..55296ba 100644 --- a/Tools/gravityModule.m +++ b/Tools/gravityModule.m @@ -176,4 +176,4 @@ data.ten.Tzz = TenRR; data.ten.Txy = (TenTL + Trans1.*VecL./r); data.ten.Txz = -1.*(-VecT./r + TenRT); % might need a multiplication with -1 -data.ten.Tyz = -1.*(VecL./r - TenRL); % might need a multiplication with -1 \ No newline at end of file +data.ten.Tyz = -1.*(VecL./r - TenRL); % might need a multiplication with -1 diff --git a/Tools/import_layer.m b/Tools/import_layer.m index b33c98a..0daf45e 100644 --- a/Tools/import_layer.m +++ b/Tools/import_layer.m @@ -81,4 +81,4 @@ else error('No specific file is listed for the linear gradient in the layer') end -end \ No newline at end of file +end diff --git a/Tools/layer_SH_analysis.m b/Tools/layer_SH_analysis.m index 53f0a48..469a22e 100644 --- a/Tools/layer_SH_analysis.m +++ b/Tools/layer_SH_analysis.m @@ -1,352 +1,353 @@ -function [V] = layer_SH_analysis(nmax,geoid,Re,rhoE,max_bin,fupper,flower,fdens,falpha) -% -% This function converts a variable density layer into SH-coefficients. The -% following input variables are used: -% -% - nmax: maximum degree SH-coefficient (nmax=length(lat)-1) -% - geoid: string describing geoid used -% - 'none' (spherical geoid) -% - 'WGS84' -% - Re: value of radius [m] -% - rhoE: value of mean density [kg/m3] (linking to specific GM) -% - max_bin: maximum binomial terms that can be used (now max=3) -% -%%%%%%% Layer input values (format: matrix m*n [m=lat,n=lon (block grid definition)]): -% -% - fupper: value for upper boundary varying from Re [m] -% - flower: value for lower boundary varying from Re [m] -% - fdens: value for density [kg/m3] -% - falpha: value for linear density gradient [1/m] -% default=0 -% -%%%%%% output values: -% -% - V: Spherical harmonic coefficients (4pi-normalization) -% format: [l,m,Clm,Slm] l = 0,1,2,3,... m = 0,0,0,0,... etc. -% -% Made by: Bart Root, TUDelft, Astrodynamics and Space Missions -% First version 1.0: date: 08-04-2013 -% -% 06-feb-2014: Put in the anti-aliasing part, but don't use this with high -% degree of SH-coefficients. -% -% Matlab functions used: -% -% - GSHA -% - cs2sc -% - sc2vecml -% - geocradius -% - size -% - zeros -% -%%%%%%%%%%%%%%% Start function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% initializing certain values - -grid = 'block'; -method = 'ls'; -rhoRatio = 1./rhoE; - -% maximum binomial terms check (for now max = 3) - -if max_bin>8 - error('Maximum number of binomial terms is overriden') -end - -% No alpha stated, don't do extra the GSHA loops -if nargin<9 - radOn = 0; - falpha = 0; -else - radOn = 1; -end - -% Do all the matrices have the same format -if size(fupper)~=size(flower) - error('Boundary matrices are not the same size') -elseif size(fdens)~=size(fupper) - error('Density matrix is not the same size as the boundary matrices') -elseif ( radOn==1 ) - if size(fdens)~=size(falpha) - error('Linear density gradient matrix is not the same size as the boundary matrices') - end -else - % all is correct: do nothing -end - -% Check if nmax value is correct -if nmax+1>size(fupper,1) - error('Size input matrix is incorrect: Latitude length is to small for nmax') -elseif nmax+1*2>size(fupper,2) - error('Size input matrix is incorrect: Longitude length is to small for nmax') -end - -% Geoidetic radius reference correction -if strcmp(geoid,'none') - % do nothing: spherical case -elseif strcmp(geoid,'WGS84') - res = 180/size(fupper,1); - latV = ((-90 +res/2:res:90-res/2).*-1)'; - lat = repmat(latV,1,size(fupper,2)); - - correctionR = geocradius(lat,'WGS84') - Re; - fupper = fupper + correctionR; - correctionR = geocradius(lat,'WGS84') - Re; - flower = flower + correctionR; - clear lat - -elseif strcmp(geoid,'benchmark_Mikhail') - % not needed anymore! - disp('Benchmark Mikhail reference') - - f = 1/298.257223563; - aRe = 6371000; - - res = 180/size(fupper,1); - latV = ((-90 +res/2:res:90-res/2).*-1)'; - lat = repmat(latV,1,size(fupper,2)); - - correctionR = geocradius(lat,f,aRe) - Re; - fupper = fupper + correctionR; - correctionR = geocradius(lat,f,aRe) - Re; - flower = flower + correctionR; - clear lat - -else - % An incorrect geoid is chosen - error('An incorrect geoid reference is chosen') -end - -%%%% %%%% -%%%% Starting with the binomial terms loop in the Analysis %%%% -%%%% %%%% - -for Hi = 1:max_bin - - %disp(['binary number: ' num2str(Hi)]) - - % calibrate density if alpha is used -% if radOn ==1 -% fdens = fdens./(1+falpha .* ((fupper+flower)./2)); -% end - - % K = Hi - switch Hi - case 1 - m = (fdens).*(fupper-flower); - case 2 - m = (fdens).*(fupper.*fupper-flower.*flower); - case 3 - m = (fdens).*(fupper.*fupper.*fupper-flower.*flower.*flower); - case 4 - m = (fdens).*(fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower); - case 5 - m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower); - case 6 - m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower.*flower); - case 7 - m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower.*flower.*flower); - case 8 - m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower.*flower.*flower.*flower); - otherwise - error('Incorrect number of binomial terms') - end - -% %%%%%%% test anti-aliasing -% m_a = zeros(size(m)*4); -% -% for ii = 1:size(m,1) -% for jj = 1:size(m,2) -% -% i = 1 + (ii-1) * 4; -% j = 1 + (jj-1) * 4; -% -% D = m(ii,jj)*ones(4,4); -% -% m_a(i:i+(4-1),j:j+(4-1)) = D; -% end -% end - - - % Do the analysis - cs = GSHA(m,nmax); sc = cs2sc(cs); - %figure;plot(1:size(sc,1),log10(sqrt(sum(sc.^2,2))./(1:size(sc,1))')); - [Clm,Slm,llvec,mmvec] = sc2vecml(sc,nmax); - - %%%%%%%%%%%%%%%%%% - Vpart = [llvec' mmvec' Clm Slm]; - - % Stating the factors - fac = (3./(((2.*Vpart(:,1))+1))); - fac2 = (((Vpart(:,1)+2)/2)); - fac3 = (((Vpart(:,1)+2).*(Vpart(:,1)+1)/6)); - fac4 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1))/24)); - fac5 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1)/120)); - fac6 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1).*(Vpart(:,1)-2)/720)); - fac7 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1).*(Vpart(:,1)-2).*(Vpart(:,1)-3)/5040)); - fac8 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1).*(Vpart(:,1)-2).*(Vpart(:,1)-3).*(Vpart(:,1)-4)/40320)); - - % If needed also analyse the linear density gradient term - if radOn ==1 - if Hi==2 - % Do the SHA on the radial dens component. - % K = 2 - m = (falpha.*Re).*(fdens).*(fupper.*fupper-flower.*flower); -% %%%%%%% test anti-aliasing -% m_a = zeros(size(m)*4); -% -% for ii = 1:size(m,1) -% for jj = 1:size(m,2) -% -% i = 1 + (ii-1) * 4; -% j = 1 + (jj-1) * 4; -% -% D = m(ii,jj)*ones(4,4); -% -% m_a(i:i+(4-1),j:j+(4-1)) = D; -% end -% end - - cs = GSHA(m,nmax); sc = cs2sc(cs); - [Clm,Slm,llvec,mmvec] = sc2vecml(sc,nmax); - - %%%%%%%%%%%%%%%%%% - Vpart_rad = [llvec' mmvec' Clm Slm]; - - % Adding the radial component to the other coefficients - %Vpart(:,3:4) = Vpart(:,3:4) + Vpart_rad(:,3:4); - elseif Hi==3 - % K = 3 - m = (2.*falpha.*Re).*(fdens).*(fupper.*fupper.*fupper-flower.*flower.*flower); -% %%%%%%% test anti-aliasing -% m_a = zeros(size(m)*4); -% -% for ii = 1:size(m,1) -% for jj = 1:size(m,2) -% -% i = 1 + (ii-1) * 4; -% j = 1 + (jj-1) * 4; -% -% D = m(ii,jj)*ones(4,4); -% -% m_a(i:i+(4-1),j:j+(4-1)) = D; -% end -% end - - cs = GSHA(m,nmax); sc = cs2sc(cs); - [Clm,Slm,llvec,mmvec] = sc2vecml(sc,nmax); - - %%%%%%%%%%%%%%%%%% - Vpart_rad = [llvec' mmvec' (Vpart(:,1)+2).*Clm (Vpart(:,1)+2).*Slm]; - - % Adding the radial component to the other coefficients - %Vpart(:,3:4) = Vpart(:,3:4) + Vpart_rad(:,3:4); - end - - end - - % Constructing the correct coefficients - if Hi==1 - V = zeros(size(Vpart)); - V(:,1) = Vpart(:,1); - V(:,2) = Vpart(:,2); - end - - switch Hi - case 1 - VC = (Vpart(:,3)).*fac.*rhoRatio./Re; - VS = (Vpart(:,4)).*fac.*rhoRatio./Re; - case 2 - VC = fac2.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re; - VS = fac2.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re; - - % with lineair density variation - if radOn ==1 - VC = VC + (Vpart_rad(:,3)).*fac.*rhoRatio./Re./Re/2; - VS = VS + (Vpart_rad(:,4)).*fac.*rhoRatio./Re./Re/2; - end - - case 3 - VC = fac3.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re; - VS = fac3.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re; - - % with lineair density variation - if radOn ==1 - VC = VC + (Vpart_rad(:,3)).*fac.*rhoRatio./Re./Re./Re/6; - VS = VS + (Vpart_rad(:,4)).*fac.*rhoRatio./Re./Re./Re/6; - end - - case 4 - VC = fac4.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re; - VS = fac4.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re; - case 5 - VC = fac5.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re; - VS = fac5.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re; - case 6 - VC = fac6.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re; - VS = fac6.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re; - case 7 - VC = fac7.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re; - VS = fac7.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re; - case 8 - VC = fac8.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re./Re; - VS = fac8.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re./Re; - otherwise - error('Incorrect number of binomial terms') - end - - % Summate the coefficients - switch Hi - case 1 - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - case 2 - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - case 3 - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - case 4 - VC(Vpart(:,1)==0) = 0; - VS(Vpart(:,1)==0) = 0; - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - case 5 - VC(Vpart(:,1)==0|Vpart(:,1)==1) = 0; - VS(Vpart(:,1)==0|Vpart(:,1)==1) = 0; - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - case 6 - VC(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2) = 0; - VS(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2) = 0; - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - case 7 - VC(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3) = 0; - VS(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3) = 0; - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - case 8 - VC(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3|Vpart(:,1)==4) = 0; - VS(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3|Vpart(:,1)==4) = 0; - V(:,3) = V(:,3) + VC; - V(:,4) = V(:,4) + VS; - otherwise - error('Incorrect number of binomial terms') - end - -% figure -% plot(log10(abs(VC(1:180)))) -% hold on -% %plot(log10(abs(VS(1:180,4))),'r') -% hold off -end - -% figure -% plot(log10(fac(1:180).*rhoRatio./Re)) -% hold on -% plot(log10(fac2(1:180).*fac(1:180).*rhoRatio./Re./Re),'r') -% plot(log10(fac3(1:180).*fac(1:180).*rhoRatio./Re./Re./Re),'g') -% plot(log10(fac4(1:180).*fac(1:180).*rhoRatio./Re./Re./Re./Re),'k') -% plot(log10(fac5(1:180).*fac(1:180).*rhoRatio./Re./Re./Re./Re./Re),'m') -% hold off +function [V] = layer_SH_analysis(nmax,geoid,Re,rhoE,max_bin,fupper,flower,fdens,falpha) +% +% This function converts a variable density layer into SH-coefficients. The +% following input variables are used: +% +% - nmax: maximum degree SH-coefficient (nmax=length(lat)-1) +% - geoid: string describing geoid used +% - 'none' (spherical geoid) +% - 'WGS84' +% - Re: value of radius [m] +% - rhoE: value of mean density [kg/m3] (linking to specific GM) +% - max_bin: maximum binomial terms that can be used (now max=3) +% +%%%%%%% Layer input values (format: matrix m*n [m=lat,n=lon (block grid definition)]): +% +% - fupper: value for upper boundary varying from Re [m] +% - flower: value for lower boundary varying from Re [m] +% - fdens: value for density [kg/m3] +% - falpha: value for linear density gradient [1/m] +% default=0 +% +%%%%%% output values: +% +% - V: Spherical harmonic coefficients (4pi-normalization) +% format: [l,m,Clm,Slm] l = 0,1,2,3,... m = 0,0,0,0,... etc. +% +% Made by: Bart Root, TUDelft, Astrodynamics and Space Missions +% First version 1.0: date: 08-04-2013 +% +% 06-feb-2014: Put in the anti-aliasing part, but don't use this with high +% degree of SH-coefficients. +% +% Matlab functions used: +% +% - GSHA +% - cs2sc +% - sc2vecml +% - geocradius +% - size +% - zeros +% +%%%%%%%%%%%%%%% Start function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% initializing certain values + +grid = 'block'; +method = 'ls'; +rhoRatio = 1./rhoE; + +% maximum binomial terms check (for now max = 3) + +if max_bin>8 + error('Maximum number of binomial terms is overriden') +end + +% No alpha stated, don't do extra the GSHA loops +if nargin<9 + radOn = 0; + falpha = 0; +else + radOn = 1; +end + +% Do all the matrices have the same format +if size(fupper)~=size(flower) + error('Boundary matrices are not the same size') +elseif size(fdens)~=size(fupper) + error('Density matrix is not the same size as the boundary matrices') +elseif ( radOn==1 ) + if size(fdens)~=size(falpha) + error('Linear density gradient matrix is not the same size as the boundary matrices') + end +else + % all is correct: do nothing +end + +% Check if nmax value is correct +if nmax+1>size(fupper,1) + error('Size input matrix is incorrect: Latitude length is to small for nmax') +elseif nmax+1*2>size(fupper,2) + error('Size input matrix is incorrect: Longitude length is to small for nmax') +end + +% Geoidetic radius reference correction +if strcmp(geoid,'none') + % do nothing: spherical case +elseif strcmp(geoid,'WGS84') + res = 180/size(fupper,1); + latV = ((-90 +res/2:res:90-res/2).*-1)'; + lat = repmat(latV,1,size(fupper,2)); + + correctionR = geocradius(lat,'WGS84') - Re; + fupper = fupper + correctionR; + correctionR = geocradius(lat,'WGS84') - Re; + flower = flower + correctionR; + clear lat + +elseif strcmp(geoid,'benchmark_Mikhail') + % not needed anymore! + disp('Benchmark Mikhail reference') + + f = 1/298.257223563; + aRe = double(6371000); + + res = 180/size(fupper,1); + latV = ((-90 +res/2:res:90-res/2).*-1)'; + lat = repmat(latV,1,size(fupper,2)); + + correctionR = geocradius(lat,f,aRe) - Re; + fupper = fupper + correctionR; + correctionR = geocradius(lat,f,aRe) - Re; + flower = flower + correctionR; + clear lat + +else + % An incorrect geoid is chosen + error('An incorrect geoid reference is chosen') +end + +%%%% %%%% +%%%% Starting with the binomial terms loop in the Analysis %%%% +%%%% %%%% + +for Hi = 1:max_bin + + %disp(['binary number: ' num2str(Hi)]) + + % calibrate density if alpha is used +% if radOn ==1 +% fdens = fdens./(1+falpha .* ((fupper+flower)./2)); +% end + + % K = Hi + switch Hi + case 1 + m = (fdens).*(fupper-flower); + case 2 + m = (fdens).*(fupper.*fupper-flower.*flower); + case 3 + m = (fdens).*(fupper.*fupper.*fupper-flower.*flower.*flower); + case 4 + m = (fdens).*(fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower); + case 5 + m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower); + case 6 + m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower.*flower); + case 7 + m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower.*flower.*flower); + case 8 + m = (fdens).*(fupper.*fupper.*fupper.*fupper.*fupper.*fupper.*fupper.*fupper-flower.*flower.*flower.*flower.*flower.*flower.*flower.*flower); + otherwise + error('Incorrect number of binomial terms') + end + +% %%%%%%% test anti-aliasing +% m_a = zeros(size(m)*4); +% +% for ii = 1:size(m,1) +% for jj = 1:size(m,2) +% +% i = 1 + (ii-1) * 4; +% j = 1 + (jj-1) * 4; +% +% D = m(ii,jj)*ones(4,4); +% +% m_a(i:i+(4-1),j:j+(4-1)) = D; +% end +% end + + + % Do the analysis + cs = GSHA(double(m),nmax); sc = cs2sc(cs); + %figure;plot(1:size(sc,1),log10(sqrt(sum(sc.^2,2))./(1:size(sc,1))')); + + [Clm,Slm,llvec,mmvec] = sc2vecml(sc,nmax); + + %%%%%%%%%%%%%%%%%% + Vpart = [llvec' mmvec' Clm Slm]; + + % Stating the factors + fac = (3./(((2.*Vpart(:,1))+1))); + fac2 = (((Vpart(:,1)+2)/2)); + fac3 = (((Vpart(:,1)+2).*(Vpart(:,1)+1)/6)); + fac4 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1))/24)); + fac5 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1)/120)); + fac6 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1).*(Vpart(:,1)-2)/720)); + fac7 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1).*(Vpart(:,1)-2).*(Vpart(:,1)-3)/5040)); + fac8 = (((Vpart(:,1)+2).*(Vpart(:,1)+1).*(Vpart(:,1)).*(Vpart(:,1)-1).*(Vpart(:,1)-2).*(Vpart(:,1)-3).*(Vpart(:,1)-4)/40320)); + + % If needed also analyse the linear density gradient term + if radOn ==1 + if Hi==2 + % Do the SHA on the radial dens component. + % K = 2 + m = (falpha.*Re).*(fdens).*(fupper.*fupper-flower.*flower); +% %%%%%%% test anti-aliasing +% m_a = zeros(size(m)*4); +% +% for ii = 1:size(m,1) +% for jj = 1:size(m,2) +% +% i = 1 + (ii-1) * 4; +% j = 1 + (jj-1) * 4; +% +% D = m(ii,jj)*ones(4,4); +% +% m_a(i:i+(4-1),j:j+(4-1)) = D; +% end +% end + + cs = GSHA(m,nmax); sc = cs2sc(cs); + [Clm,Slm,llvec,mmvec] = sc2vecml(sc,nmax); + + %%%%%%%%%%%%%%%%%% + Vpart_rad = [llvec' mmvec' Clm Slm]; + + % Adding the radial component to the other coefficients + %Vpart(:,3:4) = Vpart(:,3:4) + Vpart_rad(:,3:4); + elseif Hi==3 + % K = 3 + m = (2.*falpha.*Re).*(fdens).*(fupper.*fupper.*fupper-flower.*flower.*flower); +% %%%%%%% test anti-aliasing +% m_a = zeros(size(m)*4); +% +% for ii = 1:size(m,1) +% for jj = 1:size(m,2) +% +% i = 1 + (ii-1) * 4; +% j = 1 + (jj-1) * 4; +% +% D = m(ii,jj)*ones(4,4); +% +% m_a(i:i+(4-1),j:j+(4-1)) = D; +% end +% end + + cs = GSHA(m,nmax); sc = cs2sc(cs); + [Clm,Slm,llvec,mmvec] = sc2vecml(sc,nmax); + + %%%%%%%%%%%%%%%%%% + Vpart_rad = [llvec' mmvec' (Vpart(:,1)+2).*Clm (Vpart(:,1)+2).*Slm]; + + % Adding the radial component to the other coefficients + %Vpart(:,3:4) = Vpart(:,3:4) + Vpart_rad(:,3:4); + end + + end + + % Constructing the correct coefficients + if Hi==1 + V = zeros(size(Vpart)); + V(:,1) = Vpart(:,1); + V(:,2) = Vpart(:,2); + end + + switch Hi + case 1 + VC = (Vpart(:,3)).*fac.*rhoRatio./Re; + VS = (Vpart(:,4)).*fac.*rhoRatio./Re; + case 2 + VC = fac2.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re; + VS = fac2.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re; + + % with lineair density variation + if radOn ==1 + VC = VC + (Vpart_rad(:,3)).*fac.*rhoRatio./Re./Re/2; + VS = VS + (Vpart_rad(:,4)).*fac.*rhoRatio./Re./Re/2; + end + + case 3 + VC = fac3.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re; + VS = fac3.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re; + + % with lineair density variation + if radOn ==1 + VC = VC + (Vpart_rad(:,3)).*fac.*rhoRatio./Re./Re./Re/6; + VS = VS + (Vpart_rad(:,4)).*fac.*rhoRatio./Re./Re./Re/6; + end + + case 4 + VC = fac4.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re; + VS = fac4.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re; + case 5 + VC = fac5.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re; + VS = fac5.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re; + case 6 + VC = fac6.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re; + VS = fac6.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re; + case 7 + VC = fac7.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re; + VS = fac7.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re; + case 8 + VC = fac8.*(Vpart(:,3)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re./Re; + VS = fac8.*(Vpart(:,4)).*fac.*rhoRatio./Re./Re./Re./Re./Re./Re./Re./Re; + otherwise + error('Incorrect number of binomial terms') + end + + % Summate the coefficients + switch Hi + case 1 + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + case 2 + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + case 3 + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + case 4 + VC(Vpart(:,1)==0) = 0; + VS(Vpart(:,1)==0) = 0; + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + case 5 + VC(Vpart(:,1)==0|Vpart(:,1)==1) = 0; + VS(Vpart(:,1)==0|Vpart(:,1)==1) = 0; + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + case 6 + VC(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2) = 0; + VS(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2) = 0; + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + case 7 + VC(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3) = 0; + VS(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3) = 0; + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + case 8 + VC(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3|Vpart(:,1)==4) = 0; + VS(Vpart(:,1)==0|Vpart(:,1)==1|Vpart(:,1)==2|Vpart(:,1)==3|Vpart(:,1)==4) = 0; + V(:,3) = V(:,3) + VC; + V(:,4) = V(:,4) + VS; + otherwise + error('Incorrect number of binomial terms') + end + +% figure +% plot(log10(abs(VC(1:180)))) +% hold on +% %plot(log10(abs(VS(1:180,4))),'r') +% hold off +end + +% figure +% plot(log10(fac(1:180).*rhoRatio./Re)) +% hold on +% plot(log10(fac2(1:180).*fac(1:180).*rhoRatio./Re./Re),'r') +% plot(log10(fac3(1:180).*fac(1:180).*rhoRatio./Re./Re./Re),'g') +% plot(log10(fac4(1:180).*fac(1:180).*rhoRatio./Re./Re./Re./Re),'k') +% plot(log10(fac5(1:180).*fac(1:180).*rhoRatio./Re./Re./Re./Re./Re),'m') +% hold off diff --git a/Tools/model_SH_analysis.m b/Tools/model_SH_analysis.m index 4216130..b233b3a 100644 --- a/Tools/model_SH_analysis.m +++ b/Tools/model_SH_analysis.m @@ -33,20 +33,22 @@ %G = 6.67384e-11; % Universal gravity constant %G = 6.6732e-11; %G = 6.67428e-11; -G = 6.673e-11; + +Model.Re = double(Model.Re); +Model.GM = double(Model.GM); +G = double(6.673e-11); Re3 = (Model.Re).^3; -rhoE = 3.*Model.GM./(4*pi*G*Re3); +rhoE = 3*Model.GM/(4*pi*G*Re3); % Default three binomial series terms (This can be modified after modifying layer_SH_analysis.m) -max_bin = 3; +max_bin = double(3); % starting the loop over the different layers for nlayer = 1:Model.number_of_layers layer_name = ['l' num2str(nlayer)]; - disp(['The layer number ' layer_name ' is starting']) % Costruct the coefficients for that particular layer [U,L,R] = import_layer(Model,nlayer); @@ -86,5 +88,4 @@ V(:,3:4) = V(:,3:4) + Vlayer(:,3:4); end - diff --git a/Tools/model_SH_synthesis.m b/Tools/model_SH_synthesis.m index a8b21f4..aa6353c 100644 --- a/Tools/model_SH_synthesis.m +++ b/Tools/model_SH_synthesis.m @@ -22,7 +22,7 @@ Lon = repmat(lon,length(lat),1); Lat = repmat(lat',1,length(lon)); -r = Model.Re + height; +r = Model.Re + double(height); % make the Stokes coefficients independend of sorting diff --git a/Tools/sc2vecml.m b/Tools/sc2vecml.m index 7048f2d..27ff5ce 100644 --- a/Tools/sc2vecml.m +++ b/Tools/sc2vecml.m @@ -23,8 +23,8 @@ Clm = zeros( (lmax+2)*(lmax+1)/2 ,1 ); Slm = zeros( (lmax+2)*(lmax+1)/2 ,1 ); -for mm = 0:lmax - for ll = mm:lmax +for mm = double(0:lmax) + for ll = double(mm:lmax) l1 = l1+1; if l1<= length(Clm) Clm(l1) = sc(ll+1,lmax+1+mm); diff --git a/Tools/topo2crust.m b/Tools/topo2crust.m index 4ece92f..4c2e011 100644 --- a/Tools/topo2crust.m +++ b/Tools/topo2crust.m @@ -1,4 +1,4 @@ -function [CM_interface,lon_CM,lat_CM] = topo2crust(h,N,CM_type,Model) +function [CM_interface, lon_CM, lat_CM] = topo2crust(h,N,CM_type,Model) % this function initiates the global spherical harmonics analysis and synthesis % for general purpose | inserted functions are written by Bart Root, and @@ -6,17 +6,17 @@ %% define parameters and constants for Mars -Te = Model.Te; % Effective Elastic thickness of lithosphere[m] -T = Model.D_c; % Crustal thickness with zero topography [m] -E = Model.E; % Young's modulus -v = Model.v; % Poisson's ratio -rho_m = Model.rho_m; % mantle density [kg/m^3] -rho_c = Model.rho_c; % crustal density [kg/m^3] -g = Model.GM./Model.Re^2; % gravity m/s^2 -R = Model.Re; % radius in [m] +Te = double(Model.Te); % Effective Elastic thickness of lithosphere[m] +T = double(Model.D_c); % Crustal thickness with zero topography [m] +E = double(Model.E); % Young's modulus +v = double(Model.v); % Poisson's ratio +rho_m = double(Model.rho_m); % mantle density [kg/m^3] +rho_c = double(Model.rho_c); % crustal density [kg/m^3] +g = double(Model.GM./Model.Re^2); % gravity m/s^2 +R = double(Model.Re); % radius in [m] %% Calculate some variables -d_rho = rho_m - rho_c; % density difference [kg/m^3]f +d_rho = double(rho_m - rho_c); % density difference [kg/m^3]f D = E*Te^3/(12*(1-v^2)); % Flexural rigidity disp(['The Flexural Rigidity used is: ' num2str(D./1e9,3) ' GPa m^3']) @@ -27,7 +27,7 @@ %% Run GSHA (Analysis) % Clm & Slm in |C\S| format [matrix] -cs = GSHA(root,N); +cs = GSHA(root,N); % converts spherical harmonics coefficients in |C\S| storage format into % a rectangular (L+1)x(2L+1) matrix in /S|C\format. @@ -69,8 +69,8 @@ % is the disturbing potential and any derivative up to the fourth. res = 180/(N+1); -lam = [res/2:res:360-res/2]; % lam [n x 1] longitude [deg] -th = [res/2:res:180-res/2]; % th [m x 1] co-latitude [deg] 【90-lat】 +lam = double([res/2:res:360-res/2]); % lam [n x 1] longitude [deg] +th = double([res/2:res:180-res/2]); % th [m x 1] co-latitude [deg] 【90-lat】 % call the function CM_root = GSHS(sc_CM,lam,th,N); @@ -79,3 +79,4 @@ lon_CM = lam; lat_CM = 90-th; + diff --git a/Tools/vecml2sc.m b/Tools/vecml2sc.m index 6cbff2d..2a8f25c 100644 --- a/Tools/vecml2sc.m +++ b/Tools/vecml2sc.m @@ -18,8 +18,8 @@ l1 = 0; sc = zeros(lmax+1,2*lmax+1); -for mm = 0:lmax - for ll = mm:lmax +for mm = double(0:lmax) + for ll = double(mm:lmax) l1 = l1+1; if l1<= length(Clm) sc(ll+1,lmax+1+mm) = Clm(l1);