%% covert mba from figure 5a to crustal thickness
% 2022-06-18 Tingting Zheng

clc
close all
clear all

%% plot distance vs. thermal & mba, left panel 
%%%%%%%%%%%%% original code in this section is from Hans %%%%%%%%%
da=0.01;
age=-1000:1000;
abat=sqrt(abs(da*age));
dt=10;
dg=110;

fid_f1=figure
figname=strcat('Fig4_filter_model.pdf');
%%%%%%%%%%%%%%%%%%%%%%%
subplot(441)

lon=-100;

N=141;
hgate1=fix(N/2)/20;

h=70;

hold on
gate_scale=[-70,-70,30,30];
fill(hgate1*[-1,1,1,-1],gate_scale,.95*[1,1,1],'LineStyle','none')
plot(age/dt,abat*dg-h,'-r','LineWidth',1)
gate=1/(N)*ones(1,N);
convabat=conv((abat),gate,'same');
plot(age/dt,convabat*dg*1.5-h-35,'-b','LineWidth',1)
title(strcat('a. ',num2str(lon)));
set(gca,'linewidth',1,'xticklabel',[]);
% legend('filter width','observed thermal','observed MBA','thermal model ',' flattened mba')
% axis ij
hold off
xlim([-5 5])
xticks(-5:1:5)
yticks(-60:20:20)
ylabel('Filtered model');

y_scale=[-60 20];
ylim(y_scale)

model_thermal=abat*dg-h;
model_mba=convabat*dg-h;
model_rmba1=model_mba-model_thermal;

%%
%%%% plot rmba
subplot(445)
hold on
fill(hgate1*[-1,1,1,-1],[-20,-20,60,60],.95*[1,1,1],'LineStyle','none');
plot(age/dt,model_rmba1,'-b','LineWidth',1)
set(gca,'linewidth',1,'xticklabel',[]);
ylabel('RMBA (mGal)')
xlim([-5 5])
xticks(-5:1:5)
ylim([-20 60])

%%%%%%%%%%%%%%%%%%%%%
subplot(442)
lon=-99.7;

N=101;
hgate2=fix(N/2)/20;

hold on
fill(hgate2*[-1,1,1,-1],gate_scale,.95*[1,1,1],'LineStyle','none')
plot(age/dt,abat*dg-h,'-r','LineWidth',1)
gate=1/(N)*ones(1,N);
convabat=conv(abat,gate,'same');
plot(age/dt,1.5*convabat*dg-h-30,'-b','LineWidth',1)
title(strcat('b. ',num2str(lon)));
set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);
xlim([-5 5])
xticks(-5:1:5)
yticks(-60:20:20)
ylim(y_scale)
hold off

model_thermal=abat*dg-h;
model_mba=convabat*dg-h;
model_rmba2=model_mba-model_thermal;

subplot(446)
hold on
fill(hgate2*[-1,1,1,-1],[-20,-20,60,60],.95*[1,1,1],'LineStyle','none')
plot(age/dt,model_rmba2,'-b','LineWidth',1)
set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);
xlim([-5 5])
xticks(-5:1:5)
ylim([-20 60])

%%%%%%%%%%%%%%%%%%%%%%
subplot(443)
lon=-98.5;

N=61;
hgate3=fix(N/2)/20;


fill(hgate3*[-1,1,1,-1],gate_scale,.95*[1,1,1],'LineStyle','none')
hold on
plot(age/dt,abat*dg-h,'-r','LineWidth',1)
gate=1/(N)*ones(1,N);
convabat=conv(abat,gate,'same');
plot(age/dt,convabat*dg-h+1,'-b','LineWidth',1)
hold off
xlim([-5 5])
xticks(-5:1:5)
yticks(-60:20:20)
ylim(y_scale)
title(strcat('c. ',num2str(lon)));
set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);

model_thermal=abat*dg-h;
model_mba=convabat*dg-h;
model_rmba3=model_mba-model_thermal;

subplot(447)
hold on
fill(hgate3*[-1,1,1,-1],[-20,-20,60,60],.95*[1,1,1],'LineStyle','none')
plot(age/dt,model_rmba3,'-b','LineWidth',1)
set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);
xlim([-5 5])
xticks(-5:1:5)
ylim([-20 60])

