clear all
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%importing gravity anomalies
filename_true_grav = fullfile('.', 'tested_data', 'grav_anomaly3.txt');
true_grav=importdata(filename_true_grav);

filename_pred_grav = fullfile('.', 'tested_data', 'grav_pred_data_3.txt');
pred_grav=importdata(filename_pred_grav);

filename_inverted_grav = fullfile('.', 'tested_data', 'grav_inverted_data_3.txt');
inverted_grav=importdata(filename_inverted_grav);


%importing inverted depth
filename_true_depth = fullfile('.', 'tested_data', 'depth_true3.txt');
true_depth=importdata(filename_true_depth);

filename_pred_depth = fullfile('.', 'tested_data', 'depth_predicted3.txt');
pred_depth=importdata(filename_pred_depth);

filename_inv_depth = fullfile('.','tested_data', 'inv_depth3.txt');
depth_inv_all=importdata(filename_inv_depth);
dd_inv=(squeeze(depth_inv_all(94,:,:))).*10^-3;
B = ones(1,1)/1^2;
inverted_depth = conv2(dd_inv,B,'same');
   
filename_lon = fullfile('.', 'tested_data', 'XX_predicted3.txt');
lon=importdata(filename_lon);

filename_lat = fullfile('.', 'tested_data', 'YY_predicted3.txt');
lat=importdata(filename_lat);
%%
figure(1)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
%ccmap2=makecolormap({'white','magenta','red','gold','lightseagreen','darkblue'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
ccmap2=makecolormap({'white','goldenrod','limegreen','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,true_grav,[min(true_grav(:))+1:0.5:max(true_grav(:))-1],'edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Gravity anomalies (mGal)');
m_grid('box','fancy','grid','none','fontsize',14);
title('True Gravity anomalies')

figure(2)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
%ccmap2=makecolormap({'white','magenta','red','gold','lightseagreen','darkblue'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
ccmap2=makecolormap({'white','goldenrod','limegreen','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,pred_grav,[min(true_grav(:))+1:0.5:max(true_grav(:))-1],'edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Gravity anomalies (mGal)');
m_grid('box','fancy','grid','none','fontsize',14);
title('cGAN inverted Gravity anomalies')
figure(3)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
%ccmap2=makecolormap({'white','magenta','red','gold','lightseagreen','darkblue'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
ccmap2=makecolormap({'white','goldenrod','limegreen','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,inverted_grav,[min(true_grav(:))+1:0.5:max(true_grav(:))-1],'edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Gravity anomalies (mGal)');
m_grid('box','fancy','grid','none','fontsize',14);
title('Bott''s inverted Gravity anomalies')


%dd=importdata('C:\Users\Arka\Desktop\south_india_tesseroid\output\final_moho_depth.txt');
true_depth=true_depth.*10^-3;
pred_depth=pred_depth.*10^-3;
figure(4)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
ccmap2=makecolormap({'white','magenta','red','gold','lightseagreen','darkblue'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
%ccmap2=makecolormap({'white','goldenrod','limegreen','olivedrab','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,true_depth,[min(true_depth(:))+1:0.5:max(true_depth(:))-1],'edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Moho depth (km)');

m_grid('box','fancy','grid','none','fontsize',14);
title('True Moho Depth')

figure(5)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
ccmap2=makecolormap({'white','magenta','red','gold','lightseagreen','darkblue'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
%ccmap2=makecolormap({'white','goldenrod','limegreen','olivedrab','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,pred_depth,[min(true_depth(:))+1:0.5:max(true_depth(:))-1],'edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Moho depth (km)');

m_grid('box','fancy','grid','none','fontsize',14);
title('cGAN inverted Moho Depth')

figure(6)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
ccmap2=makecolormap({'white','magenta','red','gold','lightseagreen','darkblue'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
%ccmap2=makecolormap({'white','goldenrod','limegreen','olivedrab','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,inverted_depth,[min(true_depth(:))+1:0.5:max(true_depth(:))-1],'edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Moho depth (km)');

m_grid('box','fancy','grid','none','fontsize',14);
title('Bott''s inverted Moho Depth')

%%
dd_err1=true_depth-pred_depth;
figure(7)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
ccmap2=makecolormap({'blue','white','red'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
%ccmap2=makecolormap({'white','goldenrod','limegreen','olivedrab','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,dd_err1,[-3.5:0.05:3.5],'edgecolor','none');
%ax=m_contfbar(.85,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Error (km)');

m_grid('box','fancy','grid','none','fontsize',14);
title('Error in Moho Depth')
fprintf('predicted depth error is %f\n',(norm(dd_err1)./norm(true_depth))*100)

dd_err2=true_depth-inverted_depth;
dd_err2(end,end)=-3.1;
figure(8)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
ccmap2=makecolormap({'blue','white','red'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
%ccmap2=makecolormap({'white','goldenrod','limegreen','olivedrab','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,dd_err2,[-3.5:0.05:3.5],'edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Error (km)');

m_grid('box','fancy','grid','none','fontsize',14);
title('Error in Moho Depth')
fprintf('inverted depth error is %f\n',(norm(dd_err2)./norm(true_depth))*100)

%%
gg_err1=true_grav-pred_grav;
%gg_err2(end,end)=-3.1;
figure(9)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
ccmap2=makecolormap({'green','white','gold'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
%ccmap2=makecolormap({'white','goldenrod','limegreen','olivedrab','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,gg_err1,[-35:0.05:35],'edgecolor','none');
%ax=m_contfbar(.85,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Error (mGal)');

m_grid('box','fancy','grid','none','fontsize',14);
title('Error in Gravity Anomalies')
fprintf('predicted gravity error is %f\n',(norm(gg_err1)./norm(true_grav))*100)

gg_err2=true_grav-inverted_grav;
%gg_err2(end,end)=-3.1;

figure(10)
m_proj('miller','lon',[min(lon(:))+1.5 max(lon(:))-1.5],'lat',[min(lat(:))+1.5 max(lat(:))-1.5]);  
ccmap2=makecolormap({'green','white','gold'}, 128);
%ccmap2=makecolormap({'white','magenta','red','gold','aqua','darkblue'}, 128);
%ccmap2=makecolormap({'white','goldenrod','limegreen','olivedrab','orchid','royalblue','darkred'}, 128);
colormap(ccmap2);
[CS,CH]=m_contourf(lon,lat,gg_err2,[-35:0.05:35],'edgecolor','none');
%ax=m_contfbar(.85,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','levels','match','edgecolor','none');
ax=m_contfbar(.97,[.55 .95],CS,CH,'axfrac',.02,'endpiece','yes','edgecolor','none');
%ax=m_contfbar(.97,[.5 .9],[12.5 50.5],[15:0.05:50],'endpiece','yes','edgecolor','none');
ylabel(ax,'Error (mGal)');
m_grid('box','fancy','grid','none','fontsize',14);
title('Error in Gravity Anomalies')
fprintf('predicted gravity error is %f\n',(norm(gg_err2)./norm(true_grav))*100)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
lon1=min(lon(:)); lon2=max(lon(:));
lat1=min(lat(:)); lat2=max(lat(:));
lnn=linspace(lon1,lon2,100);
ltt=lat1+((lat1-lat2)/(lon1-lon2))*(lnn-lon1);
%loop for finding distance
latlon1=[ltt(1) lnn(1)];
for i=1:length(lnn)
    latlon2=[ltt(i) lnn(i)];
    [d1km(i), d2km(i)]=lldistkm(latlon1,latlon2);
end

moho_line1=griddata(lon,lat,true_depth,lnn,ltt);
moho_line2=griddata(lon,lat,pred_depth,lnn,ltt);
moho_line3=griddata(lon,lat,inverted_depth,lnn,ltt);
%%
figure(11)
plot(d2km(10:90)-d2km(10),moho_line1(10:90))
hold on
plot(d2km(10:90)-d2km(10),moho_line2(10:90))
plot(d2km(10:90)-d2km(10),moho_line3(10:90))
xlabel('Distance (km)')
ylabel('Moho Depth (km)')
legend('true moho','cGAN inverted moho','Bott''s inverted moho','location','best')

grav_line1=griddata(lon,lat,true_grav,lnn,ltt);
grav_line2=griddata(lon,lat,pred_grav,lnn,ltt);
grav_line3=griddata(lon,lat,inverted_grav,lnn,ltt);

figure(12)
plot(d2km(10:90)-d2km(10),grav_line1(10:90))
hold on
plot(d2km(10:90)-d2km(10),grav_line2(10:90))
plot(d2km(10:90)-d2km(10),grav_line3(10:90))
xlabel('Distance (km)')
ylabel('Gravity Anomalies (mGal)')
legend('true anomalies','cGAN inverted anomalies','Bott''s inverted anomalies','location','best')

figure(1)
m_line(lnn(10:90),ltt(10:90),'color','k')
figure(4)
m_line(lnn(10:90),ltt(10:90),'color','k')