restart:with(linalg):with(plots):

load parameter omega  ( w nondimensional)

the CW angle beta ( b in degrees) between the x-axis and the direction of force  

the CW angle alpha (a in degrees) between tangent to the beam base curve at the tip point and direction of the terminal force 

Follower load -  the load parameter w and angle a are given. This problem has one solution.

Load parameter problem - both a and b are given. This problem has infinitely  many solutions.

Conservative load - load parameter w and angle b are given. This problem has finite number of solutions.

#========================================

xi := proc( s, omega) global k, K, K1, E;  if omega = 0 then (1 - s)*(1 - 2*k**2) else  (2*E*K1-1)*(1 - s) + 2/omega*(JacobiZeta(omega,k) - JacobiZeta( omega*s, k) - k**2*(JacobiSN(omega,k)*JacobiCN(omega,k)/JacobiDN(omega,k)-JacobiSN(omega*s,k)*JacobiCN(omega*s,k)/JacobiDN(omega*s,k))) end if; end proc;

#========================================

eta := proc( s, omega) global k, k1, K; if omega = 0 then (1 - s)*2*k*k1 else 2*k*k1/omega*(JacobiSN(omega,k)/JacobiDN(omega,k) - JacobiSN(omega*s,k)/JacobiDN(omega*s,k)) end if; end proc;

#========================================

#========================================

x:=proc( s, omega) local xx, yy, beta0; global k, K; beta0:=beta(omega); xx:=xi( s, omega); yy:= eta(s, omega);xx*cos(beta0)+yy*sin(beta0);end proc;

#========================================

y:=proc( s, omega) local xx, yy, beta0; global k, K;beta0:=beta(omega); xx:=xi( s, omega); yy:= eta(s, omega);-xx*sin(beta0)+yy*cos(beta0);end proc;

#==============================================

#==============================================

beta:=proc( omega) global k; 2*arcsin(k*JacobiCN(omega,k)/JacobiDN(omega,k)) end proc;

#==============================================

#===================================================

setup := proc( alpha) global k, k1, K, K1, E; k:=evalf(sin(rad(alpha)/2)); k1:= sqrt(1 - k**2); K:=EllipticK(k); E:=EllipticE(k); if k < 1 then K1 := 1/K else K1 := 0 end if; end proc;

#====================================================

#===================================================

force:=(f,beta0)-><-f*cos(beta0),f*sin(beta0)>;

#====================================================

# Utility functions

#=======================

rad:=alpha->evalf(alpha*Pi/180): # convert deg to rad

deg:=alpha->evalf(alpha*180/Pi): # convert rad to deg

#====================================

number of points that is used for approximation of deformed beam

NumPts:=100; # number of points for beam shape

number of points used for approximation of beam tip path

NumPtsPath:=180;

the scale factor for the displaying  the applied force

ForceScale:=0.025; # Scale for force

 number of frames used for animations

NumFrames:=50; # number of frames for animations

Digits:=14;

#========================================

alpha0:=150.; #Input

setup(alpha0): # Now K is availabe

omega0:=5.3*K; # Input

#ForceScale:=0.025:   # scale for force

#NumPts:=100:  # number of points

beta0:=beta(omega0): # we need this to show beam root

Beam:=plot([xi(s,omega0),eta(s,omega0),s=0..1],numpoints = NumPts, color=black, thickness=2, scaling=constrained,caption=typeset("omega = ", omega0," alpha = ",alpha0)): # Local coordinates

ForceVect:= arrow(<xi(0,omega0),eta(0,omega0)>,<-ForceScale,0>):

BeamRoot:=pointplot( [[ForceScale*sin(beta0),-ForceScale*cos(beta0)],[-ForceScale*sin(beta0),ForceScale*cos(beta0)]],connect=true):

display( Beam, ForceVect,BeamRoot);

#-------------------------------------

#-------------------------

ForceVect:= arrow(<x(0,omega0),y(0,omega0)>,force(ForceScale,beta0)):

Beam:=plot([x(s,omega0),y( s , omega0),s=0..1],numpoints = NumPts, color=black, thickness=2, scaling=constrained,caption=typeset("omega = ", omega0," alpha = ",alpha0,"  beta = ",evalf(beta0*180/Pi))): # Space coordinates