%%%%%%%%%%%%%%%%%
subplot(444)

lon=-97;

N=21;
hgate4=fix(N/2)/20;

fill(hgate4*[-1,1,1,-1],gate_scale,.95*[1,1,1],'LineStyle','none')
hold on
plot(age/dt,abat*dg-h,'-r','LineWidth',1);
gate=1/(N)*ones(1,N);
convabat=conv(abat,gate,'same');
plot(age/dt,convabat*dg-h,'-b','LineWidth',1)
hold off
xlim([-5 5])
xticks(-5:1:5)
yticks(-60:20:20)
ylim(y_scale)
title(strcat('d. ',num2str(lon)));;

set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);

model_thermal=abat*dg-h;
model_mba=convabat*dg-h;
model_rmba4=model_mba-model_thermal;

subplot(448)
hold on
fill(hgate4*[-1,1,1,-1],[-20,-20,60,60],.95*[1,1,1],'LineStyle','none');
plot(age/dt,model_rmba4,'-b','LineWidth',1);

set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);
xlim([-5 5])
xticks(-5:1:5)
ylim([-20 60])

temp_aa=[hgate1 hgate2 hgate3 hgate4];
%% cut rmba data
temp=[age;model_rmba1;model_rmba2;model_rmba3;model_rmba4]';

%% plot RMBA vs misfit width
temp1=temp(temp(:,1)>-100 & temp(:,1)<100,:);
temp_aa=[hgate1 hgate2 hgate3 hgate4];
temp_bb=[max(temp1(:,2)) max(temp1(:,3)) max(temp1(:,4)) max(temp1(:,5))]


temp1=temp;
dis=temp1(:,1)/dt;
rmba1=temp1(:,2);
rmba2=temp1(:,3);
rmba3=temp1(:,4);
rmba4=temp1(:,5);

%% 1-D rmba to crust 
sdata1=rmba1;
sdata2=rmba2;
sdata3=rmba3;
sdata4=rmba4;
nx=length(sdata1);
dx=max(dis)-min(dis);

mm=find(isnan(sdata1));
sdata1(mm)=0;

%%%%%%%%set parameters
addm=10;            %%%%distance for downward continuation, km
density=0.6;     %%%g/cm^3
wlarge=135; % large wavelength of the taper, km
wsmall=50;  % small wavelength of the taper, km

%---Normalize data to new reference level-----

%------Mean of Min & Max----------
zmax=-1e10;zmin=1e10;
zmax=max(max(max(zmax,sdata1)));
zmin=min(min(min(zmin,sdata1)));
slev1=(zmin+zmax)/2;

zmax=max(max(max(zmax,sdata2)));
zmin=min(min(min(zmin,sdata2)));
slev2=(zmin+zmax)/2;

zmax=max(max(max(zmax,sdata3)));
zmin=min(min(min(zmin,sdata3)));
slev3=(zmin+zmax)/2;

zmax=max(max(max(zmax,sdata4)));
zmin=min(min(min(zmin,sdata4)));
slev4=(zmin+zmax)/2;

sdata1=-0.001*(sdata1-slev1);
sdata2=-0.001*(sdata2-slev2);
sdata3=-0.001*(sdata3-slev3);
sdata4=-0.001*(sdata4-slev4);

dx=100000*dx;
addm=addm*100000;

wlarge=wlarge*100000;
wsmall=wsmall*100000;
kcut1=2*pi/wlarge;
kcut2=2*pi/wsmall;
dkcut=kcut2-kcut1;

G=6.673*10^-8;
grav=2*pi*G*density;
topo=1./grav;

%%%%
dkx = 2.*pi./(nx.*dx); %%    dkx - sample interval in the kx direction.
% dky = 2.*pi./(ny.*dy);
ifsdata1 = fft(sdata1);
ifsdata2 = fft(sdata2);
ifsdata3 = fft(sdata3);
ifsdata4 = fft(sdata4);

