/*
Parametric study with zrbspan01.mac to account for elastic-perfectly-plastic zrb support.
Must run zrbspan01.mac before running this routine.  All inputs except cka and xias are those
set in zrbspan01.mac.  Additional inputs are listed below.

Inputs:  */
umas:0.002$          /* ZRB support mobilization displ (m) */
fz:0.068882912$      /*  static vertical support reaction (MN) */
muas:0.5$            /*  axial friction coefficient */

/* Referenced Notes 
 {1} Input dynamic axial displacements at ZRB support may not be less than the mobilization displacement.
     I.e. ncaud >= 1.  Otherwise the expressions for cka and xias would need to be adapted.
     When this is not the case, the problem is linear and the same result as for ncaud=1 should be obtained
     for frq_ and zeta, and the amplitudes can be any fraction of those for ncaud=1.
     To properly write the results to the .mo file, the first value of ncaud MUST be 1.
*/


/* derived inputs */
fmas:fz*muas$            /* axial slip force */
caud:float(ncaud*umas)$  /* dyn ax displ @ X=0 */
cka:float(fmas/caud)$
xias:float((2/%pi)*(1-umas/caud))$

/* required statements extracted from mzrspan01.mac */
cnf:cgg*ea/olea$  /* unit amplitude value of capital n sub f in notes */
cnd:cgg/(2/cka + olea/ea)$  /* unit amplitude value of capital n sub d in notes */
cuab:0.5*cnd^2/sqrt(ea*ka)$  /* axial seabed springs */
cuap:cl_*cnd^2/ea+cuab$  /* axial in pipe */
cuas:2*cnd^2/cka$
cua:cuab+cuas+cuap$
frqa_:float(sqrt(cua/(meff*ew_)))/(2*fpi)$
frq_:sqrt(frqu_^2+frqa_^2)$
cuk:(2*%pi*frq_)^2*meff*ew_$  /* unit-amplitude kinetic energy */
zeta_as:xias*cuas/cuk$
zeta_a:zeta_ab+zeta_as+zeta_ap$
zeta:zeta_a+zeta_u$
ud0:2*cnd/cka$   /*  modal axial displ at support */
aod:caud/(odwc*ud0)$

/*plot some results */
plot2d([parametric,subst(ncaud=t,aod),subst(ncaud=t,ncaud)],[t,1,8],
         [gnuplot_preamble,"set key bottom"],
         [xlabel,"Normalized Maximum Amplitude of Vibration, A/D"],
         [ylabel,"Amplitude of Axial Displacement at ZRB Support / Mobilization Displacement"])$
plot2d([parametric,subst(ncaud=t,aod),subst(ncaud=t,1000*caud)],[t,1,8],
         [gnuplot_preamble,"set key bottom"],
         [xlabel,"Normalized Maximum Amplitude of Vibration, A/D"],
         [ylabel,"Amplitude of Axial Displacement at ZRB Support (mm)"])$
plot2d([parametric,subst(ncaud=t,aod),subst(ncaud=t,frq_)],[t,1,8],
         [gnuplot_preamble,"set key bottom"],
         [xlabel,"Normalized Maximum Amplitude of Vibration, A/D"],
         [ylabel,"Natural Frequency (Hz)"])$
plot2d([parametric,subst(ncaud=t,aod),subst(ncaud=t,zeta)],[t,1,8],
         [gnuplot_preamble,"set key bottom"],
         [xlabel,"Normalized Maximum Amplitude of Vibration, A/D"],
         [ylabel,"Modal Damping Ratio"])$

/* write results to a zrbspan02.mo */
file_output_append:true$
wrl:["aod","caud","frq_","xias","zeta_as","zeta"]$
write_data(wrl,"C:/Users/Ralf Peek/Documents/w/spans/zrbspan02.mo")$
t:1$
wrl:[0,0,
    float(subst(ncaud=t,frq_)),
    float(subst(ncaud=t,xias)),
    float(subst(ncaud=t,zeta_as)),
    float(subst(ncaud=t,zeta))
    ]$
write_data(wrl,"C:/Users/Ralf Peek/Documents/w/spans/zrbspan02.mo")$
for t:1 step 0.02 thru 8 do block(
  wrl:[
    float(subst(ncaud=t,aod)),
    float(subst(ncaud=t,caud)),
    float(subst(ncaud=t,frq_)),
    float(subst(ncaud=t,xias)),
    float(subst(ncaud=t,zeta_as)),
    float(subst(ncaud=t,zeta))
    ],
  write_data(wrl,"C:/Users/Ralf Peek/Documents/w/spans/zrbspan02.mo")
)$