display( Beam, ForceVect);

#========================================

#  Follower load problem

#========================================

# Constant load parameter -------------------------------------------------

#omega0:=12; # Input

#NumFrames := 50;

ForceScale:=0.15:

#NumPtsPath:=180:

xx0:=array(1..NumPtsPath):yy0:=array(1..NumPtsPath):

for i from 0 to NumPtsPath-1 do; alpha:=i/NumPtsPath; setup(alpha*180); xx0[i+1]:=x(0,omega0);yy0[i+1]:=y(0,omega0);od:

TipPath:=pointplot( {seq( [xx0[i],yy0[i]],i=1..NumPtsPath)}, color=red , symbolsize=10):

NumPts:=50:

BaseCurve:=array(1..NumFrames):ForceVect:=array(1..NumFrames):

for i from 0 to NumFrames - 1 do alpha:=i/NumPts; setup(alpha*180);beta0:=beta(omega0);BaseCurve[i+1]:= plot([x(s,omega0),y(s,omega0), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=3,caption=typeset("frame #",i,"/",NumFrames," omega = ", omega0," alpha = ",alpha*180.,"  beta = ",evalf(beta0*180/Pi))); ForceVect[i+1]:= arrow(<x(0,omega0),y(0,omega0)>,force(ForceScale,beta0)):od:

Beam := display(seq(BaseCurve[i],i=1..NumFrames),insequence = true, scaling=constrained):

Force :=display(seq(ForceVect[i],i=1..NumFrames),insequence = true, scaling=constrained):

display( Beam, Force,TipPath);

#========================================

#  Load parameter problem

#========================================

# These functions calculate load parameter as function of alpha and wave number n

#-------------------------------------------------------------------------

f1:=(alpha,n)->(-InverseJacobiSN(sin(beta0*Pi/360)/sin(alpha*Pi/2),sin(alpha*Pi/2)))+(4*n+1)*EllipticK(sin(alpha*Pi/2)):

f2:=(alpha,n)->-(-InverseJacobiSN(sin(beta0*Pi/360)/sin(alpha*Pi/2),sin(alpha*Pi/2)))+(4*n+3)*EllipticK(sin(alpha*Pi/2)):

#---------------------------------------------------------

beta0:=30.; # Input

if beta0 < 0 then beta0 := -beta0 end if:

m:=2; # Input wave number

if beta0 < 0 then beta0 := -beta0 end if;p:=evalf(sin(rad(beta0)/2)): # Check data

ma:=abs(m):if modp(abs(m),2) = 1 then n:=floor((ma-1)/2);else n:=floor(ma/2-1);end if:da:=(0.999-beta0/180)/(NumFrames-1):

for i from 0 to NumFrames - 1 do if m > 0 then alpha:=evalf(beta0/180+da*i) else alpha := evalf(-beta0/180 - da*i); end if; setup(evalf(alpha*180));omega0:=Re(-InverseJacobiSN(p/k,k));if modp(n,2) = 1 then w[i+1]:=omega0+(4*n+1)*K else  w[i+1]:=-omega0 + (4*n + 3)*K end if; aa[i+1]:=alpha;od:

ss:=seq([aa[i],w[i]],i=1..NumFrames):

ww:=pointplot( {ss},color=black, symbol=circle,symbolsize = 10 ):

mx:=w[NumFrames]:

unassign('alpha');tt:=plot({seq(f1(alpha,n),n=0..2),seq(f2(alpha,n),n=0..2)},alpha=-1..1,y=0..mx,color=red,labels=["alpha/pi","omega"]):

display(tt,ww);

ForceScale:=0.2: # force scale factor

NumPts:=50:

for i from 0 to NumFrames - 1 do alpha:=aa[i+1]; setup(evalf(alpha*180));omega:=w[i+1]; BaseCurve[i+1]:= plot([x(s,omega),y(s,omega), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=3,caption=typeset("frame #",i+1,"/",NumFrames," m = ",m," omega = ", omega," alpha = ",alpha*180.,"  beta = ",evalf(beta0))); ForceVect[i+1]:= arrow(<x(0,omega),y(0,omega)>,force(ForceScale*omega/w[NumFrames],evalf(beta0*Pi/180))):od:

Beam := display(seq(BaseCurve[i],i=1..NumFrames),insequence = true, scaling=constrained):

Force :=display(seq(ForceVect[i],i=1..NumFrames),insequence = true, scaling=constrained):

display( Beam, Force);

beta0:=40; # Input

alpha0:=145; # Input alpha

NumWaves:=6; # How many waves

ForceScale:=0.15: # force scale factor

if beta0 < 0 then beta0 := -beta0 end if:p:=evalf(sin(beta0*Pi/180/2)):setup(alpha0): # check beta

for i from 1 to NumWaves do omega0:=Re(K-InverseJacobiSN(p/k,k));if modp(abs(i),2) = 1 then n:=(abs(i)-1)/2;  w[i]:=omega0+4*n*K else n:=abs(i)/2-1; w[i]:=-omega0 + 4*(n + 1)*K end if;od:

ss:=seq([alpha0/180.,w[i]],i=1..NumWaves):

ww:=pointplot( {ss},color=black, symbol=solidcircle,symbolsize = 10,labels=["alpha/pi","omega"] ):

unassign('alpha');nn:=floor(NumWaves/2)+1:graph:=plot({seq(f1(alpha,n),n=0..nn),seq(f2(alpha,n),n=0..nn)},alpha=-1..1,y=0..w[NumWaves],color=red,labels=["alpha/pi","omega"]):

alp:=pointplot( [[alpha0/180,0],[alpha0/180,w[NumWaves]]],connect=true,color=black,thickness=1,labels=["alpha/pi","omega"]):

display(graph,ww, alp);

for i from 1 to NumWaves do omega:=w[i]; BaseCurve[i]:= plot([x(s,omega),y(s,omega), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=2,caption=typeset("wave = ",i," omega = ", omega," alpha = ",alpha0,"  beta = ",evalf(beta0))); ForceVect[i]:= arrow(<x(0,omega),y(0,omega)>,force(ForceScale*omega/w[NumWaves],evalf(beta0*Pi/180))):od:

Beam := display(seq(BaseCurve[i],i=1..NumWaves),insequence = true, scaling=constrained):

Force :=display(seq(ForceVect[i],i=1..NumWaves),insequence = true, scaling=constrained):

display( Beam, Force);

Beam := display(seq(BaseCurve[i],i=1..NumWaves), scaling=constrained):

Force :=display(seq(ForceVect[i],i=1..NumWaves), scaling=constrained):

display( Beam, Force);

#========================================

#  (3) Conservative Load Problem

#========================================

gc(): # cleanup

omega0:=15; # Input

#---------------------------------------------

# This is variant of procedure beta where also k and K are calculates

f:=proc(alpha,omega) global k, K; k:=evalf(sin(alpha/2));K:=EllipticK(k); beta(omega); end proc:

#------------------------------------------

# This is first derivative of function f which is needed to calculate extrem values

df := (k,omega)->2*(JacobiSN(omega,k)*omega*EllipticK(k)*k^2-JacobiSN(omega,k)*omega*EllipticK(k)+JacobiCN(omega,k)*JacobiDN(omega,k)*EllipticK(k)+JacobiSN(omega,k)*JacobiZeta(omega,k)*EllipticK(k)+JacobiSN(omega,k)*omega*EllipticE(k))/JacobiDN(omega,k)^2/EllipticK(k)/((-1+k^2)/(-1+k^2*JacobiSN(omega,k)^2))^(1/2):

#---------------------------------------------

# These two functions enclose n-th root of function f

f1:=n->2*arcsin(sqrt(1-exp(Pi-omega0*2/(2*n-1)))):f2:=n->2*arcsin(sqrt(1-16*exp(-omega0*2/(2*n-1)))):

#----------------------------------

# Calculate roots of f in the interval -Pi..Pi

#--------------------------------------------

dg:=Digits: Digits:=22;   # Increase number of digits

N:=floor(omega0/Pi+1/2):  # half number of possible nontrivial eqilibrium configurations

M:=2*N+1;                 # max. number of equilibrium configurations

z0:=array(1..M):          # initialize array

unassign('alpha'):for n from 1 to N do alpha1:=evalf(f1(n)); alpha2:=evalf(f2(n));z0[n]:=evalf(fsolve(f(alpha,omega0),alpha=alpha1..alpha2)/Pi);od:

# obtain rest of roots by symetry

z0[N+1]:=0:for n from 1 to N do z0[N+1+n]:=-z0[N+1-n] od:

#------------------------

# Form graph

#---------------------

zero0:=(seq( pointplot( [z0[i],0], symbol = circle, symbolsize = 10),i=1..M)):

graph:=plot({evalf(f(alpha*Pi,omega0)/Pi)},alpha=-1..1, color=black, numpoints=400): # note that number of points should be relatively large

#display(graph,zero0);

#-------------------------------------

# Calculate extems of function f

#------------------------------------

s:=array(1..M+1):r:=array(1..M+1):

s[1]:=1:r[1]:=evalf(Pi):

for n from 1 to M-1 do;if modp(n,2) =1 then xx:=fsolve(df(evalf(sin(alpha*Pi/2)),omega0),alpha=z0[n+1]..z0[n]); else xx:=fsolve(df(evalf(sin(alpha*Pi/2)),omega0),alpha=z0[n+1]..z0[n]); end if; s[n+1]:=xx;r[n+1]:=f(xx*Pi,omega0) od:s[M+1]:=-1:r[M+1]:=evalf(-Pi):

#for i from 1 to M do z0[i],r[i]; od;

ext0:=(seq( pointplot( [s[i],evalf(r[i]/Pi)], symbol = circle, symbolsize = 10, color=blue),i=1..M+1)):

#display(graph,zero0,ext0);

beta0:=0; # Input <---------

#unassign('alpha'):fsolve(f(evalf(alpha*Pi),omega0)-rad(beta0),alpha=z0[1]..r[1]);z0[1];r[1];

if beta0 < 0 then beta0 := -beta0 end if; # Check

#-------------------------------------- Here the calulation phase is done !

if (beta0 > 0) then n:=0;for i from 1 to M do zz:=fsolve(f(evalf(alpha*Pi),omega0)-rad(beta0),alpha=z0[i]..s[i]): if type(zz,numeric) then n:=n+1;zzz[n]:=zz end if; zz:=fsolve(f(evalf(alpha*Pi),omega0)-rad(beta0),alpha=s[i+1]..z0[i]): if type(zz,numeric) then n:=n+1;zzz[n]:=zz end if; od: else n:=M; for i from 1 to M do zzz[i]:=z0[i] od; end if:

NumShape:=n;

for i from 1 to NumShape do i,180*zzz[i] od;

zero:=(seq( pointplot( [zzz[i],evalf(beta0/180)], symbol = solidcircle, symbolsize = 12, color=red),i=1..n)):

display(graph,zero0,ext0,zero,plot(beta0/180,t=-1..1),labels=["alpha/pi","beta/Pi"]);

NumPts:=50:for i from 1 to NumShape do alpha:=zzz[i]; setup(alpha*180); BaseCurve[i]:= plot([x(s,omega0),y(s,omega0), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=2,caption=typeset("shape = ",i,"/",NumShape," omega = ", omega0," alpha = ",alpha*180.,"  beta = ",evalf(beta0))); ForceVect[i]:= arrow(<x(0,omega0),y(0,omega0)>,force(ForceScale,evalf(beta0*Pi/180))):od:

Beam := display(seq(BaseCurve[i],i=1..NumShape), scaling=constrained, insequence = true):

Force :=display(seq(ForceVect[i],i=1..NumShape), scaling=constrained, insequence = true):

display( Beam, Force);

Beam := display(seq(BaseCurve[i],i=1..NumShape), scaling=constrained):

#Force :=display(seq(ForceVect[i],i=1..NumShape), scaling=constrained):

display( Beam);

Digits:=dg: # reset number of digits

# END