$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$     ----------     Model of Flag Pole ---------
$
$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$*********************************************      set basics
$
&channel gra_device  -page_dimen %pwidth% %pleng%
$
$*********************************************      set basics
$
&default -use @ $-nuse #buoy
&parameter -amctub 1 -drgtube 0 1e20 -f_cd_tube 0
&set diag     = .false.
&set n_nodes  = 11
&set twei     = 100
&set tfor     = 10
&set lene     = 10
&set lent     = 100
&set dia      = 1000
&set thick    = 25
&set Ch       = 2
&set g        = 9.8126                                   $ m/s^2
&set pi       = 3.14159
&set R        = &number(real .5*%dia/1000)
&set l1       = %lene
&set l2       = &number(real %lent-10)
&set EI       = &number(real 1.8188e6*%thick%/25)
&set rho_wat  = &number(real 1.025*%pi*%R**2 )
&set rho_stl  = &number(real %thick*53.63/(25*9.09*%g) )
&set rho_stl  = &number(real 1/590)
&set rho_bwat = &number(real %rho_stl+%rho_wat )
$
$*********************************************      list of periods
$
&set n_rao   = 51
&set p_incre = .25
&set p_beg   = 6.
&set liste   =
&loop i 1 %n_rao%
   &set liste = %liste% &number(real %p_beg%+((%i%-1)*%p_incre%) )
&endloop
$
$*********************************************      model of flag pole
$
&describe body f_pole
$
$*********************************************      nodes
$
&loop i 1 %n_nodes%
   &set where = &number(real -%lent%+%lene%*(%i%-1))
   *%i% 0 0 %where
   &type *%i% 0 0 %where
&endloop
$
$*********************************************      beams
$
~tube tube %dia% %thick% -buoy 0 -displ %rho_stl*%g -wtp %rho_stl*%g
&loop i 1 %n_nodes%-1
   &set j = &number(real %i%+1)
   beam ~tube *%j% *%i%
   &type beam ~tube *%i% *%j%
