`Free roll decay time domain simulation
clear all
clear

`1) Generate hydrodynamic data using SEAkeep
read fv
we 250 0 0 5 /gyr:gf,gf,gf
so

wave (spe) p1 3.5 `this really doesn't matter in this case, just need a wave defined

report sk.$$$
 sea /brief /head:90 /samp:50,0.1,2.0 /data:hy `sample 50 samples from 0.1 to 2.0 rad/sec
report off | erase sk.$$$

`2) Parse the hydrodynamic data file using SK.LIB
run sk.lib /call /quiet 
.sk.gethydrodata "sk-hydros.dat" `see .sk.libinfo

`3) Compute hydrodynamic damping and added mass at the natural frequency
vari wn44,wn44p,a44n,b44n,b44vn,b44e,c44n,am44,f1,f2,a1,a2,b1,b2,bv1,bv2
vari RAD=0.01745 `PI/180

macro NaturalFreq
 `Recall: wn = SQRT( k / m )
 `...where k is the stiffness and m is the total mass (physical + added)
 `%1 is physical mass
 `%2 is added mass
 `%3 is stiffness
 `%4 is returned variable name
 %4:=SQRT(%3/(%1+%2))
/

macro Interpolate `linear interpolation
 `%1,%2 is x1,y1
 `%3,%4 is x2,y2
 `%5 is x3, y3 returned in %6 
 vari m,b
 m:=(%4-%2)/(%3-%1)
 b:=%2-m*%1
 %6:=m*%5+b
/

macro FindABC `compute the added mass, damping, and stiffness at the natural frequency
 i:=i+1
 wn44p:=wn44
 f1:=sk.freq_{i}  | f2:=sk.freq_{i+1}
 a1:=sk.a44_{i}   | a2:=sk.a44_{i+1}
 b1:=sk.b44_{i}   | b2:=sk.b44_{i+1}
 bv1:=sk.b44v_{i} | bv2:=sk.b44v_{i+1}
 if {sk.freq_{i+1}}>={wn44p} then if {sk.freq_{i}}<={wn44p} then ;;
    .Interpolate {f1} {a1} {f2} {a2} {wn44p} a44n ;;
  | .Interpolate {f1} {b1} {f2} {b2} {wn44p} b44n ;;
  | .Interpolate {f1} {bv1} {f2} {bv2} {wn44p} b44vn ;;
  | .NaturalFreq {sk.ixx} {a44n} {sk.c44} wn44 ;;
  | if {wn44-wn44p}<0.1 then exit `iterate until convergence
 if {i+1}>{sk.nwaves} then exit else exit FindABC `exit if frequency not found in sampled range
/

vari i=0
.NaturalFreq {sk.ixx} {sk.a44_{sk.nwaves}} {sk.c44} wn44 `compute initial guess
.FindABC
am44:=sk.ixx+a44n `compute total mass at natural frequency (physical plus added mass)
b44e:=b44n+b44vn `compute total linearized damping coefficient (including viscous term)

`critical damping ratio
vari zeta,b44c
b44c:=2*am44*wn44 `critical damping
zeta:=b44e/b44c

report RollDecay.PF
\ROLL DECAY SIMULATION\
\
\Undamped Roll Natural Frequency (Rad/Sec): {wn44:3} \
\
\Added Mass (Slug-Ft^2): {a44n:1}  Total Mass (Slug-Ft^2): {am44:1}\
\Equivalent Linear Damping Coefficient (Slug-Ft^2/Sec): {b44n:1}\
\Stiffness Coefficient (Lb-Ft): {sk.c44:1}\
\
\Critical Damping Ratio: {zeta:4}\

`4) Solve the uncoupled roll equation of motion in the time domain
vari v1,v1p,v2,v2p,v2d,v2dp
vari endtime=60 `simulation time, 60 seconds
vari t,dt=0.05 `timestep

macro Integrator `implements a first order explicit Euler scheme, so watch your timestep :)
 t:=t+dt | v2p:=v2 | v1p:=v1
 v2dp:=(-1/am44)*(sk.c44*v1p+b44e*v2p) `compute acceleration
 v2:=v2p+dt*v2dp `compute velocity
 v1:=v1p+dt*v2 `compute position
 me {t} {v1/RAD}
 if {t}>{endtime} then exit else exit Integrator
/

t:=0
v1:=20*RAD `initial condition 20 degree roll
v2:=0 `zero initial velocity

\
me plotstart "Free Roll Decay" /norot /page%:50
me plotlabel "Time (s)","Roll Angle (deg)"
me {t} {v1/RAD}
.Integrator
me plotend
report close /preview