Code covered by the BSD License  

Highlights from
Ship Stability and Hydrostatic Calculation

image thumbnail
from Ship Stability and Hydrostatic Calculation by Baharuddin Ali
Plot body plan, hydrostatic, GZ calculations

HidroMatik(varargin)
function varargout = HidroMatik(varargin)
% HIDROMATIK for calculation of hydrostatic and GZ curve for ship stability
%
%      Input data of ship can be input using Microsoft Excel
%      Line 1 = Main dimensions, [ LOA  LPP  B  H  T ]
%      Line 2 - end : [ z   y   x ]
%      Coloum 1 = z , station number
%      Coloum 2 = y , ordinate of bodyplan
%      Coloum 3 = x , abscess of bodyplan
%      
%      This is may a simple program, you can develop for detail.
%
%      Baharuddin Ali
%      Jl. Rambutan 106, Cerme, Grogol, Kediri, 64151 INDONESIA
%      + 62-354-777290, b4li313@yahoo.com

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @HidroMatik_OpeningFcn, ...
                   'gui_OutputFcn',  @HidroMatik_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before HidroMatik is made visible.
function HidroMatik_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to HidroMatik (see VARARGIN)

% Choose default command line output for HidroMatik
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes HidroMatik wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = HidroMatik_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in pbBodyplan.
function pbBodyplan_Callback(hObject, eventdata, handles)
% hObject    handle to pbBodyplan (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --- Executes on button press in pb3D.
function pb3D_Callback(hObject, eventdata, handles)
% hObject    handle to pb3D (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --- Executes on button press in pbHidros.
function pbHidros_Callback(hObject, eventdata, handles)
% hObject    handle to pbHidros (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --- Executes on button press in pbProses.
function pbProses_Callback(hObject, eventdata, handles)
% hObject    handle to pbProses (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global A UK KG Phideg GZ pilih

ttk   = A(2:end,1:3);
k     = diff(ttk(:,1));
k2    = find(k ~= 0);
k2(length(k2) +1)   = length(ttk);
JumST = length(k2);
z = ttk(:,1); y = ttk(:,2); x = ttk(:,3);
h = UK(2)/10; 
d1 = linspace(0.15*UK(5),UK(5),20);

popup_sel_index = get(handles.popProses, 'Value');
switch popup_sel_index
    case 1
        pilih = 1;
    case 2
        pilih = 2;
    case 3
        pilih = 3; 
end

if pilih == 3
    KG = str2double(get(handles.eKG,'string')); % isilah KG yang diinginkan
    if KG <= 0
        errordlg('KG SHOULD BE > 0','Warning');
        return
    end
    StabGZ;
end

for k=1:length(d1)
    d=d1(k);j=0;
    for i = 1: JumST
        x1 = x(find(z == z(k2(i)))); % mencari data tiap station
        y1 = y(find(z == z(k2(i))));
        z1 = z(find(z == z(k2(i))));
        %------------------------------ Gambar body plan -------
        if pilih == 1
            figure(1);
            if (z(k2(i)) <=5)
                plot(-x1,y1); hold on
            else
                plot(x1,y1); hold on
            end
        end
        %-------------------------------------------------------
        y2 = y1(find(y1 <= d));
        x2 = x1(find(y1 <= d));
        z2 = z1(find(y1 <= d));
    
        if length(y2) > 1  % hanya titik yang di bawah sarat
            if y2(end) == d
                x2 = x2; y2 = y2; z2= z2;
            else
                y1f=diff(y1);
                y1p=find(y1f==0);
                for p = 1:length(y1p)
                    y1(y1p(p)+1)=y1(y1p(p))+0.001;
                end
            
                xn = interp1(y1,x1,d);      
                x2(end+1) = xn;            
                y2(end+1) = d;
                z2(end+1) = z2(1);
            end
        
            if pilih == 1
                if (z2(1) <= 5)
                    plot(-x2,y2,'r'); hold on
                else
                    plot(x2,y2,'r'); hold on
                end
            end
        
            z3=z2;y3=y2;x3=x2;
        else
            z3=z1(i);y3=0;x3=0;
        end
        if pilih == 2 & length(y3) > 1 % hanya station yang punya dua titik atau lebih
            j=j+1;
            x31 = -sort(-x3);
            y31 = -sort(-y3);
            X   = [x31;-x3(2:end)];
            Y   = [y31;y3(2:end)];
     
            xa  = [-x3(end) x3(end) X'];
            ya  = [d d Y'];
            Area1   = TtkBerat( xa,ya );
            LuasST(j)  = Area1(1);
            KbST(j)    = Area1(3);
            PanjGirthST(j) = Area1(4)-( 2*x3(end) );
            NoST(j)    = z3(1);
            Xb(j)      = x3(end);
       
            %plot(xa,ya,Area1(2),Area1(3),'or');
        end
    end
    
    if pilih == 1
        axis equal
        plot([-max(x)*1.2,max(x)*1.2],[UK(5),UK(5)],'--r');
    elseif pilih == 2
        NXLP = [NoST',Xb',LuasST',PanjGirthST',KbST'];
        Bd      = 2*NXLP( find(NXLP(:,1)==5),2);          % Lebar midship 
        Amd     = NXLP( find(NXLP(:,1)==5),3);            % Luasan midship
        Vold    = trapz(NoST*h,LuasST);                   % Volume carene


        Displd  = 1.025 * Vold;                          % Displacement 
        Awld    = 2 * trapz(NoST*h,Xb);                  % Luasan WL
        TPCd    = 1.025*Awld/100;                        % TPC

        xNoST   = [NoST*h,-sort(-NoST(1:end-1)*h)];
        yXb     = [Xb,-Xb(end-1:-1:1)];
        xNoSTa  = [NoST(1)*h NoST(1)*h xNoST];
        yXba    = [-Xb(1) Xb(1) yXb];
        [Area2,Area3] = TtkBerat(xNoSTa,yXba);
        
        yLuasST = [LuasST,-LuasST(end-1:-1:1)];
        yLuasSTa= [-LuasST(1) LuasST(1) yLuasST];
        [Area4,Area5] = TtkBerat(xNoSTa,yLuasSTa);
        
        %plot(xNoSTa,yXba,Area2(2),Area2(3),'or');
        KBd     = trapz(NoST*h,KbST.*LuasST)/Vold;        % KB
        Fmd     = Area2(2)-5*h;                           % F ke midship (- dibelakang midship)
        Bmd     = Area4(2)-5*h;                           % B ke midship (- dibelakang midship)
        BMd     = Area3(4)/Vold;                          % BM
        BMLd    = Area3(5)/Vold;                          % BML

        Wsad    = trapz(NoST*h,PanjGirthST);             % WSA 
        Lwld    = ( NoST(end)-NoST(1) )* h;              % Panjang LWL
        MTCd    = Displd * BMLd / (Lwld*100);            % MTC
        Cbd     = Vold / ( Lwld * Bd * d);               % Cb 
        Cmd     = Amd / ( Bd * d);                       % Cm 
        Cwd     = Awld / ( Lwld * Bd );                  % Cw 
        Cpd     = Vold / ( Lwld * Amd );                 % Cp

        Cb(k) = Cbd;
        Cm(k) = Cmd;
        Cw(k) = Cwd;
        Cp(k) = Cpd;
        Wsa(k)= Wsad;
        Vol(k)= Vold;
        Displ(k)= Displd;
        Awl(k) = Awld;
        TPC(k) = TPCd;
        KB(k)  = KBd;
        BM(k)  = BMd;
        KM(k)  = KBd+BMd;
        BML(k) = BMLd;
        KML(k) = KBd+BMLd;
        Fm(k)  = Fmd;
        MTC(k) = MTCd;
        Bm(k)  = Bmd;
    end
end

if pilih == 2
    clc;
    d_Cb_Cm_Cw_Cp=[d1',Cb',Cm',Cw',Cp'];
    d_Wsa_Awl_Vol_Disp_TPC=[d1',Wsa',Awl',Vol',Displ',TPC'];
    d_KB_BM_KM_BML_KML_Fm=[d1',KB',BM',KM',BML',KML',Fm'];
    
    
    d2  = linspace(0.01,d1(end),30);
    Cb2 = interp1(d1,Cb,d2,'linear','extrap');
    Cm2 = interp1(d1,Cm,d2,'linear','extrap');
    Cw2 = interp1(d1,Cw,d2,'linear','extrap');
    Cp2 = interp1(d1,Cp,d2,'linear','extrap');

    Wsa2 = interp1(d1,Wsa,d2,'linear','extrap');
    Awl2 = interp1(d1,Awl,d2,'linear','extrap');
    Vol2 = interp1(d1,Vol,d2,'linear','extrap');
    Disp2= interp1(d1,Displ,d2,'linear','extrap');


    KB2 = interp1(d1,KB,d2,'linear','extrap');
    BM2 = interp1(d1,BM,d2,'linear','extrap');
    KM2 = interp1(d1,KM,d2,'linear','extrap');
    BML2 = interp1(d1,BML,d2,'linear','extrap');
    KML2 = interp1(d1,KML,d2,'linear','extrap');
    TPC2 = interp1(d1,TPC,d2,'linear','extrap');
    MTC2 = interp1(d1,MTC,d2,'linear','extrap');

    Fm2 = interp1(d1,Fm,d2,'linear','extrap');
    Bm2 = interp1(d1,Bm,d2,'linear','extrap');

    figure(2);

    plot(Cb2,d2,'k');grid
    text(max(Cb2),max(d2)*0.95,'Cb','FontSize',8);
    cfig = get(gcf,'color');
    ax(1)=gca;
    pos = [0.12 0.45 0.75 0.5];
    offset = pos(4)/5.4;
    set(ax(1),'Position',pos,'FontSize',8)
    line(Cm2,d2,'Color','k','Parent',ax(1));text(max(Cm2),max(d2)*0.93,'Cm','FontSize',8);
    line(Cw2,d2,'Color','k','Parent',ax(1));text(max(Cw2),max(d2)*0.91,'Cw','FontSize',8);
    line(Cp2,d2,'Color','k','Parent',ax(1));text(max(Cp2),max(d2)*0.89,'Cp','FontSize',8);
    xlabel('Cb, Cm, Cw, Cp');ylabel('d (m)');

    pos2 = [pos(1) pos(2)-offset pos(3) pos(4)+offset];
    pos3 = [pos(1) pos(2)-2*offset pos(3) pos(4)+2*offset];
    pos4 = [pos(1) pos(2)-3*offset pos(3) pos(4)+3*offset];
    pos5 = [pos(1) pos(2)-4*offset pos(3) pos(4)+4*offset];

    scale2 = pos2(4)/pos(4);
    scale3 = pos3(4)/pos(4);
    scale4 = pos4(4)/pos(4);
    scale5 = pos5(4)/pos(4);
    limy1 = get(ax(1),'ylim');
    limy2 = [limy1(2)-scale2*(limy1(2)-limy1(1)) limy1(2)];
    limy3 = [limy1(2)-scale3*(limy1(2)-limy1(1)) limy1(2)];
    limy4 = [limy1(2)-scale4*(limy1(2)-limy1(1)) limy1(2)];
    limy5 = [limy1(2)-scale5*(limy1(2)-limy1(1)) limy1(2)];

    ax(2)=axes('Position',pos2,'Box','off',...
        'Color','none','YColor',cfig,...
        'XColor','r','Ytick',[],'Ylim',limy2,'Xlim',[0 inf],...
        'XAxisLocation','bottom','FontSize',8);
    lokWsa = interp1(d2,Wsa2,max(d2)*0.9);
    lokAwl = interp1(d2,Awl2,max(d2)*0.6);
    lokVol = interp1(d2,Vol2,max(d2)*0.95);
    lokDisp = interp1(d2,Disp2,max(d2)*0.8);
    line(Wsa2,d2,'Color','r','Parent',ax(2));text(lokWsa,max(d2)*0.9,'Wsa','FontSize',8,'Color','r');
    line(Awl2,d2,'Color','r','Parent',ax(2));text(lokAwl,max(d2)*0.6,'Awl','FontSize',8,'Color','r');
    line(Vol2,d2,'Color','r','Parent',ax(2));text(lokVol,max(d2)*0.95,'\nabla','FontSize',8,'Color','r');
    line(Disp2,d2,'Color','r','Parent',ax(2));text(lokDisp,max(d2)*0.8,'\Delta','FontSize',8,'Color','r');
    xlabel('Wsa ( m^2 ), Awl ( m^2 ), \nabla ( m^3 ), \Delta ( ton )');

    ax(3)=axes('Position',pos3,'Box','off',...
        'Color','none','YColor',cfig,...
        'XColor','b','Ytick',[],'Ylim',limy3,...
        'XAxisLocation','bottom','FontSize',8);
    lokKB = interp1(d2,KB2,max(d2)*0.9);
    lokBM = interp1(d2,BM2,max(d2)*0.85);
    lokKM = interp1(d2,KM2,max(d2)*0.6);
    lokTPC = interp1(d2,TPC2,max(d2)*0.75);
    lokMTC = interp1(d2,MTC2,max(d2)*0.4);
    line(KB2,d2,'Color','b','Parent',ax(3));text(lokKB,max(d2)*0.9,'KB','FontSize',8,'Color','b');
    line(BM2,d2,'Color','b','Parent',ax(3));text(lokBM,max(d2)*0.85,'BM','FontSize',8,'Color','b');
    line(KM2,d2,'Color','b','Parent',ax(3));text(lokKM,max(d2)*0.6,'KM','FontSize',8,'Color','b');
    line(TPC2,d2,'Color','b','Parent',ax(3));text(lokTPC,max(d2)*0.75,'TPC','FontSize',8,'Color','b');
    line(MTC2,d2,'Color','b','Parent',ax(3));text(lokMTC,max(d2)*0.4,'MTC','FontSize',8,'Color','b');
    xlabel('KB ( m ), BM ( m ), KM ( m ), TPC ( ton ), MTC ( ton.m )');

    ax(4)=axes('Position',pos4,'Box','off',...
        'Color','none','YColor',cfig,...
        'XColor','m','Ytick',[],'Ylim',limy4,...
        'XAxisLocation','bottom','FontSize',8);
    lokBML = interp1(d2,BML2,max(d2)*0.8);
    lokKML = interp1(d2,KML2,max(d2)*0.6);
    line(BML2,d2,'Color','m','Parent',ax(4));text(lokBML,max(d2)*0.8,'BML','FontSize',8,'Color','m');
    line(KML2,d2,'Color','m','Parent',ax(4));text(lokKML,max(d2)*0.6,'KML','FontSize',8,'Color','m');
    xlabel('BML ( m ), KML ( m )');

    warna = [ 0.502 0.000 1.000 ];
    ax(5)=axes('Position',pos5,'Box','off',...
        'Color','none','YColor',cfig,...
        'XColor',warna,'Ytick',[],'Ylim',limy5,'Xlim',[-5*h 5*h],...
        'XAxisLocation','bottom','FontSize',8);
    lokFm = interp1(d2,Fm2,max(d2)*0.4);
    lokBm = interp1(d2,Bm2,max(d2)*0.7);
    line(Fm2,d2,'Color',warna,'Parent',ax(5)); text(lokFm+0.5,max(d2)*0.4,'Fm','FontSize',8,'Color',warna);
    line(Bm2,d2,'Color',warna,'Parent',ax(5)); text(lokBm+0.5,max(d2)*0.7,'Bm','FontSize',8,'Color',warna);
    line([-5*h -5*h],[0 25],'Color','k','Parent',ax(5));
    xlabel('Fm ( m ), Bm ( m )');     
end


% --- Executes on button press in pbFile.
function pbFile_Callback(hObject, eventdata, handles)
% hObject    handle to pbFile (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global A UK

NamaFile = uigetfile('*.xls');
A    = xlsread(NamaFile);
UK  = A(1,:);

NamaFile = NamaFile(1:length(NamaFile)-4);
set(handles.text2,'String',NamaFile);
fclose('all');
set(handles.edit1,'String',strvcat(['Main Dimensions'],'-------------------------',...
    ['Loa (m) = ',blanks(3),num2str(UK(1))],...
    ['Lpp (m) = ',blanks(3),num2str(UK(2))],['B     (m) = ',blanks(3),num2str(UK(3))],...
    ['H     (m) = ',blanks(3),num2str(UK(4))],['T     (m) = ',blanks(3),num2str(UK(5))]));

% --- Executes on button press in pbHapus.
function pbHapus_Callback(hObject, eventdata, handles)
global Phideg pilih
% hObject    handle to pbHapus (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%if pilih == 3
%    for i= 1:length(Phideg)
%        close (figure(i));
%    end
%end
close (figure(1));close (figure(2));close (figure(3));
% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double


% --- Executes during object creation, after setting all properties.
function popProses_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popProses (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- Executes on selection change in popProses.
function popProses_Callback(hObject, eventdata, handles)
% hObject    handle to popProses (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popProses contents as cell array

%        contents{get(hObject,'Value')} returns selected item from popProses

function [ geom, iner, cpmo ] = TtkBerat( x, y ) 

[ x, ns ] = shiftdim( x );
[ y, ns ] = shiftdim( y );
[ n, c ] = size( x );

xm = mean(x);
ym = mean(y);
x = x - xm*ones(n,1);
y = y - ym*ones(n,1);


dx = x( [ 2:n 1 ] ) - x;
dy = y( [ 2:n 1 ] ) - y;


A = sum( y.*dx - x.*dy )/2;
Axc = sum( 6*x.*y.*dx -3*x.*x.*dy +3*y.*dx.*dx +dx.*dx.*dy )/12;
Ayc = sum( 3*y.*y.*dx -6*x.*y.*dy -3*x.*dy.*dy -dx.*dy.*dy )/12;
Ixx = sum( 2*y.*y.*y.*dx -6*x.*y.*y.*dy -6*x.*y.*dy.*dy ...
          -2*x.*dy.*dy.*dy -2*y.*dx.*dy.*dy -dx.*dy.*dy.*dy )/12;
Iyy = sum( 6*x.*x.*y.*dx -2*x.*x.*x.*dy +6*x.*y.*dx.*dx ...
          +2*y.*dx.*dx.*dx +2*x.*dx.*dx.*dy +dx.*dx.*dx.*dy )/12;
Ixy = sum( 6*x.*y.*y.*dx -6*x.*x.*y.*dy +3*y.*y.*dx.*dx ...
          -3*x.*x.*dy.*dy +2*y.*dx.*dx.*dy -2*x.*dx.*dy.*dy )/24;
P = sum( sqrt( dx.*dx +dy.*dy ) );

if A < 0,
  A = -A;
  Axc = -Axc;
  Ayc = -Ayc;
  Ixx = -Ixx;
  Iyy = -Iyy;
  Ixy = -Ixy;
end

xc = Axc / A;
yc = Ayc / A;
Iuu = Ixx - A*yc*yc;
Ivv = Iyy - A*xc*xc;
Iuv = Ixy - A*xc*yc;
J = Iuu + Ivv;

x_cen = xc + xm;
y_cen = yc + ym;
Ixx = Iuu + A*y_cen*y_cen;
Iyy = Ivv + A*x_cen*x_cen;
Ixy = Iuv + A*x_cen*y_cen;

I = [ Iuu  -Iuv ;
     -Iuv   Ivv ];
[ eig_vec, eig_val ] = eig(I);
I1 = eig_val(1,1);
I2 = eig_val(2,2);
ang1 = atan2( eig_vec(2,1), eig_vec(1,1) );
ang2 = atan2( eig_vec(2,2), eig_vec(1,2) );

geom = [ A  x_cen  y_cen  P ];
iner = [ Ixx  Iyy  Ixy  Iuu  Ivv  Iuv ];
cpmo = [ I1  ang1  I2  ang2  J ];


% --- Executes during object creation, after setting all properties.
function eKG_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eKG (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function eKG_Callback(hObject, eventdata, handles)
% hObject    handle to eKG (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of eKG as text
%        str2double(get(hObject,'String')) returns contents of eKG as a double

function StabGZ
global A UK Phideg KG GZ

ttk  = A(2:end,1:3);
k    = diff(ttk(:,1));
k2   = find(k ~= 0);
k2(length(k2) +1)   = length(ttk);
JumST   = length(k2);
z  = ttk(:,1); y = ttk(:,2); x = ttk(:,3);
h  = UK(2)/10;
d  = UK(5);

Phideg=linspace(-0.01,-89.9,19);

for pp = 1:length(Phideg)
j=0;
Phi= Phideg(pp)*pi/180;%figure(pp);
for i = 1: JumST
    x1 = x(find(z == z(k2(i)))); % mencari data tiap station
    y1 = y(find(z == z(k2(i))));
    z1 = z(find(z == z(k2(i))));
    
    y1(end+1) = y1(end) + 0.1;   % membuat chamber tiap station
    x1(end+1) = 0;
    z1(end+1) = z1(1);
    
    x2=x1(-sort(-find(x1(1:end-1)))); % mengurutkan data
    y2=-sort(-y1(2:end-1));
    
    x3=[x1(end);x2;x1(1);-x1(2:end)]; % data kanan dan kiri
    y3=[y1(end);y2;y1];
    z3=[z1(1)*ones(length(y3),1)];
    
    y6=-1:UK(4)*1.2;
    x6=zeros(1,length(y6));
    
    R = [cos(Phi),-sin(Phi),0; % data diputar
         sin(Phi),cos(Phi),0;
         0,0,1];
    x4=R(1,1)*x3+R(1,2)*(y3-d);
    y4=R(2,1)*x3+R(2,2)*(y3-d)+d;
    
    x7=R(1,1)*x6+R(1,2)*(y6-d);
    y7=R(2,1)*x6+R(2,2)*(y6-d)+d;
    
    x9=R(1,1)*0+R(1,2)*(KG-d);
    y9=R(2,1)*0+R(2,2)*(KG-d)+d;
    
    x51 = x4(find(y4 <= d));
    y51 = y4(find(y4 <= d));
    z51 = z3(find(y4 <= d));
    
    x5  = zeros(length(x51)+2,1);
    y5  = zeros(length(y51)+2,1);
    z5  = ones (length(z51)+2,1);
    
    lenky4 = length(y4);       % mencari perpotongan data ke garis air
    if Phideg == 90|Phideg == -90
        x31 = [x3(1);x3(find(x3 >= 0 & x3 ~= x3(end)))];
        y31 = [y3(1);y3(find(x3 >= 0 & x3 ~= x3(end)))];
        ny41 = x31(1);
        ny42 = x31(end);
    else     
        ny41    = find( (y4(1:lenky4-1) > d) & (y4(2:lenky4) < d) );
        ny42    = find( (y4(1:lenky4-1) < d) & (y4(2:lenky4) > d) );
    end
    ny=[ny41,ny42];
    if  ~isempty(ny)
        if Phideg == 90
            nx51 = find( (y4 == d) & (x4 ~= x4(1)) & (x4 ~= x4(end)));
            x5  = x4(nx51:end);
            y5  = y4(nx51:end);z5=z5*z51(1);
        elseif Phideg == -90
            nx51 = find( (y4 == d) & (x4 ~= x4(1)) & (x4 ~= x4(end)));
            x5  = x4(1:nx51);
            y5  = y4(1:nx51);z5=z5*z51(1);
        else
            xn1 = interp1(y4(ny41:ny41+1),x4(ny41:ny41+1),d);
            xn2 = interp1(y4(ny42:ny42+1),x4(ny42:ny42+1),d);
            x5(1)   = xn1; y5(1)   = d;
            x5(end) = xn2; y5(end) = d;z5=z5*z51(1);
            x5(2:end-1) = x51; y5(2:end-1) = y51;
        end
        
        %plot(x4,y4,x7,y7,'r'); hold on
        %[x5,y5]
        %plot(x5,y5,'r');hold on
    end
    
    if (length(y5) > 2 & length(x5) > 2) % hanya station yang punya dua titik atau lebih
        j=j+1;
        Area1   = TtkBerat( x5,y5 );
        LuasST(j)  = Area1(1);
        BxST(j)    = Area1(2);
        ByST(j)    = Area1(3);
        NoST(j)    = z5(1);
    end
    
end
%axis equal
%plot([-UK(3)*0.6,UK(3)*0.6],[UK(5),UK(5)],'--r');

Vold    = trapz(NoST*h,LuasST);                   % Volume carene
Bxd     = trapz(NoST*h,BxST.*LuasST)/Vold;  
Byd     = trapz(NoST*h,ByST.*LuasST)/Vold;
M       = interp1(x7,y7,Bxd);

%plot(x9,y9,'og',Bxd,y9,'ob'); hold on;            % Titik GZ
%plot(x6+Bxd,y6,'k',Bxd,Byd,'ob',Bxd,M,'or');      

%text(Bxd+0.1,Byd,'B1','FontSize',8,'Color','b');
%text(Bxd+0.1,M,'M','FontSize',8,'Color','r');

%text(x9+0.1,y9,'G','FontSize',8,'Color','g');
%text(Bxd+0.1,y9,'Z','FontSize',8,'Color','b');
GZ(pp) = Bxd - x9;
end
m=(GZ(2)-GZ(1))/(-Phideg(2));
xm = linspace(0,180/pi,20);
ym = m*xm;
figure(3);
plot(-Phideg,GZ,'-o',[0 90],[0 0],'k',xm,ym,'r',...
    [180/pi 180/pi],[0 m*180/pi],'r',180/pi,m*180/pi,'or');grid;
text(180/pi+ 3,m*180/pi,['GM = ',sprintf('%5.3f',m*180/pi),' m'],'FontSize',8,'Color','r');
pbaspect([2 1 1]);
xlabel('\phi ( deg )');
ylabel('GZ ( m )');

Contact us at files@mathworks.com