Index: ../trunk-jpl/test/NightlyRun/test2085.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test2085.m (revision 25302) +++ ../trunk-jpl/test/NightlyRun/test2085.m (revision 25303) @@ -10,10 +10,10 @@ md.materials.numlayers = 10; md.love.forcing_type = 9; - md.materials.density= (1:md.materials.numlayers)'*0+5511; - md.materials.lame_mu= (1:md.materials.numlayers)'*0+0.75e11; - md.materials.viscosity=(1:md.materials.numlayers)'*0+1e21; - md.materials.lame_lambda=md.materials.lame_mu*0+5e17; + md.materials.density=zeros(md.materials.numlayers,1)+5511; + md.materials.lame_mu=zeros(md.materials.numlayers,1)+0.75e11; + md.materials.viscosity=zeros(md.materials.numlayers,1)+1e21; + md.materials.lame_lambda=zeros(md.materials.numlayers,1)+5e17; md.materials.issolid=ones(md.materials.numlayers,1); md.materials.isburgers=zeros(md.materials.numlayers,1); md.materials.burgers_mu=md.materials.lame_mu/3; @@ -24,7 +24,7 @@ md.love.allow_layer_deletion=1; md.love.frequencies=0; - md.love.nfreq=length(md.love.frequencies); + md.love.nfreq=length(md.love.frequencies); % TODO: Why are we grabbing length of a scalar here? md.love.sh_nmin = 2; md.love.sh_nmax = 2; @@ -37,7 +37,7 @@ md=solve(md,'lv'); % extract love kernels {{{ - y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1))); + y1=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1))); y2=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2))); y3=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3))); y4=squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,4))); @@ -147,18 +147,18 @@ 'y6_load_degree2_interface1','y6_load_degree2_interface2','y6_load_degree2_interface7','y6_load_degree2_interface10','y6_load_degree2_interface11',... }; field_tolerances={... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... - 1e-10,1e-10,1e-10,1e-10,1e-10,... + 9e-8, 1e-7, 3e-7, 3e-7, 3e-7,... + 9e-8, 1e-7, 3e-7, 3e-7, 1e-10,... + 9e-8, 8e-8, 2e-8, 2e-7, 4e-7,... + 9e-8, 9e-8, 2e-7, 4e-7, 1e-10,... + 4e-7, 4e-7, 2e-7, 3e-8, 2e-8,... + 2e-5, 2e-6, 2e-6, 1e-6, 2e-7,... + 3e-6, 3e-6, 3e-6, 4e-6, 4e-6,... + 3e-6, 3e-6, 2e-6, 6e-7, 1e-10,... + 3e-6, 3e-6, 5e-7, 3e-6, 5e-6,... + 3e-6, 3e-6, 9e-7, 7e-7, 1e-10,... + 4e-6, 4e-6, 3e-6, 5e-7, 3e-7,... + 3e-6, 3e-6, 2e-6, 7e-7, 2e-7... }; field_values={... y1_tidal_degree2_interface1, y1_tidal_degree2_interface2, y1_tidal_degree2_interface7, y1_tidal_degree2_interface10,y1_tidal_degree2_interface11,... Index: ../trunk-jpl/test/NightlyRun/test2085.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test2085.py (revision 25302) +++ ../trunk-jpl/test/NightlyRun/test2085.py (revision 25303) @@ -1,49 +1,188 @@ -#Test Name: LovenumberstAtDepth. -#Same as test #1 of test2084.m +# Test Name: FourierLoveKernels +# Homogenous Earth, for which analytic solutions exist. +# Love kernels for degree 2 (tested against analytic solns). -from model import * + from socket import gethostname -from solve import * -from numpy import * + +import numpy as np + from generic import generic from materials import * +from model import * +from solve import * +# Fro volumetric potential md = model() -md.cluster = generic('name', gethostname(), 'np', 1) md.materials = materials('litho') -md.miscellaneous.name = 'FourierLoveTest' -md.groundingline.migration = 'None' -md.verbose = verbose('111111101') -cst = 365.25 * 24 * 3600 * 1000 +md.materials.numlayers = 10 +md.love.forcing_type = 9 -md.materials.numlayers = 6 -md.materials.radius = np.array([10, 1222.5, 3.4800e+03, 5.7010e+03, 5.9510e+03, 6.3010e+03, 6.3710e+03]).reshape(-1, 1) * 1e3 -md.materials.density = np.array([1.0750e4, 1.0750e+04, 4.9780e+03, 3.8710e+03, 3.4380e+03, 3.0370e+03]).reshape(-1, 1) -md.materials.lame_mu = np.array([1e-5, 0, 2.2834e+00, 1.0549e+00, 7.0363e-01, 5.0605e-01]).reshape(-1, 1) * 1e11 -md.materials.viscosity = np.array([0, 0, 2.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+25]).reshape(-1, 1) * 1e21 -md.materials.lame_lambda = np.array(md.materials.lame_mu) * 0 + 5e14 -md.materials.issolid = np.array([1, 0, 1, 1, 1, 1]).reshape(-1, 1) -md.materials.isburgers = np.zeros((md.materials.numlayers)).reshape(-1, 1) +md.materials.density = np.zeros((md.materials.numlayers, 1)) + 5511 +md.materials.lame_mu = np.zeros((md.materials.numlayers, 1)) + 0.75e11 +md.materials.viscosity = np.zeros((md.materials.numlayers, 1)) + 1e21 +md.materials.lame_lambda = np.zeros((md.materials.numlayers, 1)) + 0.5e17 +md.materials.issolid = np.ones((md.materials.numlayers, 1)) +md.materials.isburgers = np.zeros((md.materials.numlayers, 1)) +md.materials.burgers_mu = md.materials.lame_mu / 3 +md.materials.burgers_viscosity = md.materials.viscosity / 10 -md.love.love_kernels = 1 +md.materials.radius = np.linspace(10e3, 6371e3, md.materials.numlayers + 1).reshape(-1, 1) +md.love.g0 = 9.8134357285509388 # directly grabbed from fourierlovesolver for this particular case + md.love.allow_layer_deletion = 1 -md.love.frequencies = (np.array([0]) * 2 * pi).reshape(-1, 1) / cst -md.love.nfreq = len(md.love.frequencies) +md.love.frequencies = 0 +md.love.nfreq = 1 + +md.love.sh_nmin = 2 md.love.sh_nmax = 2 +md.love.love_kernels = 1 -md.materials.burgers_mu = md.materials.lame_mu -md.materials.burgers_viscosity = md.materials.viscosity +md.miscellaneous.name = 'kernels' +md.cluster = generic('name', gethostname(), 'np', 1) +md.verbose = verbose('111111101') -md = solve(md, 'lv') +md = solve(md,'lv') +# Extract love kernels #{{{ +y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze() +y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze() +y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze() +y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze() +y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze() +y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze() + +y1_tidal_degree2_interface1 = y1[2, 0] +y1_tidal_degree2_interface2 = y1[2, 1] +y1_tidal_degree2_interface7 = y1[2, 6] +y1_tidal_degree2_interface10 = y1[2, 9] +y1_tidal_degree2_interface11 = y1[2, 10] + +y2_tidal_degree2_interface1 = y2[2, 0] +y2_tidal_degree2_interface2 = y2[2, 1] +y2_tidal_degree2_interface7 = y2[2, 6] +y2_tidal_degree2_interface10 = y2[2, 9] +y2_tidal_degree2_interface11 = y2[2, 10] + +y3_tidal_degree2_interface1 = y3[2, 0] +y3_tidal_degree2_interface2 = y3[2, 1] +y3_tidal_degree2_interface7 = y3[2, 6] +y3_tidal_degree2_interface10 = y3[2, 9] +y3_tidal_degree2_interface11 = y3[2, 10] + +y4_tidal_degree2_interface1 = y4[2, 0] +y4_tidal_degree2_interface2 = y4[2, 1] +y4_tidal_degree2_interface7 = y4[2, 6] +y4_tidal_degree2_interface10 = y4[2, 9] +y4_tidal_degree2_interface11 = y4[2, 10] + +y5_tidal_degree2_interface1 = y5[2, 0] +y5_tidal_degree2_interface2 = y5[2, 1] +y5_tidal_degree2_interface7 = y5[2, 6] +y5_tidal_degree2_interface10 = y5[2, 9] +y5_tidal_degree2_interface11 = y5[2, 10] + +y6_tidal_degree2_interface1 = y6[2, 0] +y6_tidal_degree2_interface2 = y6[2, 1] +y6_tidal_degree2_interface7 = y6[2, 6] +y6_tidal_degree2_interface10 = y6[2, 9] +y6_tidal_degree2_interface11 = y6[2, 10] +#}}} + +# For surface load +md.love.forcing_type = 11 + +md = solve(md,'lv') + +# Extract love kernels #{{{ +y1 = md.results.LoveSolution.LoveKernelsReal[:,0,:,0].squeeze() +y2 = md.results.LoveSolution.LoveKernelsReal[:,0,:,1].squeeze() +y3 = md.results.LoveSolution.LoveKernelsReal[:,0,:,2].squeeze() +y4 = md.results.LoveSolution.LoveKernelsReal[:,0,:,3].squeeze() +y5 = md.results.LoveSolution.LoveKernelsReal[:,0,:,4].squeeze() +y6 = md.results.LoveSolution.LoveKernelsReal[:,0,:,5].squeeze() + +y1_load_degree2_interface1 = y1[2, 0] +y1_load_degree2_interface2 = y1[2, 1] +y1_load_degree2_interface7 = y1[2, 6] +y1_load_degree2_interface10 = y1[2, 9] +y1_load_degree2_interface11 = y1[2, 10] + +y2_load_degree2_interface1 = y2[2, 0] +y2_load_degree2_interface2 = y2[2, 1] +y2_load_degree2_interface7 = y2[2, 6] +y2_load_degree2_interface10 = y2[2, 9] +y2_load_degree2_interface11 = y2[2, 10] + +y3_load_degree2_interface1 = y3[2, 0] +y3_load_degree2_interface2 = y3[2, 1] +y3_load_degree2_interface7 = y3[2, 6] +y3_load_degree2_interface10 = y3[2, 9] +y3_load_degree2_interface11 = y3[2, 10] + +y4_load_degree2_interface1 = y4[2, 0] +y4_load_degree2_interface2 = y4[2, 1] +y4_load_degree2_interface7 = y4[2, 6] +y4_load_degree2_interface10 = y4[2, 9] +y4_load_degree2_interface11 = y4[2, 10] + +y5_load_degree2_interface1 = y5[2, 0] +y5_load_degree2_interface2 = y5[2, 1] +y5_load_degree2_interface7 = y5[2, 6] +y5_load_degree2_interface10 = y5[2, 9] +y5_load_degree2_interface11 = y5[2, 10] + +y6_load_degree2_interface1 = y6[2, 0] +y6_load_degree2_interface2 = y6[2, 1] +y6_load_degree2_interface7 = y6[2, 6] +y6_load_degree2_interface10 = y6[2, 9] +y6_load_degree2_interface11 = y6[2, 10] +#}}} + + #Fields and tolerances to track changes #loading love numbers -field_names = ['LoveH_loading_elastic', 'LoveK_loading_elastic', 'LoveL_loading_elastic', 'LoveKernels_degree1', 'LoveKernels_degree2'] -field_tolerances = [1e-10, 1e-10, 1e-10, 1e-10, 1e-10] -field_values = [md.results.LoveSolution.LoveHr[:, 0], - md.results.LoveSolution.LoveKr[:, 0], - md.results.LoveSolution.LoveLr[:, 0], - md.results.LoveSolution.LoveKernelsReal[1, 0, :, :], - md.results.LoveSolution.LoveKernelsReal[2, 0, :, :]] +field_names = [ + 'y1_tidal_degree2_interface1', 'y1_tidal_degree2_interface2', 'y1_tidal_degree2_interface7', 'y1_tidal_degree2_interface10', 'y1_tidal_degree2_interface11', + 'y2_tidal_degree2_interface1', 'y2_tidal_degree2_interface2', 'y2_tidal_degree2_interface7', 'y2_tidal_degree2_interface10', 'y2_tidal_degree2_interface11', + 'y3_tidal_degree2_interface1', 'y3_tidal_degree2_interface2', 'y3_tidal_degree2_interface7', 'y3_tidal_degree2_interface10', 'y3_tidal_degree2_interface11', + 'y4_tidal_degree2_interface1', 'y4_tidal_degree2_interface2', 'y4_tidal_degree2_interface7', 'y4_tidal_degree2_interface10', 'y4_tidal_degree2_interface11', + 'y5_tidal_degree2_interface1', 'y5_tidal_degree2_interface2', 'y5_tidal_degree2_interface7', 'y5_tidal_degree2_interface10', 'y5_tidal_degree2_interface11', + 'y6_tidal_degree2_interface1', 'y6_tidal_degree2_interface2', 'y6_tidal_degree2_interface7', 'y6_tidal_degree2_interface10', 'y6_tidal_degree2_interface11', + 'y1_load_degree2_interface1', 'y1_load_degree2_interface2', 'y1_load_degree2_interface7', 'y1_load_degree2_interface10', 'y1_load_degree2_interface11', + 'y2_load_degree2_interface1', 'y2_load_degree2_interface2', 'y2_load_degree2_interface7', 'y2_load_degree2_interface10', 'y2_load_degree2_interface11', + 'y3_load_degree2_interface1', 'y3_load_degree2_interface2', 'y3_load_degree2_interface7', 'y3_load_degree2_interface10', 'y3_load_degree2_interface11', + 'y4_load_degree2_interface1', 'y4_load_degree2_interface2', 'y4_load_degree2_interface7', 'y4_load_degree2_interface10', 'y4_load_degree2_interface11', + 'y5_load_degree2_interface1', 'y5_load_degree2_interface2', 'y5_load_degree2_interface7', 'y5_load_degree2_interface10', 'y5_load_degree2_interface11', + 'y6_load_degree2_interface1', 'y6_load_degree2_interface2', 'y6_load_degree2_interface7', 'y6_load_degree2_interface10', 'y6_load_degree2_interface11' + ] +field_tolerances = [ + 9e-8, 1e-7, 3e-7, 3e-7, 3e-7, + 9e-8, 1e-7, 3e-7, 3e-7, 1e-10, + 9e-8, 8e-8, 2e-8, 2e-7, 4e-7, + 9e-8, 9e-8, 2e-7, 4e-7, 1e-10, + 4e-7, 4e-7, 2e-7, 3e-8, 2e-8, + 2e-5, 2e-6, 2e-6, 1e-6, 2e-7, + 3e-6, 3e-6, 3e-6, 4e-6, 4e-6, + 3e-6, 3e-6, 2e-6, 6e-7, 1e-10, + 3e-6, 3e-6, 5e-7, 3e-6, 5e-6, + 3e-6, 3e-6, 9e-7, 7e-7, 1e-10, + 4e-6, 4e-6, 3e-6, 5e-7, 3e-7, + 3e-6, 3e-6, 2e-6, 7e-7, 2e-7 + ] +field_values = [ + y1_tidal_degree2_interface1, y1_tidal_degree2_interface2, y1_tidal_degree2_interface7, y1_tidal_degree2_interface10,y1_tidal_degree2_interface11, + y2_tidal_degree2_interface1, y2_tidal_degree2_interface2, y2_tidal_degree2_interface7, y2_tidal_degree2_interface10,y2_tidal_degree2_interface11, + y3_tidal_degree2_interface1, y3_tidal_degree2_interface2, y3_tidal_degree2_interface7, y3_tidal_degree2_interface10,y3_tidal_degree2_interface11, + y4_tidal_degree2_interface1, y4_tidal_degree2_interface2, y4_tidal_degree2_interface7, y4_tidal_degree2_interface10,y4_tidal_degree2_interface11, + y5_tidal_degree2_interface1, y5_tidal_degree2_interface2, y5_tidal_degree2_interface7, y5_tidal_degree2_interface10,y5_tidal_degree2_interface11, + y6_tidal_degree2_interface1, y6_tidal_degree2_interface2, y6_tidal_degree2_interface7, y6_tidal_degree2_interface10,y6_tidal_degree2_interface11, + y1_load_degree2_interface1, y1_load_degree2_interface2, y1_load_degree2_interface7, y1_load_degree2_interface10,y1_load_degree2_interface11, + y2_load_degree2_interface1, y2_load_degree2_interface2, y2_load_degree2_interface7, y2_load_degree2_interface10,y2_load_degree2_interface11, + y3_load_degree2_interface1, y3_load_degree2_interface2, y3_load_degree2_interface7, y3_load_degree2_interface10,y3_load_degree2_interface11, + y4_load_degree2_interface1, y4_load_degree2_interface2, y4_load_degree2_interface7, y4_load_degree2_interface10,y4_load_degree2_interface11, + y5_load_degree2_interface1, y5_load_degree2_interface2, y5_load_degree2_interface7, y5_load_degree2_interface10,y5_load_degree2_interface11, + y6_load_degree2_interface1, y6_load_degree2_interface2, y6_load_degree2_interface7, y6_load_degree2_interface10,y6_load_degree2_interface11 + ] Index: ../trunk-jpl/test/NightlyRun/test2052.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test2052.m (revision 25302) +++ ../trunk-jpl/test/NightlyRun/test2052.m (revision 25303) @@ -30,8 +30,6 @@ %% solve for GIA deflection md=solve(md,'Gia'); -pause - %Test Name: GiaIvinsBenchmarksAB2dC1 U_AB2dC1 = md.results.GiaSolution.UGia; URate_AB2dC1 = md.results.GiaSolution.UGiaRate; Index: ../trunk-jpl/src/m/classes/geometry.m =================================================================== --- ../trunk-jpl/src/m/classes/geometry.m (revision 25302) +++ ../trunk-jpl/src/m/classes/geometry.m (revision 25303) @@ -89,11 +89,8 @@ end % }}} function marshall(self,prefix,md,fid) % {{{ - disp(md.geometry.thickness) WriteData(fid,prefix,'object',self,'fieldname','surface','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - disp(md.geometry.thickness) - pause WriteData(fid,prefix,'object',self,'fieldname','base','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','bed','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','hydrostatic_ratio','format','DoubleMat','mattype',1); Index: ../trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- ../trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 25302) +++ ../trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 25303) @@ -247,34 +247,34 @@ mu0 = md.love.mu0 rr = md.materials.radius rho = md.materials.density - rho_avg = (rho * np.diff(np.power(rr, 3), n=1, axis=0)) / np.diff(np.power(rr, 3).sum()).sum() + rho_avg_partial = np.diff(np.power(rr, 3), n=1, axis=0) + rho_avg = ((rho * rho_avg_partial) / rho_avg_partial.sum()).sum() temp_field = np.empty((degmax + 1, nfreq, nlayer + 1, 6)) temp_field.fill(0.0) for ii in range(degmax + 1): for jj in range(nfreq): for kk in range(nlayer + 1): - if kk < nlayer + 1: + if kk < nlayer: # NOTE: Upper bound of range is non-inclusive (compare to src/m/solve/parseresultsfromdisk.m) ll = ii * (nlayer + 1) * 6 + (kk * 6 + 1) + 3 - temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0 # mm = 4 - temp_field[ii, jj, kk, 1] = field[ll + (2 - 1), jj] * mu0 # mm = 5 - temp_field[ii, jj, kk, 2] = field[ll + (3 - 1), jj] * r0 # mm = 6 - temp_field[ii, jj, kk, 3] = field[ll + (4 - 1), jj] * mu0 # mm = 1 - temp_field[ii, jj, kk, 4] = field[ll + (5 - 1), jj] * r0 * g0 # mm = 2 - temp_field[ii, jj, kk, 5] = field[ll + (6 - 1), jj] * g0 # mm = 3 - print(temp_field) + temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0 # mm = 4 + temp_field[ii, jj, kk, 1] = field[ll + (1 - 1), jj] * mu0 # mm = 5 + temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0 # mm = 6 + temp_field[ii, jj, kk, 3] = field[ll + (3 - 1), jj] * mu0 # mm = 1 + temp_field[ii, jj, kk, 4] = field[ll + (4 - 1), jj] * r0 * g0 # mm = 2 + temp_field[ii, jj, kk, 5] = field[ll + (5 - 1), jj] * g0 # mm = 3 else: # surface - ll = ii * (nlayer + 1) * 6 - 2 - temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0 - temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0 - temp_field[ii, jj, kk, 4] = field[ll + (3 - 1), jj] * r0 * g0 + ll = (ii + 1) * (nlayer + 1) * 6 - 2 + temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0 + temp_field[ii, jj, kk, 2] = field[ll + (1 - 1), jj] * r0 + temp_field[ii, jj, kk, 4] = field[ll + (2 - 1), jj] * r0 * g0 # surface BC temp_field[ii, jj, kk, 3] = 0 if md.love.forcing_type == 9: temp_field[ii, jj, kk, 1] = 0 - temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0 + temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0 elif md.love.forcing_type == 11: - temp_field[ii, jj, kk, 1] = -(2 * (ii - 1) + 1) * rho_avg / 3 - temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0 + temp_field[ii, jj, kk, 1] = -(2 * ii + 1) * rho_avg / 3 + temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0 field = temp_field if time != -9999: Index: ../trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- ../trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 25302) +++ ../trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 25303) @@ -228,15 +228,15 @@ g0 = md.love.g0; mu0 = md.love.mu0; rr = md.materials.radius; - rho= md.materials.density; - rho_avg = sum( rho .* diff(rr.^3)/sum(diff(rr.^3)) ); - temp_field = cell(degmax+1,nfreq,nlayer+1,6); + rho= md.materials.density; + rho_avg = sum(rho.*diff(rr.^3)/sum(diff(rr.^3))); + temp_field = cell(degmax+1,nfreq,nlayer+1,6); for ii=1:degmax+1 for jj=1:nfreq for kk=1:nlayer+1 if (kk