function [minimum,position]=ResultsInversion(Xobs,Yobs,Zobs,Uxoss,Uyoss,Uzoss,P,AddPath,ForwardModel,RangeParameters,n,ni,str)



% Subroutine to analyze the results provided by the inversion algorithm. 
% The subroutine provides:
% (1) the plot of the horizontal and vertical displacements at the optimal parameters; 
% (2) the analysis of the 1D and 2D marginal distributions of the set of solutions provided by the inversion algorithm






SourceType=ForwardModel(1);


% observation points over a ground surface profile along the x-axis and with respect to the centre of the source
rho_radial_obs2=[0:100:5000]; 






%% Optimal solution 





switch SourceType
    
    
    
    case 1 % sphere               
        [minimum,position]=min(P(:,5));
        Xc_opt=P(position,1);
        Yc_opt=P(position,2);
        rho_radial_obs=sqrt((Xobs-Xc_opt).^2+(Yobs-Yc_opt).^2);   
        Zc_opt=P(position,3);
        DeltaV_opt=P(position,4);
        
        [Ux_opt,Uy_opt,Uz_opt]=deformation_sphere(Xobs,Yobs,Zobs,Xc_opt,Yc_opt,Zc_opt,DeltaV_opt,ni);             
        [Ur_opt,Uzr_opt]=deformation_sphere_radial(rho_radial_obs,Zc_opt,DeltaV_opt,ni);
        [Ur_opt2,Uzr_opt2]=deformation_sphere_radial(rho_radial_obs2,Zc_opt,DeltaV_opt,ni);           
       

        
    case 2 % rectangular prism 
        SourceRotation=ForwardModel(2);
        if SourceRotation==1 % yes rotation
            [minimum,position]=min(P(:,11));
        else
            [minimum,position]=min(P(:,8));
        end
        Xc_opt=P(position,1);
        Yc_opt=P(position,2);
        rho_radial_obs=sqrt((Xobs-Xc_opt).^2+(Yobs-Yc_opt).^2);   
        depth_opt=P(position,3);
        l_a_opt=P(position,4);
        l_b_opt=P(position,5);
        l_c_opt=P(position,6);
        if SourceRotation==1 % yes rotation
            alpha_opt=P(position,7);
            delta_opt=P(position,8);
            gamma_opt=P(position,9);
            eps_0_opt=P(position,10);
        else
            alpha_opt=0;
            delta_opt=0;
            gamma_opt=0;
            eps_0_opt=P(position,7);
        end

        [Ux_opt,Uy_opt,Uz_opt]=deformation_prism(Xobs,Yobs,Zobs,Xc_opt,Yc_opt,depth_opt,l_a_opt,l_b_opt,l_c_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni); 
        [Ur_opt,Uyr_opt,Uzr_opt]=deformation_prism_radial(rho_radial_obs,depth_opt,l_a_opt,l_b_opt,l_c_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni); 
        [Ur_opt2,Uyr_opt2,Uzr_opt2]=deformation_prism_radial(rho_radial_obs2,depth_opt,l_a_opt,l_b_opt,l_c_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni);           

        

    case 3 % prolate or oblate ellipsoid            
        SourceRotation=ForwardModel(2);
        if SourceRotation==1 % yes rotation
            [minimum,position]=min(P(:,10));
        else
            [minimum,position]=min(P(:,7));
        end
        Xc_opt=P(position,1);
        Yc_opt=P(position,2);
        rho_radial_obs=sqrt((Xobs-Xc_opt).^2+(Yobs-Yc_opt).^2);  
        depth_opt=P(position,3);
        a_opt=P(position,4);
        b_opt=P(position,5);
        if SourceRotation==1 % yes rotation
            alpha_opt=P(position,6);
            delta_opt=P(position,7);
            gamma_opt=P(position,8);
            eps_0_opt=P(position,9);
        else
            alpha_opt=0;
            delta_opt=0; 
            gamma_opt=0;
            eps_0_opt=P(position,6);
        end
           
        [Ux_opt,Uy_opt,Uz_opt]=deformation_prolate_oblate_ellipsoid(Xobs,Yobs,Zobs,Xc_opt,Yc_opt,depth_opt,a_opt,b_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni);             
        [Ur_opt,Uyr_opt,Uzr_opt]=deformation_prolate_oblate_ellipsoid_radial(rho_radial_obs,depth_opt,a_opt,b_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni); 
        [Ur_opt2,Uyr_opt2,Uzr_opt2]=deformation_prolate_oblate_ellipsoid_radial(rho_radial_obs2,depth_opt,a_opt,b_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni);



    case 4 % triaxial ellipsoid         
        SourceRotation=ForwardModel(2);
        if SourceRotation==1 % yes rotation
            [minimum,position]=min(P(:,11));
        else          
            [minimum,position]=min(P(:,8));
        end
        Xc_opt=P(position,1);
        Yc_opt=P(position,2);
        rho_radial_obs=sqrt((Xobs-Xc_opt).^2+(Yobs-Yc_opt).^2); 
        depth_opt=P(position,3);
        a_opt=P(position,4);
        b_opt=P(position,5);
        c_opt=P(position,6);
        if SourceRotation==1 % yes rotation
            alpha_opt=P(position,7);
            delta_opt=P(position,8);
            gamma_opt=P(position,9);
            eps_0_opt=P(position,10);
        else
            alpha_opt=0;
            delta_opt=0; 
            gamma_opt=0;
            eps_0_opt=P(position,7);
        end

        [Ux_opt,Uy_opt,Uz_opt]=deformation_triaxial_ellipsoid(Xobs,Yobs,Zobs,Xc_opt,Yc_opt,depth_opt,a_opt,b_opt,c_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni);            
        [Ur_opt,Uyr_opt,Uzr_opt]=deformation_triaxial_ellipsoid_radial(rho_radial_obs,depth_opt,a_opt,b_opt,c_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni);
        [Ur_opt2,Uyr_opt2,Uzr_opt2]=deformation_triaxial_ellipsoid_radial(rho_radial_obs2,depth_opt,a_opt,b_opt,c_opt,alpha_opt,delta_opt,gamma_opt,eps_0_opt,ni);                    

        
        
    case 5 % cylinder                
        [minimum,position]=min(P(:,7));     
        Xc_opt=P(position,1);
        Yc_opt=P(position,2);
        rho_radial_obs=sqrt((Xobs-Xc_opt).^2+(Yobs-Yc_opt).^2);  
        depth_opt=P(position,3);
        R_opt=P(position,4);
        h_opt=P(position,5);
        eps_0_opt=P(position,6); 
                        
        [Ur_opt,Uzr_opt]=deformation_cylinder_2(rho_radial_obs,Zobs,depth_opt,R_opt,h_opt,eps_0_opt,ni); 
        Uz_opt=Uzr_opt;
        Ux_opt=Ur_opt.*(Xobs-Xc_opt)./rho_radial_obs; 
        Uy_opt=Ur_opt.*(Yobs-Yc_opt)./rho_radial_obs;
           
        Zobs2=zeros(1,length(rho_radial_obs2));
        [Ur_opt2,Uzr_opt2]=deformation_cylinder_2(rho_radial_obs2,Zobs2,depth_opt,R_opt,h_opt,eps_0_opt,ni);    

                              
