%Test Name: GiaCaron
%Forward Love number solution for a viscoelastic earth, model M3-L70-V01 from
%Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M., Martinec, Z.,
%Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A. and King, M. A.
%(2011), A benchmark study for glacial isostatic adjustment codes. Geophysical
%Journal International, 185: 106--132. doi:10.1111/j.1365-246X.2011.04952.x

md=model();
md.cluster=generic('name',oshostname(),'np',1);

% set validation=1 for comparing against the Spada benchark.   
validation=0; 

md.materials=materials('litho');
md.miscellaneous.name='FourierLoveTest';
md.groundingline.migration='None';

md.verbose=verbose('all');
cst=365.25*24*3600*1000;

md.materials.numlayers=6;
md.materials.radius =  [10 1222.5 3.4800e+03   5.7010e+03   5.9510e+03   6.3010e+03   6.3710e+03]'*1e3;
md.materials.density=  [1.0750e4 1.0750e+04   4.9780e+03   3.8710e+03   3.4380e+03   3.0370e+03]';
md.materials.lame_mu=  [1e-5         0   2.2834e+00   1.0549e+00   7.0363e-01   5.0605e-01]'*1e11;
md.materials.viscosity=[0            0   2.0000e+00   1.0000e+00   1.0000e+00   1.0000e+25]'*1e21;
md.materials.lame_lambda=md.materials.lame_mu*0+5e14;
md.materials.issolid=[1 0 1 1 1 1]';
md.materials.isburgers=zeros(md.materials.numlayers,1);

md.love.allow_layer_deletion=1;
md.love.frequencies=([0]*2*pi)'/cst;
md.love.nfreq=length(md.love.frequencies);
md.love.sh_nmax=256;

md.materials.burgers_mu=md.materials.lame_mu;
md.materials.burgers_viscosity=md.materials.viscosity;

md=solve(md,'lv');

%Fields and tolerances to track changes
%loading love numbers
field_names     ={'LoveH_loading_elastic','LoveK_loading_elastic','LoveL_loading_elastic'};
field_tolerances={2.0e-8,2.0e-8,2.0e-8};
field_values={...
	(md.results.LoveSolution.LoveHr(:,1)),...
	(md.results.LoveSolution.LoveKr(:,1)),...
	(md.results.LoveSolution.LoveLr(:,1)),...
	};

% validate elastic loading solutions against the Spada benchmark. {{{ 
if validation 
	spada_solutions = load('spada_elastic_loading_deg_h_l_k'); 
	spada_d = spada_solutions(:,1); 
	spada_h = spada_solutions(:,2); 
	spada_l = spada_solutions(:,3); 
	spada_k = spada_solutions(:,4); 

	%rename ISSM solutions.  
	issm_d = [1:md.love.sh_nmax];  
	issm_h = md.results.LoveSolution.LoveHr(2:end,1); 
	issm_l = md.results.LoveSolution.LoveLr(2:end,1); 
	issm_k = md.results.LoveSolution.LoveKr(2:end,1); 

	% relative difference for each degree, except for zero.  
	diff_h = 1 - issm_h./spada_h;  
	diff_l = 1 - issm_l./spada_l;  
	diff_k = 1 - issm_k./spada_k;  

	figure 
	plot(spada_d,[diff_h diff_l diff_k]); grid on; 
	legend('h','l','k'); title('loading love'); 

else
	% 
end 
% }}} 

