# S.A. Hires, L. Pammer, K. Svoboda and D. Golomb
# Tapered whiskers are required for active tactile sensation.
# eLIFE 2:e01350, 2013.
#
# Quasi-static solution of the bending of an isolated whisker.
# File for bifurcation calculation.
#
# For Fig. 2, stable solution, use: theta=-0.174533 rad.
# Use "Bndryval -> Show" in XPPAUT.
#
# p xcen=15.1286239388569 ycen=4.28562676857186 
# p xcen=14.8713760611431 ycen=4.71437323142814
p thetap=0.0
p dd=15.6604597633658
# CB: positive rpole, negative theta
p rpole=0.25
# CF: negative rpole, positive theta
##p rpole=-0.25
p g0=0.0 g1=0.0 g2=0.02 g3=0.0 g4=0.0
#
p EE=1.0,rbase=0.03,rtip=0.0015
p LL=20.0
#
gg(vv)=((((g4*vv)+g3)*vv+g2)*vv+g1)*vv+g0
ggp(vv)=(((4.0*g4*vv)+3.0*g3)*vv+2.0*g2)*vv+g1
ggpp(vv)=((12.0*g4*vv)+6.0*g3)*vv+2.0*g2
sqgpxy(vv)=sqrt(1.0+(ggp(vv))^2.0)
#
rr=(1.0-(sobj*sighat/LL)*(1.0-(rtip/rbase)))*rbase
IIarea=0.25*Pi*(rr^4.0)
coefaux=FF/(EE*IIarea)
#
xxaux'=xlen
yyaux'=xlen*(2.0*g2*xxaux)
xlen'=0
xcen'=0
ycen'=0
#
xu'=sobj/sqgpxy(xu)
xx'= sobj*cos(theta)
yy'= sobj*sin(theta)
theta'=(sobj*ggpp(xu)/(sqgpxy(xu)^3))+sobj*coefaux*((xcen-xx-rpole*sin(thobj))*cos(thobj)+(ycen-yy+rpole*cos(thobj))*sin(thobj))
FF'=0
thobj'=0
sobj'=0.0
sighat'=1
#
boundary xxaux
boundary yyaux
boundry xcen'-(xxaux'+rpole*2.0*g2*xxaux'*sqrt(1.0/(1+(4.0*g2*g2*xxaux'*xxaux'))))
boundary ycen'-(yyaux'-rpole*sqrt(1.0/(1+(4.0*g2*g2*xxaux'*xxaux'))))
boundary xxaux'*xxaux'+yyaux'*yyaux'-dd*dd
#
boundary xu
boundary xx
boundary yy
boundary theta-thetap
boundary theta'-thobj'
boundary xx'-xcen+rpole*sin(thobj')
boundary yy'-ycen-rpole*cos(thobj')
boundary sighat
#
xxaux(0)=0.0
yyaux(0)=0.0
xlen(0)=15.0
xcen(0)=15.0
ycen(0)=4.5
#
xu(0)=0.0
xx(0)=0.0
yy(0)=0.0
theta(0)=thetap
sobj(0)=18.0
sighat(0)=0.0
FF(0)=0.0e-8
thobj(0)=0.54042
#
aux kapud=sobj*ggpp(xu)/(sqgpxy(xu)^3)
aux kapd=sobj*coefaux*((xcen-xx-rpole*sin(thobj))*cos(thobj)+(ycen-yy+rpole*cos(thobj))*sin(thobj))
#
@ total=1
@ DT=0.001
@ MAXSTOR=800000
@ BACK=Black
@ XP=xx
@ YP=yy
@ XLO=0.0, XHI=16.0, YLO=-1.0, YHI=5.0
@ NTST=50,NMAX=2000,NPR=50
@ DS=0.002,DSMIN=0.0001,DSMAX=0.05
@ PARMIN=-5.0,PARMAX=5.0,NORMMIN=0.0,NORMMAX=10000.0
@ AUTOVAR=ycen,AUTOXMIN=0.0,AUTOXMAX=5.0,AUTOYMIN=0.0,AUTOYMAX=4.0
done