function RunInversion(Xobs,Yobs,Zobs,Uxoss,Uyoss,Uzoss,ForwardModel,RangeParameters,AddPath,n,MaxGenerations,PopulationSize,ni)



% Subroutine to perform the inversion of the deformation data using a single-objective genetic algorithm






SourceType=ForwardModel(1);

U_obs=[Uxoss, Uyoss, Uzoss];  % row vector of the observed deformation
U_obs=U_obs'; % column vector of the observed deformation




fid = fopen([AddPath 'parameters.txt'],'wt'); 


for i=1:n


    Y=['Number of simulations run: ', num2str(i)];
    disp(Y);
    

    [X,FVAL]=RunGA(Xobs,Yobs,Zobs,Uxoss,Uyoss,Uzoss,ForwardModel,RangeParameters,MaxGenerations,PopulationSize,ni);
    

    P=[X FVAL]; 


    switch SourceType
        
        
        case 1 % sphere            
            Xc=X(1);
            Yc=X(2);
            Zc=X(3);
            DeltaV=1;      
            [Ux_calc,Uy_calc,Uz_calc]=deformation_sphere(Xobs,Yobs,Zobs,Xc,Yc,Zc,DeltaV,ni); % row vectors of the computed deformation
       
            U_calc=[Ux_calc, Uy_calc, Uz_calc]; % computed deformation
            U_calc=U_calc'; 
 
            DeltaV=inv(U_calc'*U_calc)*U_calc'*U_obs;  

            P=[X DeltaV FVAL];

            P=P';
            P=P(:);

            fprintf(fid,'%f %f %f %f %f \n',P);



        case 2 % rectangular prism 
            SourceRotation=ForwardModel(2);                       
            Xc=X(1);
            Yc=X(2);
            depth=X(3);
            l_a=X(4);
            l_b=X(5);
            l_c=X(6);
            if SourceRotation==1 % yes rotation
                alpha=X(7);
                delta=X(8);
                gamma=X(9);
            else
                alpha=0;
                delta=0;
                gamma=0;
            end
            eps_0=1;    
            [Ux_calc,Uy_calc,Uz_calc]=deformation_prism(Xobs,Yobs,Zobs,Xc,Yc,depth,l_a,l_b,l_c,alpha,delta,gamma,eps_0,ni); % row vectors of the computed deformation

            U_calc=[Ux_calc, Uy_calc, Uz_calc]; 
            U_calc=U_calc';        

            eps_0=inv(U_calc'*U_calc)*U_calc'*U_obs; 

            P=[X eps_0 FVAL];

            P=P';
            P=P(:);

            if SourceRotation==1
                fprintf(fid,'%f %f %f %f %f %f %f %f %f %f %f \n',P);
            else
                fprintf(fid,'%f %f %f %f %f %f %f %f \n',P);
            end



        case 3 % prolate or oblate ellipsoid 
            SourceRotation=ForwardModel(2);           
            Xc=X(1);
            Yc=X(2);
            depth=X(3);
            a=X(4);
            b=X(5);
            if SourceRotation==1 % yes rotation
                alpha=X(6);
                delta=X(7);
                gamma=X(8);
            else
                alpha=0;
                delta=0;
                gamma=0;
            end
            eps_0=1; 
            [Ux_calc,Uy_calc,Uz_calc]=deformation_prolate_oblate_ellipsoid(Xobs,Yobs,Zobs,Xc,Yc,depth,a,b,alpha,delta,gamma,eps_0,ni); % row vectors of the computed deformation
       
            U_calc=[Ux_calc, Uy_calc, Uz_calc]; 
            U_calc=U_calc';      
        
            eps_0=inv(U_calc'*U_calc)*U_calc'*U_obs;   

            P=[X eps_0 FVAL];

            P=P';
            P=P(:);

            if SourceRotation==1
                fprintf(fid,'%f %f %f %f %f %f %f %f %f %f \n',P);
            else
                fprintf(fid,'%f %f %f %f %f %f %f \n',P);
            end




        case 4 % triaxial ellipsoid
            SourceRotation=ForwardModel(2);
            Xc=X(1);
            Yc=X(2);
            depth=X(3);
            a=X(4);           
            b=X(5);
            c=X(6);
            if SourceRotation==1 % yes rotation
                alpha=X(7);
                delta=X(8);
                gamma=X(9);
            else
                alpha=0;
                delta=0;
                gamma=0;
            end
            eps_0=1;      
            [Ux_calc,Uy_calc,Uz_calc]=deformation_triaxial_ellipsoid(Xobs,Yobs,Zobs,Xc,Yc,depth,a,b,c,alpha,delta,gamma,eps_0,ni);  % row vectors of the computed deformation
      
            U_calc=[Ux_calc, Uy_calc, Uz_calc]; 
            U_calc=U_calc';        

            eps_0=inv(U_calc'*U_calc)*U_calc'*U_obs; %% minimi quadrati   

            P=[X eps_0 FVAL];

            P=P';
            P=P(:);

            if SourceRotation==1
                fprintf(fid,'%f %f %f %f %f %f %f %f %f %f %f \n',P);
            else
                fprintf(fid,'%f %f %f %f %f %f %f %f \n',P);
            end

        
        
               
        case 5 % cylinder
            Xc=X(1);
            Yc=X(2);
            rho_radial_obs=sqrt((Xc-Xobs).^2+(Yc-Yobs).^2); % radial distance between source and obs points
            depth=X(3);
            R=X(4);
            h=X(5);
            eps_0=1;
            [Ur_calc,Uz_calc]=deformation_cylinder(rho_radial_obs,Zobs,depth,R,h,eps_0,ni); 
            Ux_calc=Ur_calc.*(Xobs-Xc)./rho_radial_obs;   
            Uy_calc=Ur_calc.*(Yobs-Yc)./rho_radial_obs;

            U_calc=[Ux_calc, Uy_calc, Uz_calc];  
            U_calc=U_calc';         
        
            eps_0=inv(U_calc'*U_calc)*U_calc'*U_obs;

            P=[X eps_0 FVAL];

            P=P';
            P=P(:);
 
            fprintf(fid,'%f %f %f %f %f %f %f \n',P);


    end



end