function Model = wecModel1Body(paramFile)
% wecModel1Body
%   Model = wecModel1Body(paramFile)
%   This function creates a model in frequency domain for the heave motion
%   for WECs type 1 or type 2.
%   It returns the Model in a structure with the fields mag, phase and
%   frequency
%
%   Project: IMPACT
%   Title of Code: wecModel1Body
%   Relevant additional documentation:
%   Developer: Carlos Gonzalez
%   External Dependencies: NA
%
%   Disclaimer
%   CAUTION: It is strongly recommended that changes to this code
%   should be avoided without first consulting with authors.

% Run parameters file
if ~exist(paramFile,'file')
    errordlg(sprintf('File %s cannot be found',paramFile),'Error')
    Model=[];
    return;
end

[iSuccess, stErrorMsg] = readFile(paramFile);
if ~iSuccess
    errordlg(stErrorMsg,'Error');
    Model=[];
    return;
end

% Read BEM parameters
if ~exist(stH5File,'file')
    errordlg(sprintf('File %s cannot be found',stH5File),'Error')
    Model=[];
    return;
end
h5Data = getParamH5(stH5File,1);
[n,m,p] = size(h5Data(1).addedMassC);

arModelFreq = h5Data(1).BEMFreq;
numFreq = numel(arModelFreq);

% Body 1
A1BEMnorm = h5Data(1).addedMassC;
B1BEMnorm = h5Data(1).radDampC;
C1BEMnorm = h5Data(1).hydroStiffC;

A1BEM = A1BEMnorm * rDENSITY;
B1BEM = B1BEMnorm * rDENSITY;
C1BEM = C1BEMnorm * rDENSITY*rGRAVITY;

A1_1DOF = reshape(A1BEM(:,3,3),1,n);
B1_1DOF = reshape(B1BEM(:,3,3),1,n);
C1 = C1BEM(3,3);


% Model from pto force to relative displacement
Z = (-arModelFreq.^2.*(rFloaterMass+A1_1DOF)+ C1) + 1i*(arModelFreq.*(B1_1DOF + rAddedDampingBody1));

magZ = abs(Z);
phaseZ = angle(Z);
[PTOBodeMag,PTOBodePhase,freqPTO] = bode(tfPto,arModelFreq);
PTOMag =reshape(PTOBodeMag,1,numFreq);
PTOPhase = reshape(PTOBodePhase,1,numFreq);

modelPto.arMagnitude = 1./magZ.*PTOMag;
modelPto.arPhase = wrapAngle(-phaseZ*180/pi+PTOPhase);
modelPto.arFreq = arModelFreq;


% Assign output
Model.zModelPto = modelPto;
Model.zModelFloater = modelPto;