end 
  




figure(1) 

set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% QUIVER PLOT ON THE OPTIMAL HORIZONTAL SOLUTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot('Position',[0.1 0.15 0.35 0.7])

scale=5*10^4;

% conversion from UTM to lat/lon
mstruct = defaultm('utm');
mstruct.zone= '33S';
mstruct.geoid = almanac('earth','wgs84','meters');
mstruct = defaultm(mstruct);

[Yc_opt_lat,Xc_opt_lon] = projinv(mstruct,Xc_opt,Yc_opt);
[Yobs_lat,Xobs_lon] = projinv(mstruct,Xobs,Yobs);

geoplot(Yc_opt_lat,Xc_opt_lon,'Marker','pentagram','MarkerSize',14,'MarkerFaceColor','c','MarkerEdgeColor','k'); % plot of the position of the source

geobasemap satellite

hold on 
geoscatter(Yobs_lat,Xobs_lon,[18],'filled','MarkerFaceColor','black'); % plot of the observation points

for i=1:length(Xobs)
    begin_point=[Yobs_lat(i),Xobs_lon(i)];
    [end_point_lat,end_point_lon] = projinv(mstruct,Xobs(i)+Uxoss(i)*scale,Yobs(i)+Uyoss(i)*scale);
    end_point=[end_point_lat,end_point_lon];
    [end_point_lat2,end_point_lon2] = projinv(mstruct,real(Xobs(i)+Ux_opt(i)*scale),real(Yobs(i)+Uy_opt(i)*scale));
    end_point2=[real(end_point_lat2),real(end_point_lon2)]; % the imaginary part cancel to zero
    
    hold on 
    plot_arrow_geoplot(begin_point,end_point,'color',[0.6350 0.0780 0.1840],'linewidth',2.0) % observed deformation
    hold on
    plot_arrow_geoplot(begin_point,end_point2,'color','c','linewidth',2.0) % computed deformation

    text(begin_point(1),begin_point(2),str(i),'Color','w','FontWeight','bold','FontSize',14)
