cb=cbind(c(-50,-35,-5,12,42,57),c(-10,-10,-10,-10,-10,-10))
ch=cbind(c(-50,-42,-10,17,49,57),c(-4,-5,-5.5,-5.5,-5,-4))

moorings=function(vw=92.6,sidearea=950,alpha=pi/2,cw=1.3,lwl=117,coordboll=cb,
coordch=ch,crossarea=90,ei=125,ropetension=25)
{
#Computation of environmental forcing and mooring forces
#for a berthed ship

#Procedure is followed as in https://www.albertomontanari.it/?q=node/256

#Computation of wind force
py=cw*sidearea*sin(alpha)^2*0.01225*(vw/3.6)^2/(2*9.81)
#Computation of rotational moment positive clockwise
mw=-0.04*lwl*py

#Computation of angles theta
theta=rep(0,6)
for (i in 1:6)
{
theta[i]=atan((cb[i,1]-ch[i,1])/(ch[i,2]-cb[i,2]))
}

#Computation of the length of the mooring lines
lbc=rep(0,6)
for (i in 1:6)
{
#Length from bollard to chock
lbc[i]=(ch[i,2]-cb[i,2])/cos(theta[i])
#Total length mooring line
ltot=lbc*1.15
}

#Computation of spring constant
ki=crossarea*ei/ltot

#Determine spring component along y
kyi=ki*cos(theta)

#Determine kyi*xch,i and kyi*xch,i(squared)
kyixch=kyi*ch[,1]
kyixch2=kyi*ch[,1]^2

#Compute A, B and C
A=sum(kyi)
B=sum(kyixch)
C=sum(kyixch2)

#Calculate ship's translation in the y direction
sigmay=(py*C-mw*B)/(A*C-B^2)

#Calculate rotation of ship around the vertical axis
rotation=(py*B-mw*A)/(B^2-A*C)

#For each line calculate its force component
pyi=kyi*(sigmay+ch[,1]*rotation)

#Compute forces in the direction of the mooring line
ti=pyi/cos(theta)

#Determine safety factors for each mooring line, FS = BS / Ti.
FS=ropetension*crossarea/ti

if(min(FS)>10) checksafety="T"

#Check the computation by verifying that ∑iPy,i = Py.
if ((sum(pyi)-py)<0.1) checkforces="T"

output=list(py=py,mw=mw,pyi=pyi,FS=FS,checkforces=checkforces,checksafety=checksafety)
return(output)
}