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

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

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

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

filename_lat = fullfile('.', 'tested_data', 'YY_predicted1.txt');
lat=importdata(filename_lat);

figure(1)
m_proj('miller','lon',[min(lon(:)) max(lon(:))],'lat',[min(lat(:)) max(lat(:))]);   % Projection
%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(.85,[.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 disturbance (mGal)');
m_grid('box','fancy','grid','none','fontsize',14);

%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(2)
m_proj('miller','lon',[min(lon(:)) max(lon(:))],'lat',[min(lat(:)) max(lat(:))]);   % Projection
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(.85,[.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(3)
m_proj('miller','lon',[min(lon(:)) max(lon(:))],'lat',[min(lat(:)) max(lat(:))]);   % Projection
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(.85,[.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 Moho Depth')

figure(4)
m_proj('miller','lon',[min(lon(:)) max(lon(:))],'lat',[min(lat(:)) max(lat(:))]);   % Projection
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(.85,[.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('inverted Moho Depth')

dd_err1=true_depth-pred_depth;
figure(5)
m_proj('miller','lon',[min(lon(:)) max(lon(:))],'lat',[min(lat(:)) max(lat(:))]);   % Projection
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,[-2:0.05:2],'edgecolor','none');
ax=m_contfbar(.85,[.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('predicted depth error is %f\n',(norm(dd_err1)./norm(true_depth))*100)

dd_err2=true_depth-inverted_depth;
figure(6)
m_proj('miller','lon',[min(lon(:)) max(lon(:))],'lat',[min(lat(:)) max(lat(:))]);   % Projection
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,[-2:0.05:2],'edgecolor','none');
ax=m_contfbar(.85,[.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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
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);

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(7)
plot(moho_line1(10:90))
hold on
plot(moho_line2(10:90))
plot(moho_line3(10:90))
legend('true moho','predicted moho','inverted moho','location','best')