end

% % % % % % % % % % % %  unit of measurement of the displacement -- you can change its position depending on your needs
begin_x=Xobs(end)+2550;
begin_y=Yobs(end)-900;
[begin_point_lat3,begin_point_lon3] = projinv(mstruct,begin_x,begin_y);
[end_point_lat3,end_point_lon3] = projinv(mstruct,begin_x+0.01*scale,begin_y); % unit of measurement: 1 cm
end_point3=[end_point_lat3,end_point_lon3];
begin_point3=[begin_point_lat3,begin_point_lon3];
plot_arrow_geoplot(begin_point3,end_point3,'color','w','linewidth',1.5)

text(begin_point3(1),begin_point3(2),'1 cm','Color','white','FontSize',14)
% % % % % % % % % % % %

set(gca,'fontsize',14);
lgd=legend('','','','','Observed displacement','Computed displacement','','','','','','','','','','','','','','','','','','','color','black','textcolor','white','interpreter','latex');
lgd.FontSize=18;
set(gca,'TickDir','out');

title('Horizontal displacement')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RADIAL AND VERTICAL DISPLACEMENTS ON THE OPTIMAL SOLUTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot('Position',[0.55 0.6 0.35 0.25])
plot(rho_radial_obs2/10^3,Ur_opt2,'-','LineWidth',1.0,'Color','k') % analytical solution (forward model)
hold on
plot(rho_radial_obs/10^3,Ur_opt,'s','LineWidth',0.5,'MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',12)  % computed displacement
text(rho_radial_obs/10^3,Ur_opt,str,'FontWeight','bold','color','k','FontSize',14)
hold on 
plot(0,0,'Marker','pentagram','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',14)
set(gca,'TickDir','out');
set(gca,'fontsize',14);
ylabel('$\textbf{U}_r$ [m]','interpreter','latex','FontSize',20);
xlabel('$\rho$ [km]','Interpreter','latex','FontSize',20)
title('Radial displacement')
if SourceType==1||SourceType==2||SourceType==3
    lgd=legend('Analytical solution','Computed $\textbf{U}_r$','','interpreter','latex');
else 
    lgd=legend('Semi-Analytical solution','Computed $\textbf{U}_r$','','interpreter','latex');    
end
lgd.FontSize=18;
set(gca,'Color',[.8 .8 .8])


subplot('Position',[0.55 0.15 0.35 0.25])
plot(rho_radial_obs2/10^3,Uzr_opt2,'-','LineWidth',1.0,'Color','k') % analytical solution (forward model)
hold on
plot(rho_radial_obs/10^3,Uzr_opt,'s','LineWidth',0.5,'MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',12)  % computed displacement
hold on 
plot(rho_radial_obs/10^3,Uzoss,'s','LineWidth',0.5,'MarkerEdgeColor','k','MarkerFaceColor',[0.6350 0.0780 0.1840],'MarkerSize',12) % observed deformation
hold on 
plot(0,0,'Marker','pentagram','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',14)
set(gca,'TickDir','out');
set(gca,'fontsize',14);
ylabel('$\textbf{U}_z$ [m]','Interpreter','latex','FontSize',20);
xlabel('$\rho$ [km]','Interpreter','latex','FontSize',20)
lgd=legend('','Computed $\textbf{U}_z$','Observed $\textbf{U}_z$','interpreter','latex');
lgd.FontSize=18;
title('Vertical displacement')
set(gca,'Color',[.8 .8 .8])

