Index: /issm/trunk-jpl/src/m/boundaryconditions/love_numbers.m =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/love_numbers.m (revision 25137) +++ /issm/trunk-jpl/src/m/boundaryconditions/love_numbers.m (revision 25138) @@ -1,3 +1,3 @@ -function series=love_numbers(type,reference_frame) +function series=love_numbers(varargin) %LOVE_NUMBERS: provide love numbers % retrieved from: http://www.srosat.com/iag-jsg/loveNb.php @@ -5,34 +5,35 @@ % % Usage: -% series = love_numbers(type) -% series = love_numbers(type,reference_frame) +% series = love_numbers('type','loadingverticaldisplacement','referenceframe','CM','maxdeg',1001) % -% - type = one of 'h','k','l','gamma' and 'lambda'. +% - type = one of 'loadingverticaldisplacement', 'loadinggravitationalpotential', 'loadinghorizontaldisplacement', ... +% 'tidalverticaldisplacement', 'tidalgravitationalpotential', 'tidalhorizontaldisplacement' % - reference_frame = one of 'CM' (default) and 'CF'. +% - maxdeg = default 1001 % % Example: -% love_k = love_numbers('k'); -% love_k = love_numbers('k','CF'); -% +% h=love_number('type','loadingverticaldisplacement','referenceframe','CM','maxdeg',maxdeg); +% k=love_number('type','loadinggravitationalpotential','referenceframe','CM','maxdeg',maxdeg); +% l=love_number('type','loadinghorizontaldisplacement','referenceframe','CM','maxdeg',maxdeg); +% th=love_number('type','tidalverticaldisplacement','referenceframe','CM','maxdeg',maxdeg); +% tk=love_number('type','tidalgravitationalpotential','referenceframe','CM','maxdeg',maxdeg); +% tl=love_number('type','tidalhorizontaldisplacement','referenceframe','CM','maxdeg',maxdeg); -% some checks: -if nargin==1, - frame='CM'; - disp('Info: computation is done in Center of Mass (CM) reference frame'); -elseif nargin==2, - frame=reference_frame; - if ~( strcmpi(reference_frame,'CM') | strcmpi(reference_frame,'CF')), - error('reference_frame should be one of ''CM'' or ''CF'''); - end -else - help love_numbers - error('bad usage'); +%recover options: +options=pairoptions(varargin{:}); +type=getfieldvalue(options,'type'); +frame=getfieldvalue(options,'referenceframe','CM'); +maxdeg=getfieldvalue(options,'maxdeg',1000); + +if(maxdeg>10000), + error('PREM love numbers computed only for deg < 10,000. Request lower maxdeg' ); end -if ~( strcmpi(type,'h') | strcmpi(type,'k') | strcmpi(type,'l') | strcmpi(type,'gamma') | strcmpi(type,'lambda') ), - error('type should be one of ''h'',''k'',''l'',''gamma'' and ''lambda'''); +if ~( strcmpi(type,'loadingverticaldisplacement') | strcmpi(type,'loadinggravitationalpotential') | strcmpi(type,'loadinghorizontaldisplacement') ... + | strcmpi(type,'tidalverticaldisplacement') | strcmpi(type,'tidalgravitationalpotential') | strcmpi(type,'tidalhorizontaldisplacement')), + error('type should be one of ''loadingverticaldisplacement'', ''loadinggravitationalpotential'', ''loadinghorizontaldisplacement'',''tidalverticaldisplacement'', ''tidalgravitationalpotential'', ''tidalhorizontaldisplacement'''); end - love_numbers=[... + love_numbers=[... 0 0 0 0 0 0 0 -1.28740059 -1.00000000 -0.89858519 1.28740059 0.42519882 0.89858519 0.00000000 @@ -10037,14 +10038,20 @@ -6.27342778 -0.00030945 0.00018956 7.27311833 0.99905480 0.99950099 0.49327194]; - if type=='h', + %cut off series at maxdeg: + love_numbers(maxdeg+2:end,:)=[]; + + %retrive right type: + if strcmpi(type,'loadingverticaldisplacement'), series=love_numbers(:,1); - elseif type=='k', + elseif strcmpi(type,'loadinggravitationalpotential'), series=love_numbers(:,2); - elseif type=='l', + elseif strcmpi(type,'loadinghorizontaldisplacement'), series=love_numbers(:,3); - elseif type=='gamma', + elseif strcmpi(type,'tidalverticaldisplacement'), series=love_numbers(:,4); - elseif type=='lambda', + elseif strcmpi(type,'tidalgravitationalpotential'), series=love_numbers(:,5); + elseif strcmpi(type,'tidalhorizontaldisplacement'), + series=love_numbers(:,6); else error(['love_numbers error message: unknow type:' type]); @@ -10052,12 +10059,12 @@ % choose degree 1 term for CF reference system - if frame=='CM', + if strcmpi(frame,'CM'), return; - elseif frame=='CF', % from Blewitt, 2003, JGR - if type=='h', + elseif strcmpi(frame,'CF'), % from Blewitt, 2003, JGR + if strcmpi(type,'loadingverticaldisplacement'), series(2,1) = -0.269; - elseif type=='k', + elseif strcmpi(type,'loadinggravitationalpotential'), series(2,1) = 0.021; - elseif type=='l', + elseif strcmpi(type,'loadinghorizontaldisplacement'), series(2,1) = 0.134; end @@ -10066,2 +10073,4 @@ end + + Index: /issm/trunk-jpl/src/m/classes/love.m =================================================================== --- /issm/trunk-jpl/src/m/classes/love.m (revision 25138) +++ /issm/trunk-jpl/src/m/classes/love.m (revision 25138) @@ -0,0 +1,100 @@ +%LOVE numbers class definition +% +% Usage: +% love=love() %will setup love numbers deg 1001 by default +% love=love('maxdeg',10001); %supply numbers of degrees required (here, 10001) + +classdef love + properties (SetAccess=public) + + %regular love numbers: + h = []; %provided by PREM model + k = []; %idem + l = []; %idem + + %tidal love numbers for computing rotational feedback: + th = []; + tk = []; + tl = []; + tk2secular = 0; %deg 2 secular number. + + end + methods + function self = love(varargin) % {{{ + options=pairoptions(varargin{:}); + maxdeg=getfieldvalue(options,'maxdeg',1000); + self=setdefaultparameters(self,maxdeg); + end % }}} + function self = setdefaultparameters(self,maxdeg) % {{{ + + %initialize love numbers: + self.h=love_numbers('type','loadingverticaldisplacement','referenceframe','CM','maxdeg',maxdeg); + self.k=love_numbers('type','loadinggravitationalpotential','referenceframe','CM','maxdeg',maxdeg); + self.l=love_numbers('type','loadinghorizontaldisplacement','referenceframe','CM','maxdeg',maxdeg); + self.th=love_numbers('type','tidalverticaldisplacement','referenceframe','CM','maxdeg',maxdeg); + self.tk=love_numbers('type','tidalgravitationalpotential','referenceframe','CM','maxdeg',maxdeg); + self.tl=love_numbers('type','tidalhorizontaldisplacement','referenceframe','CM','maxdeg',maxdeg); + + %secular fluid love number: + self.tk2secular=0.942; + + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.islove==0), + return; + end + + md = checkfield(md,'fieldname','solidearth.love.h','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','solidearth.love.k','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','solidearth.love.l','NaN',1,'Inf',1); + + + md = checkfield(md,'fieldname','solidearth.love.th','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','solidearth.love.tk','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','solidearth.love.tk2secular','NaN',1,'Inf',1); + + %check that love numbers are provided at the same level of accuracy: + if (size(self.h,1)~=size(self.k,1) | size(self.h,1)~=size(self.l,1)), + error('love error message: love numbers should be provided at the same level of accuracy'); + end + + + end % }}} + function list=defaultoutputs(self,md) % {{{ + list = {}; + end % }}} + function disp(self) % {{{ + disp(sprintf(' love parameters:')); + + fielddisplay(self,'h','load Love number for radial displacement'); + fielddisplay(self,'k','load Love number for gravitational potential perturbation'); + fielddisplay(self,'l','load Love number for horizontal displacements'); + + fielddisplay(self,'th','tidal load Love number (deg 2)'); + fielddisplay(self,'tk','tidal load Love number (deg 2)'); + fielddisplay(self,'tk2secular','secular fluid Love number'); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + WriteData(fid,prefix,'object',self,'fieldname','h','name','md.solidearth.love.h','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','k','name','md.solidearth.love.k','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','l','name','md.solidearth.love.l','format','DoubleMat','mattype',1); + + WriteData(fid,prefix,'object',self,'fieldname','th','name','md.solidearth.love.th','format','DoubleMat'); + WriteData(fid,prefix,'object',self,'fieldname','tk','name','md.solidearth.love.tk','format','DoubleMat'); + WriteData(fid,prefix,'object',self,'fieldname','tl','name','md.solidearth.love.tl','format','DoubleMat'); + WriteData(fid,prefix,'object',self,'fieldname','tk2secular','name','md.solidearth.love.tk2secular','format','Double'); + + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + writejs1Darray(fid,[modelname '.love.h'],self.h); + writejs1Darray(fid,[modelname '.love.k'],self.k); + writejs1Darray(fid,[modelname '.love.l'],self.l); + end % }}} + function self = extrude(self,md) % {{{ + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25137) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25138) @@ -9,5 +9,5 @@ settings = solidearthsettings(); surfaceload = surfaceload(); - love = loadinglove(); + love = love(); rotational = rotational(); planetradius = planetradius('earth'); Index: sm/trunk-jpl/src/m/classes/tidallove.m =================================================================== --- /issm/trunk-jpl/src/m/classes/tidallove.m (revision 25137) +++ (revision ) @@ -1,67 +1,0 @@ -%TIDALLOVE class definition -% -% Usage: -% tidallove=tidallove(); - -classdef tidallove - properties (SetAccess=public) - h = 0; - k = 0; - k2secular = 0; - end - methods - function self = tidallove(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - otherwise - error('constructor not supported'); - end - end % }}} - function self = setdefaultparameters(self) % {{{ - - %tidal love numbers: - self.h=0.6149; %degree 2 - self.k=0.3055; % degree 2 - - %secular fluid love number: - self.k2secular=0.942; - - end % }}} - function md = checkconsistency(self,md,solution,analyses) % {{{ - - if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.istidallove==0), - return; - end - - md = checkfield(md,'fieldname','tidallove.h','NaN',1,'Inf',1); - md = checkfield(md,'fieldname','tidallove.k','NaN',1,'Inf',1); - md = checkfield(md,'fieldname','tidallove.k2secular','NaN',1,'Inf',1); - - end % }}} - function list=defaultoutputs(self,md) % {{{ - list = {}; - end % }}} - function disp(self) % {{{ - disp(sprintf(' tidallove parameters:')); - - fielddisplay(self,'h','tidal load Love number (deg 2)'); - fielddisplay(self,'k','tidal load Love number (deg 2)'); - fielddisplay(self,'k2secular','secular fluid Love number'); - - end % }}} - function marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'name','md.solidearth.love.model','data',2,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','h','name','md.solidearth.love.h','format','Double','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','k','name','md.solidearth.love.k','format','Double','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','k2secular','name','md.solidearth.love.k2secular','format','DoubleMat','mattype',1); - end % }}} - function savemodeljs(self,fid,modelname) % {{{ - writejs1Darray(fid,[modelname '.tidallove.h'],self.h); - writejs1Darray(fid,[modelname '.tidallove.k'],self.k); - writejs1Darray(fid,[modelname '.tidallove.k2secular'],self.k2secular); - end % }}} - function self = extrude(self,md) % {{{ - end % }}} - end -end