&endloop
$
$*********************************************      restraints
$
~spr1  spr x 1e10 y 1e10  z 1e10 rz 1e10
~spr2  spr x 1e10 y 1e10
$
$*********************************************      weight @ top
$
&describe load_group top_weig *%n_nodes%
#weight *%n_nodes% %twei% 5 5 5
$
$*********************************************      force  @ top
$
&describe load_group top_for *%n_nodes%
#f_top -%tfor -cen 0 0 0
$
$*********************************************      interest points
$
&describe interest -ass *%n_nodes
&macro c_integrals
$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$              ----------     c_integrals ---------
$
$ This macro computes the constants so the a two piece cubic equation
$ fits the boundary conditions. It also computes the "energy" integrals
$ x^2, x^2, and xpp^2, where xp is the derivative of the cubic, and
$ xpp is its second derivative.
$
$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$***************************************       zeta
$
&set zeta   = &number(real -1/(2*(%l2**2)*(%l2+%l1)) )
$
$***************************************       alpha
$
&set alpha  = &number(real %zeta*%l2*%l1 )
$
$***************************************       beta
$
&set beta   = &number(real -%zeta*%l2/%l1 )
$
$***************************************       gamma
$
&set gamma = &number(real (1/%l2)+(2*%zeta*(%l2**2)) )
$
$***************************************       delta
$
&set delta  = &number(real -3*%zeta*%l2 )
$
$***************************************       int x^2
$
&set int_x1   = &number(real (%alpha**2)*(%l1**3)/3 )
&set int_x2   = &number(real 2*%alpha*%beta*(%l1**5)/5 )
&set int_x3   = &number(real (%beta**2)*(%l1**7)/7 )
&set int_x4   = &number(real (%gamma**2)*(%l2**3)/3 )
&set int_x5   = &number(real 2*%gamma*%delta*(%l2**4)/4 )
&set int_x6   = &number(real (1/5)*(2*%gamma*%zeta+(%delta**2))*(%l2**5) )
&set int_x7   = &number(real (2/6)*%delta*%zeta*(%l2**6) )
&set int_x8   = &number(real (%zeta**2)*(%l2**7)/7 )
&set int_xa   = &number(real %int_x1+%int_x2+%int_x3 )
&set int_xb   = &number(real %int_x4+%int_x5+%int_x6+%int_x7+%int_x8 )
&set int_x    = &number(real %int_xa+%int_xb )
$
$***************************************       int xp^2
$
&set int_xp1  = &number(real (%alpha**2)*%l1 )
&set int_xp2  = &number(real (6/3)*%alpha*%beta*(%l1**3) )
&set int_xp3  = &number(real (9/5)*(%beta**2)*(%l1**5) )
&set int_xp4  = &number(real (%gamma**2)*%l2 )
&set int_xp5  = &number(real (4/2)*%gamma*%delta*(%l2**2) )
&set int_xp6  = &number(real (1/3)*(6*%gamma*%zeta+4*(%delta**2))*(%l2**3) )
&set int_xp7  = &number(real (12/4)*%delta*%zeta*(%l2**4) )
&set int_xp8  = &number(real (9/5)*(%zeta**2)*(%l2**5) )
&set int_xpa  = &number(real %int_xp1+%int_xp2+%int_xp3 )
&set int_xpb  = &number(real %int_xp4+%int_xp5+%int_xp6+%int_xp7+%int_xp8 )
&set int_xp   = &number(real %int_xpa+%int_xpb )
$
$***************************************       int xpp^2
$
&set int_xpp1 = &number(real (36/3)*(%beta**2)*(%l1**3) )
&set int_xpp2 = &number(real 4*(%delta**2)*%l2 )
&set int_xpp3 = &number(real (24/2)*%delta*%zeta*(%l2**2) )
&set int_xpp4 = &number(real (36/3)*(%zeta**2)*(%l2**3) )
&set int_xppa = &number(real %int_xpp1 )
&set int_xppb = &number(real %int_xpp2+%int_xpp3+%int_xpp4 )
&set int_xpp  = &number(real %int_xppa+%int_xppb )
$
$***************************************       diagnostic
$
&if %diag  &then
   &type int_x1 = %int_x1
   &type int_x2 = %int_x2
   &type int_x3 = %int_x3
   &type int_x4 = %int_x4
   &type int_x5 = %int_x5
   &type int_x6 = %int_x6
   &type int_x7 = %int_x7
   &type int_x8 = %int_x8
   &type ====================
   &type int_xp1 = %int_xp1
   &type int_xp2 = %int_xp2
   &type int_xp3 = %int_xp3
   &type int_xp4 = %int_xp4
   &type int_xp5 = %int_xp5
   &type int_xp6 = %int_xp6
   &type int_xp7 = %int_xp7
   &type int_xp8 = %int_xp8
   &type ====================
   &type int_xpp1 = %int_xpp1
   &type int_xpp2 = %int_xpp2
   &type int_xpp3 = %int_xpp3
   &type int_xpp4 = %int_xpp4
   &type ====================
   &type alpha = %alpha ; beta = %beta ; gamma = %gamma
   &type delta = %delta ; zeta = %zeta
   &type ====================
   &type int_x = %int_x ;int_xp = %int_xp ;int_xpp = %int_xpp ;
&endif
$
$***************************************       all done
$
&endmacro
&macro e_estimate
$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$              ----------     e_estimate   ---------
$
$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$***************************************       evaulate integrals
$
c_integrals
$
$***************************************       solve for A B and D
$
&set A   = &number(real (1/2)*%EI*%int_xpp )
&set B   = 0
&set Dba = &number(real (1/2)*%rho_stl*%int_x )
&set Dbw = &number(real (1/2)*%rho_bwat*%int_x )
&set Dw  = &number(real %twei/(2*%g) )
$
$***************************************       solve for freq
$
&report
&set rho = %rho_stl
type =========================================================
type |                                                       |
type |      Estimates of Natural Frequency, Deflection       |
type |      ==========================================       |
type |                                                       |
$
$***************************************       deflection
$
&set defl = &number(real 10/(2*%a) )
&set defl = &format(f7.4 %defl)
type | Deflection Under 10 kn Horizontal force = %defl m    |
type |                                                       |
$
$***************************************       nonlinear defl
$
&set B      = &number(real (1/2)*%twei*%int_xp )
&set defl = &number(real 10/(2*%a-2*%b) )
&set defl = &format(f7.4 %defl)
type | Defl. - 10 kn Horiz force, 100kn Top W .= %defl m    |
type |                                                       |
$
$***************************************       omega in air
$
&set omega2 = &number(real %A%/%Dba% )
&set omega  = &number(sqrt  %omega2  )
&set period = &number(real (2*%pi%)/%omega )
&set period = &format(f8.4 %period)
type | Natural Period in Air                   = %period Sec. |
type |                                                       |
$
$***************************************       omega in water
$
&set omega2 = &number(real %A%/%Dbw% )
&set omega  = &number(sqrt  %omega2  )
&set period = &number(real (2*%pi%)/%omega )
&set period = &format(f8.4 %period)
type | Natural Period in water                 = %period Sec.|
type |                                                       |
$
$***************************************       in water, weight, linear
$
&set B      = 0
&set omega2 = &number(real (%A%-%b)/(%Dbw%+%dw%) )
&set omega  = &number(sqrt  %omega2  )
&set period = &number(real (2*%pi%)/%omega )
&set period = &format(f8.4 %period)
type | Period in water, 100kn Weight, Linear   = %period Sec.|
type |                                                       |
$
$***************************************       in water, weight, nonlin
$
&set B      = &number(real (1/2)*%twei*%int_xp )
&set omega2 = &number(real (%A%-%b)/(%Dbw%+%dw%) )
&set omega  = &number(sqrt  %omega2  )
&set period = &number(real (2*%pi%)/%omega )
&set period = &format(f8.4 %period)
type | Period in water, 100kn Weight, Nonlin.  = %period Sec.|
$
$***************************************       finish box
$
type |                                                       |
type =========================================================
end
$
$***************************************       all done
$
&endmacro
&macro mor_work weight -nonlinear nonlin
$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$              ----------     mor_work   ---------
$
$
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$
$
$***************************************       initialize
$
&set l_drao = 1
$
$***************************************       constants
$
&set J   = &number(real %Ch*%rho_wat )
&set A   = &number(real (1/2)*%EI*%int_xpp )
&if %nonlin &then
   &set B   = &number(real (1/2)*%twei*%int_xp )