md.love.frequencies=([1e-3 1e-2 1e-1 1 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst;
md.love.nfreq=length(md.love.frequencies);
md.love.sh_nmax=256;
md.materials.burgers_mu=md.materials.lame_mu;
md.materials.burgers_viscosity=md.materials.viscosity;

md=solve(md,'lv');

%Fields and tolerances to track changes
field_names     ={field_names{:},'LoveH_loading_realpart','LoveK_loading_realpart','LoveL_loading_realpart','LoveH_loading_imagpart','LoveK_loading_imagpart','LoveL_loading_imagpart'};
field_tolerances={field_tolerances{:},5e-7,5e-7,5e-7,5e-7,5e-7,5e-7};
field_values={field_values{:},...
	(md.results.LoveSolution.LoveHr(:,:)),...
	(md.results.LoveSolution.LoveKr(:,:)),...
	(md.results.LoveSolution.LoveLr(:,:)),...
	(md.results.LoveSolution.LoveHi(:,:)),...
	(md.results.LoveSolution.LoveKi(:,:)),...
	(md.results.LoveSolution.LoveLi(:,:)),...
	};

md.love.forcing_type=9;
md.love.sh_nmin=2;
md.love.frequencies=([0 1e-3 1e-2 1e-1 1 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst;
md.love.nfreq=length(md.love.frequencies);
md=solve(md,'lv');

% validate elastic tidal solutions against the Spada benchmark. {{{ 
if validation 
	spada_solutions = load('spada_elastic_tidal_deg_h_l_k'); 
	spada_d = spada_solutions(:,1); 
	spada_h = spada_solutions(:,2); 
	spada_l = spada_solutions(:,3); 
	spada_k = spada_solutions(:,4); 

	%rename ISSM solutions.  
	issm_d = [2:md.love.sh_nmax];  
	issm_h = md.results.LoveSolution.LoveHr(3:end,1); 
	issm_l = md.results.LoveSolution.LoveLr(3:end,1); 
	issm_k = md.results.LoveSolution.LoveKr(3:end,1); 

	% relative difference for each degree, except for zero and one. 
	diff_h = 1 - issm_h./spada_h;  
	diff_l = 1 - issm_l./spada_l;  
	diff_k = 1 - issm_k./spada_k;  

	% k seems to have larger relative difference at higher degrees, although values approach zero. 
	figure 
	plot(spada_d,issm_k./spada_k);

	figure 
	plot(spada_d,[diff_h diff_l diff_k]); grid on; 
	legend('h','l','k'); title('tidal love'); 
else
	% 
end 
% }}} 

%tidal love numbers
field_names     ={field_names{:},'LoveH_tidal_elastic','LoveK_tidal_elastic','LoveL_tidal_elastic','LoveH_tidal_realpart','LoveK_tidal_realpart','LoveL_tidal_realpart','LoveH_tidal_imagpart','LoveK_tidal_imagpart','LoveL_tidal_imagpart'};
field_tolerances={field_tolerances{:},8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6};
field_values={field_values{:},...
	(md.results.LoveSolution.LoveHr(:,1)),...
	(md.results.LoveSolution.LoveKr(:,1)),...
	(md.results.LoveSolution.LoveLr(:,1)),...
	(md.results.LoveSolution.LoveHr(:,2:end)),...
	(md.results.LoveSolution.LoveKr(:,2:end)),...
	(md.results.LoveSolution.LoveLr(:,2:end)),...
	(md.results.LoveSolution.LoveHi(:,2:end)),...
	(md.results.LoveSolution.LoveKi(:,2:end)),...
	(md.results.LoveSolution.LoveLi(:,2:end)),...
	};

%Many layers PREM-based model
%data=load('../Data/PREM_500layers');
%md.love.sh_nmin=1;
%md.materials.radius=data(2:end-2,1);
%md.materials.density=data(3:end-2,2);
%md.materials.lame_lambda=data(3:end-2,3);
%md.materials.lame_mu=data(3:end-2,4);
%md.materials.issolid=data(3:end-2,4)>0;
%ind=find(md.materials.issolid==0);
%md.materials.density(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.density(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)+1).^3);
%md.materials.lame_lambda(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.lame_lambda(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)+1).^3);
%md.materials.lame_mu(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.lame_mu(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)).^3);
%md.materials.radius(ind(2:end)+1)=[];
%md.materials.density(ind(2:end))=[];
%md.materials.lame_lambda(ind(2:end))=[];
%md.materials.lame_mu(ind(2:end))=[];
%md.materials.issolid(ind(2:end))=[];
%md.materials.viscosity=10.^interp1([0 3479e3 3480e3 3680e3 5720e3 5800e3 6270e3 6371e3], log10([1e8 1e8 5e21 1e23 1e22 1e20 1e21 1e40]), md.materials.radius(2:end),'PCHIP');
%md.materials.viscosity=md.materials.viscosity.*md.materials.issolid;
%md.materials.burgers_mu=md.materials.lame_mu;
%md.materials.burgers_viscosity=md.materials.viscosity;
%md.materials.isburgers=md.materials.issolid*0;
%md.love.forcing_type=11;
%md.materials.numlayers=length(md.materials.viscosity);
%md=solve(md,'lv');
%
%field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
%field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
%field_values={field_values{:},...
%	(md.results.LoveSolution.LoveHr(:,1)),...
%	(md.results.LoveSolution.LoveKr(:,1)),...
%	(md.results.LoveSolution.LoveLr(:,1)),...
%	(md.results.LoveSolution.LoveHr(:,2:end)),...
%	(md.results.LoveSolution.LoveKr(:,2:end)),...
%	(md.results.LoveSolution.LoveLr(:,2:end)),...
%	(md.results.LoveSolution.LoveHi(:,2:end)),...
%	(md.results.LoveSolution.LoveKi(:,2:end)),...
%	(md.results.LoveSolution.LoveLi(:,2:end)),...
%	};

%Model VSS96 from Vermeersen, L.L.A., Sabadini, R. & Spada, G., 1996a. Analytical visco-elastic relaxation models, Geophys. Res. Lett., 23, 697–700.
md.materials.radius=[10, 1222.5, 3480., 3600., 3630.5, 3700., 3900., 4000., 4200., 4300., 4500., 4600., 4800., 4900., 5100., 5200., 5400., 5500., 5600.5, 5650., 5701., 5736., 5771.5, 5821., 5951., 5970.5, 6016., 6061., 6150.5, 6151.5, 6251., 6371.]'*1e3;
md.materials.lame_mu=[1e-5, 0., 2.933, 2.8990002, 2.8550003, 2.7340002, 2.675, 2.559, 2.502, 2.388, 2.331, 2.215, 2.157, 2.039, 1.979, 1.8560001, 1.794, 1.73, 1.639, 1.2390001, 1.224, 1.21, 1.128, 0.97700006, 0.906, 0.79, 0.773, 0.741, 0.656, 0.665, 0.602]'*1e11;
md.materials.density=[10925., 10925., 5506.42, 5491.45, 5456.57, 5357.06, 5307.24, 5207.13, 5156.69, 5054.69, 5002.99, 4897.83, 4844.22, 4734.6, 4678.44, 4563.07, 4503.72, 4443.16, 4412.41, 3992.14, 3983.99, 3975.84, 3912.82, 3786.78, 3723.78, 3516.39, 3489.51, 3435.78, 3359.5, 3367.1, 3184.3]';
md.materials.viscosity=[0., 0., 7.999999999999999E+21, 8.5E+21, 8.999999999999999E+21, 3.E+22, 4.E+22, 5.0000000000000004E+22, 6.E+22, 5.0000000000000004E+22, 4.5E+22, 3.E+22, 2.5000000000000002E+22, 1.7999999999999998E+22, 1.3E+22, 7.999999999999999E+21, 6.999999999999999E+21, 6.5E+21, 6.E+21, 5.5E+21, 5.E+21, 4.4999999999999995E+21, 3.9999999999999995E+21, 2.5E+21, 1.9999999999999997E+21, 1.5E+21, 9.999999999999999E+20, 6.E+20, 5.5000000000000007E+20, 2.E+20, 1.E40]';
md.materials.lame_lambda=md.materials.lame_mu*0+5e14;
md.materials.issolid=md.materials.lame_mu>0;
md.materials.numlayers=length(md.materials.lame_mu);
md.materials.burgers_mu=md.materials.lame_mu;
md.materials.burgers_viscosity=md.materials.viscosity;
md.materials.isburgers=md.materials.issolid*0;
md.love.forcing_type=11;
md.love.sh_nmin=1;
md.love.sh_nmax=100;
md=solve(md,'lv');
md.love.frequencies=([0 1e-3 1e-2 1 -1e-3 -1e-2 -1]*2*pi)'/cst;
md.love.nfreq=length(md.love.frequencies);

field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
field_values={field_values{:},...
	(md.results.LoveSolution.LoveHr(:,1)),...
	(md.results.LoveSolution.LoveKr(:,1)),...
	(md.results.LoveSolution.LoveLr(:,1)),...
	(md.results.LoveSolution.LoveHr(:,2:end)),...
	(md.results.LoveSolution.LoveKr(:,2:end)),...
	(md.results.LoveSolution.LoveLr(:,2:end)),...
	(md.results.LoveSolution.LoveHi(:,2:end)),...
	(md.results.LoveSolution.LoveKi(:,2:end)),...
	(md.results.LoveSolution.LoveLi(:,2:end)),...
	};