for j = 1:nx
    %    nx - dimension of grid in ky direction (a power of two).
    %  Output parameters:
    %    kx - the wavenumber coordinate in the kx direction.
    
    nyqx=nx/2+1;
    
    if j <= nyqx
        kx=(j-1)*dkx; %kx - the wavenumber coordinate in the kx direction.
    else
        kx=(j-nx-1)*dkx;
    end
    k = sqrt(kx.^2);
    
    weight=1;
    if k>kcut2
        weight=0;
    elseif k>kcut1
        t=(k-kcut1)/dkcut*pi;
        weight=(cos(t)+1)/2;
    end
    csdata1(j) = ifsdata1(j).*topo.*(exp(k.*addm)).*weight;
    csdata2(j) = ifsdata2(j).*topo.*(exp(k.*addm)).*weight;
    csdata3(j) = ifsdata3(j).*topo.*(exp(k.*addm)).*weight;
    csdata4(j) = ifsdata4(j).*topo.*(exp(k.*addm)).*weight;
    
end
gsdata1 = ifft(csdata1)./100000;
gsdata2 = ifft(csdata2)./100000;
gsdata3 = ifft(csdata3)./100000;
gsdata4 = ifft(csdata4)./100000;

gsdata(mm)=NaN;
%% plot crustal thickness
d=8.5;
y_lim=[0 10];
box=[0,0,10,10];


subplot(449)
hold on
fill(hgate1*[-1,1,1,-1],box,.95*[1,1,1],'LineStyle','none');
plot(dis,gsdata1+d,'-r','LineWidth',1);
xlim([-5 5])
xticks(-5:1:5)
ylim(y_lim)
yticks(0:2:10)
set(gca,'linewidth',1,'xticklabel',[]);
ylabel('Crust (km)')
axis ij
hold off

subplot(4,4,10)
hold on
fill(hgate2*[-1,1,1,-1],box,.95*[1,1,1],'LineStyle','none');
plot(dis,gsdata2+d,'-r','LineWidth',1);
xlim([-5 5])
xticks(-5:1:5)
ylim(y_lim)
yticks(0:2:10)
set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);
axis ij
hold off

subplot(4,4,11)
hold on
fill(hgate3*[-1,1,1,-1],box,.95*[1,1,1],'LineStyle','none');
plot(dis,gsdata3+d,'-r','LineWidth',1);
xlim([-5 5])
xticks(-5:1:5)
ylim(y_lim)
yticks(0:2:10)
set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);
axis ij
hold off

subplot(4,4,12)
hold on
fill(hgate4*[-1,1,1,-1],box,.95*[1,1,1],'LineStyle','none');
plot(dis,gsdata4+d,'-r','LineWidth',1);
xlim([-5 5])
xticks(-5:1:5)
ylim(y_lim)
yticks(0:2:10)
set(gca,'linewidth',1,'xticklabel',[],'yticklabel',[]);
axis ij
hold off

%%%%%%%%%%% plot 4th array %%%%%%%%%%%%%%%
subplot(4,4,13)
hold on
fill(hgate1*[-1,1,1,-1],[0,0,100,100],.95*[1,1,1],'LineStyle','none');
xlim([-5 5])
xticks(-5:1:5)
ylim([1 100])
ylabel('Mantle (km)');
axis ij
hold off

subplot(4,4,14)
hold on
fill(hgate2*[-1,1,1,-1],[0,0,100,100],.95*[1,1,1],'LineStyle','none');
xlim([-5 5])
xticks(-5:1:5)
ylim([1 100])
set(gca,'linewidth',1,'yticklabel',[]);
axis ij
xlabel('Age (Myr, south to north)')


subplot(4,4,15)
hold on
fill(hgate3*[-1,1,1,-1],[0,0,100,100],.95*[1,1,1],'LineStyle','none');
xlim([-5 5])
xticks(-5:1:5)
ylim([1 100])
set(gca,'linewidth',1,'yticklabel',[]);
axis ij
hold off

subplot(4,4,16)
hold on
fill(hgate4*[-1,1,1,-1],[0,0,100,100],.95*[1,1,1],'LineStyle','none');
xlim([-5 5])
xticks(-5:1:5)
ylim([1 100])
set(gca,'linewidth',1,'yticklabel',[]);
axis ij
hold off

saveas(fid_f1,figname,'pdf');