&else
   &set B   = 0
&endif
&set Dba = &number(real (1/2)*%rho_stl*%int_x )
&set Dbw = &number(real (1/2)*%rho_bwat*%int_x )
&set Dw  = &number(real %weight/(2*%g) )
&set DU  = &number(real %dbw+%dw)
&set Ds  = &number(real %a-%b)
$
$***************************************       message
$
&report  'Raos with Top Weight = %weight'
type
type   Period          RAO Deflection
type ========================================
$
$***************************************       loop periods
$
&loop t_rao ( &reverse(%liste) )
$
$***************************************       factors
$
   &set w_rao    = &number(real 2*%pi/%t_rao )
   &set w_rao2   = &number(real %w_rao*%w_rao )
   &set k        = &number(real %w_rao2/%g )
   &set u1       = &number(real %k*%l1 )
   &set u2       = &number(real %k*%l2 )
   &set exp_u1   = &number(exp  %u1    )
   &set exp_u2   = &number(exp  %u2    )
   &set exp_nu2  = &number(exp  -1*%u2    )
   &set exp_nu3  = &number(exp  -1*%k*%lent )
   &set fac11    = &number(real %j*%w_rao2*%exp_nu3*(%l1/%u1)**2)
   &set fac12    = &number(real %fac11*(%l1/%u1))
   &set fac21    = &number(real %j*%w_rao2*%exp_nu2*(%l2/%u2)**2)
   &set fac22    = &number(real %fac21*(%l2/%u2))
   &set fac23    = &number(real %fac22*(%l2/%u2))
$
$***************************************       integrals
$
   &set I1 = &number(real ((%u1-1)*%exp_u1+1))
   &set I2 = &number(real (((%u1-2)*%u1+2)*%exp_u1-2) )
   &set I3 = &number(real ((%u2-1)*%exp_u2+1) )
   &set I4 = &number(real (((%u2-2)*%u2+2)*%exp_u2-2) )
   &set I5 = &number(real (((%u2**2-3*%u2+6)*%u2-6)*%exp_u2+6) )
$
   &set az  = &number(real %alpha*%fac11*%I1   )
   &set bz  = &number(real %beta*%fac12*%I2   )
   &set gz  = &number(real %gamma*%fac21*%I3   )
   &set dz  = &number(real %delta*%fac22*%I4   )
   &set zz  = &number(real %zeta*%fac23*%I5   )
   &set num_rao = &number(real %az+%bz+%gz+%dz+%zz )
$
$***************************************       deflection rao
$
   &set k   = %a-%b
   &set den_rao  = &number(real 2*(%a-%b-%w_rao2*%DU) )
   &set den_rao  = &number(sqrt %den_rao%**2+.0001*4*%du*%ds)
   &set d_rao    = &number(abs %num_rao/%den_rao )
   &set l_drao   = %l_drao %w_rao %d_rao
$
$***************************************       display results
$
&set per_whi = &format(f11.3 %t_rao  )
&set dra_whi = &format(f11.3 %d_rao  )
type     %per_whi         %dra_whi
$
$***************************************       end loop
$
&endloop
end
$
$***************************************       all done
$
&endmacro
$
$***************************************       evaluate
$