function []=mainForward

% Subroutine 
% (1) to compute the horizontal and vertical displacements,
% (2) for the visualization of the results. 
% You can choose between different forward models that describe the deformations induced 
% by thermo-pore-elastic sources of spherical, ellipsoidal, cylindrical and rectangular prism shapes 


%% INPUT PARAMETERS OF THE SOFTWARE
[AddPath,ForwardModel,ValueParameters,rho_radial_obs2,ni]=InputParametersForward;


SourceType=ForwardModel(1);



%% DISPLAY
disp(' ')
disp(' ')
disp('           FORWARD ANALYSIS FOR MODELING THE THERMO-PORO-ELASTIC GROUND DISPLACEMENTS           ')
disp(' ')
disp(' ')



% Computed displacement

switch SourceType
        
    case 1 % sphere

        Xc=ValueParameters(1);
        Yc=ValueParameters(2);
        Zc=ValueParameters(3);
        DeltaV=ValueParameters(4);

        [Ur2,Uzr2]=deformation_sphere_radial(rho_radial_obs2,Zc,DeltaV,ni); 


    case 2 % rectangular prism 

        Xc=ValueParameters(1);
        Yc=ValueParameters(2);
        depth=ValueParameters(3);
        l_a=ValueParameters(4);
        l_b=ValueParameters(5);
        l_c=ValueParameters(6); 
        eps_0=ValueParameters(7);

        SourceRotation=ForwardModel(2);
        if SourceRotation==1
            alpha=ValueParameters(8);
            delta=ValueParameters(9);
            gamma=ValueParameters(10);
        else
            alpha=0;
            delta=0;
            gamma=0;            
        end

        [Ux,Uy,Uz]=deformation_prism_radial(rho_radial_obs2,depth,l_a,l_b,l_c,alpha,delta,gamma,eps_0,ni);
						
    


    case 3 % prolate or oblate ellipsoid

        Xc=ValueParameters(1);
        Yc=ValueParameters(2);
        depth=ValueParameters(3);
        a=ValueParameters(4);
        b=ValueParameters(5);
        eps_0=ValueParameters(6);

        SourceRotation=ForwardModel(2);
        if SourceRotation==1
            alpha=ValueParameters(7);
            delta=ValueParameters(8);
            gamma=ValueParameters(9);
        else
            alpha=0;
            delta=0;
            gamma=0;            
        end
        
        [Ux,Uy,Uz]=deformation_prolate_oblate_ellipsoid_radial(rho_radial_obs2,depth,a,b,alpha,delta,gamma,eps_0,ni);

    


    case 4 % triaxial ellipsoid

        Xc=ValueParameters(1);
        Yc=ValueParameters(2);
        depth=ValueParameters(3);
        a=ValueParameters(4);
        b=ValueParameters(5);
        c=ValueParameters(6); 
        eps_0=ValueParameters(7);

        SourceRotation=ForwardModel(2);
        if SourceRotation==1
            alpha=ValueParameters(8);
            delta=ValueParameters(9);
            gamma=ValueParameters(10);
        else
            alpha=0;
            delta=0;
            gamma=0;            
        end

        [Ux,Uy,Uz]=deformation_triaxial_ellipsoid_radial(rho_radial_obs2,depth,a,b,c,alpha,delta,gamma,eps_0,ni); 

    
    case 5 % cylinder 
       
        Xc=ValueParameters(1);
        Yc=ValueParameters(2);
        depth=ValueParameters(3);
        R=ValueParameters(4);
        h=ValueParameters(5);
        eps_0=ValueParameters(6); 
           
        Zobs2=zeros(1,length(rho_radial_obs2));
        [Ux,Uz]=deformation_cylinder_2(rho_radial_obs2,Zobs2,depth,R,h,eps_0,ni); 


end





figure(1) 

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


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

% subplot('Position',[0.55 0.6 0.35 0.25])
subplot(1,2,1)
plot(rho_radial_obs2/10^3,Ux,'-','LineWidth',1.0,'Color','k') % analytical solution (forward model)
hold on
plot(rho_radial_obs2/10^3,Ux,'Marker','square','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',12)
hold on
plot(0,0,'Marker','pentagram','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',14)
set(gca,'TickDir','out');
set(gca,'fontsize',18);
ylabel('$\textbf{U}_x$ [m]','interpreter','latex','FontSize',24);
xlabel('$\rho$ [km]','Interpreter','latex','FontSize',24)
title('Horizontal displacement','FontSize',24)
if SourceType==1||SourceType==2||SourceType==3
    lgd=legend('Analytical solution','','','interpreter','latex');
else 
    lgd=legend('Semi-Analytical solution','','','interpreter','latex');    
end
lgd.FontSize=18;
set(gca,'Color',[.8 .8 .8])


% subplot('Position',[0.55 0.15 0.35 0.25])
subplot(1,2,2)
plot(rho_radial_obs2/10^3,Uz,'-','LineWidth',1.0,'Color','k') % analytical solution (forward model)
hold on
plot(rho_radial_obs2/10^3,Uz,'Marker','square','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',12)
hold on 
plot(0,0,'Marker','pentagram','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',14)
set(gca,'TickDir','out');
set(gca,'fontsize',18);
ylabel('$\textbf{U}_z$ [m]','Interpreter','latex','FontSize',24);
xlabel('$\rho$ [km]','Interpreter','latex','FontSize',24)
title('Vertical displacement','FontSize',24)
set(gca,'Color',[.8 .8 .8])

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




%%%%%%%%%%%%%%
%% SOURCE PLOT 
%%%%%%%%%%%%%%

if SourceType~=1

    figure(2)

    DrawOptimalSourceForward(ValueParameters,ForwardModel)

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

end