function DrawOptimalSourceInversion(P,position,ForwardModel)

SourceType=ForwardModel(1);

switch SourceType


    case 2 % rectangular prism

        SourceRotation=ForwardModel(2);
     
        Xc=0;
        Yc=0;
        Zc=0;
        % Optimal parameters
        depth=P(position,3)/10^3; % [km]
        l_a=P(position,4)/10^3;
        l_b=P(position,5)/10^3;
        l_c=P(position,6)/10^3;
        if SourceRotation==1 % yes rotation
            alpha=P(position,7);
            delta=P(position,8);
            gamma=P(position,9);
        else
            alpha=0;
            delta=0;
            gamma=0;
        end
        

        z=Zc+l_c/2;
        height=l_c;

        z2=[z;z;z;z;z]; % upper base
        z1=z2-height; % lower base

        x=[Xc-l_a/2, Xc+l_a/2, Xc+l_a/2, Xc-l_a/2, Xc-l_a/2];
        y=[Yc-l_b/2, Yc-l_b/2, Yc+l_b/2, Yc+l_b/2, Yc-l_b/2];

        h=surf([x;x].', [y;y].', [z2,z1],'FaceColor',[0.8500 0.3250 0.0980],'FaceAlpha',.3,'EdgeColor',[0.8500 0.3250 0.0980]);

        rotate(h, [0 0 1], alpha);
        rotate(h, [0 1 0], delta);
        rotate(h, [1 0 0], gamma);


        [R] = rotation(alpha,delta,gamma);    
        P0 = [0 0 0]'; 
        P1 = P0+l_a*R(:,1)/2+l_c*R(:,3)/2-l_b/2*R(:,2);
        P2 = P1-l_c*R(:,3);
        P3 = P2-l_a*R(:,1);
        P4 = P1-l_a*R(:,1);
        Q1 = P0+l_a*R(:,1)/2+l_c*R(:,3)/2+l_b/2*R(:,2);
        Q2 = Q1-l_c*R(:,3);
        Q3 = Q2-l_a*R(:,1);
        Q4 = Q1-l_a*R(:,1);
        Zc=max([P1(3),P2(3),P3(3),P4(3),Q1(3),Q2(3),Q3(3),Q4(3)]);
        Zc=-Zc+depth;

        a_point=P0+l_a*R(:,1)/2;
        b_point=P0-l_b/2*R(:,2);
        c_point=P0+l_c*R(:,3)/2;



        Xc=P(position,1)/10^3;
        Yc=P(position,2)/10^3;

       
        X = get(h, 'XData');
        Y = get(h, 'YData');
        Z = get(h, 'ZData');


        surf(X+Xc,Y+Yc,Z+Zc,'FaceColor',[0.8500 0.3250 0.0980],'FaceAlpha',.3,'EdgeColor',[0.8500 0.3250 0.0980])

        
        hold on
        plot3([Xc,a_point(1)+Xc],[Yc,a_point(2)+Yc],[Zc,a_point(3)+Zc],'k-','LineWidth',1.5)
        hold on
        plot3([Xc,b_point(1)+Xc],[Yc,b_point(2)+Yc],[Zc,b_point(3)+Zc],'r-','LineWidth',1.5)
        hold on
        plot3([Xc,c_point(1)+Xc],[Yc,c_point(2)+Yc],[Zc,c_point(3)+Zc],'b-','LineWidth',1.5)
        hold on
        plot3([Xc,Xc],[Yc,Yc],[0,depth],'linestyle','--','LineWidth',0.5,'Color',[.7 .7 .7])      
        
              
        hold on
        plot3(Xc,Yc,Zc,'Marker','pentagram','MarkerFaceColor','c','MarkerEdgeColor','k','MarkerSize',10) % center of the source 


        % major base
        hold on
        p1=fill3( [Xc+P1(1) Xc+P4(1) Xc+Q4(1) Xc+Q1(1)], [Yc+P1(2) Yc+P4(2) Yc+Q4(2) Yc+Q1(2)], [Zc+P1(3) Zc+P4(3) Zc+Q4(3) Zc+Q1(3)],[0.8500 0.3250 0.0980]);
        p1.FaceAlpha=0.3;
        p1.EdgeColor=[0.8500 0.3250 0.0980];


 
        % minor base
        hold on
        p2=fill3( [Xc+P2(1) Xc+P3(1) Xc+Q3(1) Xc+Q2(1)], [Yc+P2(2) Yc+P3(2) Yc+Q3(2) Yc+Q2(2)], [Zc+P2(3) Zc+P3(3) Zc+Q3(3) Zc+Q2(3)],[0.8500 0.3250 0.0980]);
        p2.FaceAlpha=0.3;
        p2.EdgeColor=[0.8500 0.3250 0.0980];


        set(gca,'TickDir','out');
        set(gca,'Fontsize',8)
        axis equal

        grid on
        xlabel('$x$ [km]','Interpreter','latex','FontSize',10)
        ylabel('$y$ [km]','Interpreter','latex','FontSize',10)
        zlabel('$z$ [km]','Interpreter','latex','FontSize',10)

        lgd=legend('','$l_a/2$','$l_b/2$','$l_c/2$','$d$','Interpreter','latex');
        % legend size
        lgd.FontSize=10;






    case 3 % prolate or oblate ellipsoid

        SourceRotation=ForwardModel(2);

        % Optimal parameters
        depth=P(position,3)/10^3; % [km]
        a=P(position,4)/10^3;
        b=P(position,5)/10^3;
        c=b;
        if SourceRotation==1 % yes rotation
            alpha=P(position,6);
            delta=P(position,7);
            gamma=P(position,8);
        else
            alpha=0;
            delta=0; 
            gamma=0;
        end
        alpha2=0;        
        delta2=180; 
        gamma2=0;


        [x,y,z] = ellipsoid(0,0,0,a,b,c);
        h=surf(x,y,z,[0.8500 0.3250 0.0980],'FaceAlpha',.3,'EdgeColor',[0.8500 0.3250 0.0980]);
 

        rotate(h, [0 0 1], -alpha);
        rotate(h, [0 1 0], -delta);
        rotate(h, [1 0 0], gamma);

        % rotate(h, [0 0 1], -alpha2);
        % rotate(h, [0 1 0], -delta2);
        % rotate(h, [1 0 0], gamma2);
        
        [R] = rotation(-alpha,-delta,gamma);                 
        
        % [R1] = rotation(-alpha,-delta,gamma); 
        % [R2] = rotation(-alpha2,-delta2,gamma2);                  
        % R=R2*R1;

        P0 = [0 0 0]';
        P1 = P0+b*R(:,2);
        P2 = P0-b*R(:,2);    
        Q1 = P0+c*R(:,3);
        Q2 = P0-c*R(:,3);
        R1 = P0+a*R(:,1);
        R2 = P0-a*R(:,1);   
        Zc=max([P1(3),P2(3),Q1(3),Q2(3),R1(3),R2(3)]);
        Zc=-Zc+depth;

        Xc=P(position,1)/10^3;
        Yc=P(position,2)/10^3;

       
        X = get(h, 'XData');
        Y = get(h, 'YData');
        Z = get(h, 'ZData');


        surf(X+Xc,Y+Yc,Z+Zc,'FaceColor',[0.8500 0.3250 0.0980],'FaceAlpha',.3,'EdgeColor',[0.8500 0.3250 0.0980]);


        hold on
        plot3([Xc,R1(1)+Xc],[Yc,R1(2)+Yc],[Zc,R1(3)+Zc],'k-','LineWidth',1.5)
        hold on
        plot3([Xc,P1(1)+Xc],[Yc,P1(2)+Yc],[Zc,P1(3)+Zc],'r-','LineWidth',1.5)
        hold on
        plot3([Xc,Q1(1)+Xc],[Yc,Q1(2)+Yc],[Zc,Q1(3)+Zc],'b-','LineWidth',1.5)
        hold on
        plot3([Xc,Xc],[Yc,Yc],[0,depth],'linestyle','--','LineWidth',0.5,'Color',[.7 .7 .7])



        hold on
        plot3(Xc,Yc,Zc,'Marker','pentagram','MarkerFaceColor','c','MarkerEdgeColor','k','MarkerSize',10) % center of the source 

        set(gca,'TickDir','out');
        set(gca,'Fontsize',8)
        axis equal

        grid on
        xlabel('$x$ [km]','Interpreter','latex','FontSize',10)
        ylabel('$y$ [km]','Interpreter','latex','FontSize',10)
        zlabel('$z$ [km]','Interpreter','latex','FontSize',10)

        lgd=legend('','$a$','$b$','$c$','$d$','Interpreter','latex');
        % legend size
        lgd.FontSize=10;




       






    case 4 % triaxial ellipsoid

        SourceRotation=ForwardModel(2);

        % Optimal parameters
        depth=P(position,3)/10^3; % [km]
        a=P(position,4)/10^3;
        b=P(position,5)/10^3;
        c=P(position,6)/10^3;
        if SourceRotation==1 % yes rotation
            alpha=P(position,7);
            delta=P(position,8);
            gamma=P(position,9);
        else
            alpha=0;
            delta=0; 
            gamma=0;
        end


        [x,y,z] = ellipsoid(0,0,0,a,b,c);
        h=surf(x,y,z,[0.8500 0.3250 0.0980],'FaceAlpha',.3,'EdgeColor',[0.8500 0.3250 0.0980]);


        rotate(h, [0 0 1], -alpha);
        rotate(h, [0 1 0], -delta);
        rotate(h, [1 0 0], gamma);


        [R] = rotation(-alpha,-delta,gamma);    
        P0 = [0 0 0]';
        P1 = P0+b*R(:,2);
        P2 = P0-b*R(:,2);    
        Q1 = P0+c*R(:,3);
        Q2 = P0-c*R(:,3);
        R1 = P0+a*R(:,1);
        R2 = P0-a*R(:,1);   
        Zc=max([P1(3),P2(3),Q1(3),Q2(3),R1(3),R2(3)]);
        Zc=-Zc+depth;

        Xc=P(position,1)/10^3;
        Yc=P(position,2)/10^3;

       
        X = get(h, 'XData');
        Y = get(h, 'YData');
        Z = get(h, 'ZData');


        surf(X+Xc,Y+Yc,Z+Zc,'FaceColor',[0.8500 0.3250 0.0980],'FaceAlpha',.3,'EdgeColor',[0.8500 0.3250 0.0980])


        hold on
        plot3([Xc,R1(1)+Xc],[Yc,R1(2)+Yc],[Zc,R1(3)+Zc],'k-','LineWidth',1.5)
        hold on
        plot3([Xc,P1(1)+Xc],[Yc,P1(2)+Yc],[Zc,P1(3)+Zc],'r-','LineWidth',1.5)
        hold on
        plot3([Xc,Q1(1)+Xc],[Yc,Q1(2)+Yc],[Zc,Q1(3)+Zc],'b-','LineWidth',1.5)
        hold on
        plot3([Xc,Xc],[Yc,Yc],[0,depth],'linestyle','--','LineWidth',0.5,'Color',[.7 .7 .7])



        hold on
        plot3(Xc,Yc,Zc,'Marker','pentagram','MarkerFaceColor','c','MarkerEdgeColor','k','MarkerSize',10) % center of the source 

        set(gca,'TickDir','out');
        set(gca,'Fontsize',8)
        axis equal

        grid on
        xlabel('$x$ [km]','Interpreter','latex','FontSize',10)
        ylabel('$y$ [km]','Interpreter','latex','FontSize',10)
        zlabel('$z$ [km]','Interpreter','latex','FontSize',10)

        lgd=legend('','$a$','$b$','$c$','$d$','Interpreter','latex');
        % legend size
        lgd.FontSize=10;



    case 5 % cylinder

        % Optimal parameters
        Xc=P(position,1)/10^3;
        Yc=P(position,2)/10^3;
        depth=P(position,3)/10^3;
        R=P(position,4)/10^3;
        h=P(position,5)/10^3;


        [X,Y,Z]=cylinder(R,100);  
        Z=Z*h;
        surf(X+Xc,Y+Yc,-Z+depth,'FaceColor',[0.8500 0.3250 0.0980],'FaceAlpha',.3,'EdgeColor','none');
       
        Zc=depth-h/2;

        hold on
        plot3(Xc,Yc,Zc,'Marker','pentagram','MarkerFaceColor','c','MarkerEdgeColor','k','MarkerSize',10) % center of the source 

        hold on
        plot3([Xc,Xc+R],[Yc,Yc],[depth,depth],'k-','LineWidth',1.5)
        hold on 
        plot3([Xc,Xc],[Yc,Yc],[depth,depth-h],'r-','LineWidth',1.5)
        hold on
        plot3([Xc,Xc],[Yc,Yc],[0,depth],'linestyle','--','LineWidth',0.5,'Color',[.7 .7 .7])

        
        set(gca,'TickDir','out');
        set(gca,'Fontsize',8)
        axis equal

        grid on
        xlabel('$x$ [km]','Interpreter','latex','FontSize',10)
        ylabel('$y$ [km]','Interpreter','latex','FontSize',10)
        zlabel('$z$ [km]','Interpreter','latex','FontSize',10)

        lgd=legend('','','$R$','$h$','$d$','Interpreter','latex');
        % legend size
        lgd.FontSize=10;







end