saveas(gcf,[AddPath 'Optimal_Displacement'],'fig');




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MARGINAL DISTRIBUTIONS PLOT 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if n>1
    
    nbins=50;

    NumberParameters=size(P,2)-1; % number of parameters

    % we insert the range of values of the least squares parameter
    RangeParameters2=[RangeParameters(1,:) min(P(:,NumberParameters)); RangeParameters(2,:) max(P(:,NumberParameters))];

    % we restrict the range of each parameter from the minimum value to the maximum value in the set of GA solutions
    for i=1:NumberParameters
        RangeParameters2(1,i)=min(P(:,i));
        RangeParameters2(2,i)=max(P(:,i));
    end


    switch SourceType 
        
        case 1 % sphere              
            label=char('$x_c$ [km]','$y_c$ [km]','$z_c$ [km]','$\Delta V$ [km$^3$]','interpreter','latex');
            label=cellstr(label);
         
        case 2 % rectangular prism         
            SourceRotation=ForwardModel(2);
            if SourceRotation==1 % yes rotation
                label=char('$x_c$ [km]','$y_c$ [km]','$d$ [km]','$l_a$ [km]','$l_b$ [km]','$l_c$ [km]','$\alpha$ [$^{\circ}$]','$\delta$ [$^{\circ}$]','$\gamma$ [$^{\circ}$]','$\varepsilon_0$');
            else
                label=char('$x_c$ [km]','$y_c$ [km]','$d$ [km]','$l_a$ [km]','$l_b$ [km]','$l_c$ [km]','$\varepsilon_0$');            
            end
            label=cellstr(label);  
 
        case 3 % prolate or oblate ellipsoid       
            SourceRotation=ForwardModel(2);
            if SourceRotation==1
                label=char('$x_c$ [km]','$y_c$ [km]','$d$ [km]','$a$ [km]','$b$ [km]','$\alpha$ [$^{\circ}$]','$\delta$ [$^{\circ}$]','$\gamma$ [$^{\circ}$]','$\varepsilon_0$');
            else
                label=char('$x_c$ [km]','$y_c$ [km]','$d$ [km]','$a$ [km]','$b$ [km]','$\varepsilon_0$');            
            end
            label=cellstr(label);

       
        case 4 % triaxial ellipsoid         
            SourceRotation=ForwardModel(2);
            if SourceRotation==1
                label=char('$x_c$ [km]','$y_c$ [km]','$d$ [km]','$a$ [km]','$b$ [km]','$c$ [km]','$\alpha$ [$^{\circ}$]','$\delta$ [$^{\circ}$]','$\gamma$ [$^{\circ}$]','$\varepsilon_0$');
            else
                label=char('$x_c$ [km]','$y_c$ [km]','$d$ [km]','$a$ [km]','$b$ [km]','$c$ [km]','$\varepsilon_0$');            
            end
            label=cellstr(label);

        case 5 % cylinder        
            label=char('$x_c$ [km]','$y_c$ [km]','$d$ [km]','$R$ [km]', '$h$ [km]','$\varepsilon_0$');
            label=cellstr(label);

    end


    figure(2)

    set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);

    DrawMarginalDistributions(P,RangeParameters2,nbins,NumberParameters,label,position,ForwardModel)

    saveas(gcf,[AddPath 'Marginal_Distributions'],'fig');
    
end





%%%%%%%%%%%%%%%%%%%%%%
%% OPTIMAL SOURCE PLOT 
%%%%%%%%%%%%%%%%%%%%%%

if SourceType~=1

    figure(3)

    DrawOptimalSourceInversion(P,position,ForwardModel)

    saveas(gcf,[AddPath 'Optimal_Source'],'fig');

end