function [panelsL,panelsR,panels0,panelsB,panelsF,Nfl,Nfr,Nfm] = ...
    semicircles_two(Dl,Dr,al,ar,l,r,h,Delta,scale)
% INPUT
% calculates the nodes for a single plate
% D, the extent of the plate
% a - the offsets relative to the origin of the centre of the plate
% l and r the left and right ends of the domain
% h is the water depth
% Delta is the (approximate) panel length
% Scale is an optional command which scales the circle (so that scale = 0
% corresponds to a dock)

%OUTPUT
% panelsL -- panels for left control surface
% panelsR -- panels for right control surface
% panels0 -- panels for the body
% panelsB -- body panels
% panelsF -- free-surface panels
% scale is optional
if nargin == 8
    scale = 1;
end


% first we calculate the number of nodes
Nbr = max(round(pi/2*Dr/Delta),1); % nodes on right semicircle
Nbl = max(round(pi/2*Dl/Delta),1); % nodes on left semicircle
Nfl = max(round((-l + al - Dl/2)/Delta),0); % nodes on left free surface (these can be zero)
Nfr = max(round((r - ar - Dr/2)/Delta),0); % nodes on right free surface (these can be zero)
Nfm = max(round((ar - al - (Dr+Dl)/2)/Delta),1); % nodes on middle free surface
Nc  = max(round(h/Delta),1); % nodes on vertical surfaces
Nbed = max(round((r+-l)/Delta),1); % nodes on vertical surfaces

[bodypanels_r,xrl,xrr]=body_shape(Nbr,Dr,ar,scale);
[bodypanels_l,xll,xlr]=body_shape(Nbl,Dl,al,scale);


%%%%%%% Creat intermediate free surface panels %%%%%%%%%%%%%%%%%%%%%%%%%%%


dxr=(xrr-r)/(Nfr);
dxm=(xlr-xrl)/(Nfm);
dxl=(l-xll)/(Nfl);
xr=r:dxr:xrr;  % Nfr+1 nodes, Nfr panels
xm=xrl:dxm:xlr;  % Nfm+1 nodes, Nfm panels
xl=xll:dxl:l; % Nfl+1 nodes, Nfl panels
zf=0;
z1 = 0:-h/Nc:-h;
xbed=l:(r+-l)/Nbed:r;
z2=-h:h/Nc:0;

% Panels L
panelsL=[xbed(1)*ones(Nc,1),z1(1:Nc).',xbed(1)*ones(Nc,1),z1(2:Nc+1).'];
% Panels R
panelsR=[xbed(Nbed+1)*ones(Nc,1),z2(1:Nc).',xbed(Nbed+1)*ones(Nc,1),z2(2:Nc+1).'];
% Panels Rigid Surface
panels0=[xbed(1:Nbed).',z1(Nc+1)*ones(Nbed,1),xbed(2:Nbed+1).',z1(Nc+1)*ones(Nbed,1)];
panelsB=[bodypanels_r;  bodypanels_l];
% Panels F
panelsF=[ xr(1:Nfr).',zf*ones(Nfr,1),xr(2:Nfr+1).',zf*ones(Nfr,1);...
    xm(1:Nfm).',zf*ones(Nfm,1),xm(2:Nfm+1).',zf*ones(Nfm,1);...
    xl(1:Nfl).',zf*ones(Nfl,1),xl(2:Nfl+1).',zf*ones(Nfl,1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function [bodypanels,xl,xr]=body_shape(Nb,D,a,scale)

        z_body=zeros(1,Nb);    % initialisation

        %your body shape
        for i=1:Nb+1,
            z_body(i) = a+(D/2)*(cos((i-1)*pi/Nb)-j*sin((i-1)*pi/Nb));
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        bodypanels=[real(z_body(1:Nb))',scale*imag(z_body(1:Nb))',...
            real(z_body(2:Nb+1))',scale*imag(z_body(2:Nb+1))'];
        xl=bodypanels(Nb,3);
        xr=bodypanels(1,1);


    end

    function [N]=bodypanels(eqbodynode)

        %load('nodes.in', Nb, Nb1, Nb2)
        if eqbodynode==1,
            Nb=input('Enter the number of panels on each body: ')
        else
            Nbr=input('Enter Nbr, number of panels on the right body: ')
            Nbl=input('Enter Nbl, number of panels on the left body: ')
        end

        % free-surface nodes
        [Nf]=input('Enter Nfl, Nfm, Nfr free-surface panel numbers: ')

        Nfl=Nf(1);
        Nfm=Nf(2);
        Nfr=Nf(3);

        Nc=input('Enter Nc, the number of panels on the control surfaces: ')
        Nbed=input('Enter the number of panels on the sea bed: ')
        N=[Nc,Nbed,Nc,Nfl,Nb,Nfm,Nb,Nfr]
    end


end