#! /bin/bash
# Last edited on 2013-01-24 23:46:44 by stolfilocal

# Plots Hooke's law for a simple spring


prefix="out/Hookes_spring_law"

tmp=/tmp/$$

# Create a file with selected abscissas:
XD="-0.0250:-0.0125:00.0000:+0.0125:+0.0250:"
XA=( ${XD//:/ } )
echo "${XD}" | tr ':' '\012' > ${tmp}.dat
echo "${XA[0]} foo ${XA[1]} "
  
export GDFONTPATH=.

gnuplot <<EOF
xsize = 540   # Nominal canvas X size
ysize = 720   # Nominal canvas Y size
bfrac = 0.40  # Fraction of canvas reserved for spring drawings

set term svg size (xsize),(ysize) font "arialbd,16" enhanced rounded dashed
set output "${prefix}.svg" 
set nokey
set style line 10 lt 1 lw 2 lc rgb '#888888'
set border 0 ls 10
set samples 2000
set lmargin at screen 0.05
set rmargin at screen 0.95
PI = 3.14926

# Spring profile at height y=0, x spanning [0 _ L], with radius R and N turns:
SprT(N,L) = L/(N+0.5)          # Wavelength.
SprB(N,L) = 0.25*SprT(N,L)     # Length of squashed part.
SprP(N,L,dx) = dx*(5*(dx/SprB(N,L))**2 + 1) # Reparam in squashed part (dx = dist from midpart).
SprPD(N,L,dx) = (15*(dx/SprB(N,L))**2 + 1)  # Derivative of reparam in squashed part.

SprF(N,L,x) = \
  ( (x < 0) || (x > L) ? 0/0 : \
    ( x < SprB(N,L) ? -SprP(N,L,SprB(N,L)-x) : \
      ( x > L-SprB(N,L) ? N*SprT(N,L) + SprP(N,L,x-(L-SprB(N,L))) : \
        x - SprB(N,L) \
      ) \
    ) \
  ) # Reparametrization that squashes the ends of the spring.

SprRM(L0,R0,N,L) = sqrt((2*PI*R0)**2 + (SprT(N,L0))**2 - (SprT(N,L))**2)/(2*PI) # Radius of spring in middle part.
SprRP(L0,R0,N,L,dx) = SprRM(L0,R0,N,L) + (R0-SprRM(L0,R0,N,L))*(dx/SprB(N,L))**2 # Radius of spring in squashed part.

SprR(L0,R0,N,L,x) = \
  ( (x < 0) || (x > L) ? 0/0 : \
    ( x < SprB(N,L) ? SprRP(L0,R0,N,L,SprB(N,L)-x) : \
      ( x > L-SprB(N,L) ? SprRP(L0,R0,N,L,x-(L-SprB(N,L))) : \
        SprRM(L0,R0,N,L) \
      ) \
    ) \
  ) # Radius of spring at x 


Spring(L0,R0,N,L,x) = SprR(L0,R0,N,L,x)*sin(2*PI*SprF(N,L,x)/SprT(N,L))

# Spook
S(B,x) = (x <= 0 ? 0 : ( x >= 1 ? 0/0 : exp(B*(1.0/(1-x) - 1.0/x)) ) )

# Vertical line at X
V(X,x) = (abs(x-X) > 0.002 ? 0/0 : 100000*(x-X) )

# Vertical half-line with upper endpoint at X,Y
U(X,Y,x) = (abs(x-X) > 0.002 ? 0/0 : Y-100000*abs(x-X) )

W = 0.048; # Plot width

# Critical displacements (m):
X0 = -0.0420 # Fixed point of spring
X1 = -0.0360 # Min elastic displ (m)
X2 = +0.0470 # Max elastic displ (m)

# End-of-expansion term:
H(M,B,x) = S(B,x/M)

# End-of-compression term:
G(M,B,x) = - S(B,x/M)

# Cuboid term
Z(L,x) = (x/L)**3

# Non-linear term for positive x:

# Total nonlinear term:
GA = 0.005; GB = 1;
HA = 0.13; HB = 1;
ZA = 0.30; ZL = X2;
FN(x) = ZA*Z(ZL,x) + GA*G(X1,GB,x) + HA*H(X2,HB,x)
# FN(x) = GA*G(X1,GB,x) + HA*H(X2,HB,x)

# Linear term (Hooke's law):
K = 100.0 # Hooke's constant, N/m
FL(x) = K*x

# Hooke's law plus non-linear term:
FT(x) = FL(x) + FN(x)
 
Fmin = 1.8*FL(X1)
Fmax = 1.8*FL(X2)
 
set xrange [-W:+W]

set multiplot layout 2,1
# ----------------------------------------------------------------------
set bmargin at screen 0.0
set tmargin at screen (bfrac)

#Set Y range so that scale is 1:1
aspect = (ysize*bfrac)/xsize
yspan = aspect*(2*W)
set yrange [0:yspan]

dy = yspan/12.0 # Twice the spacing between springs

unset ytics
unset xtics
unset yzeroaxis
unset xzeroaxis
set clip two

# Spring parameters:
SL0 = (${XA[2]}-X0) # Relaxed length
SR0 = 0.60*dy       # Radius of relaxed spring
SN = 6              # Number of turns
ST0 = SprT(SN,SL0)  # Wavelength of relaxed spring
SLM = 0.99*(SN + 0.5)*sqrt((2*PI*SR0)**2 + (ST0)**2)

print "ST0 = ", ST0
print "SL0 = ", SL0
print "SLM = ", SLM

# Force arrows:
fracF = 0.85
ArLen(X) = ( X = 0 ? 0/0 : 0.5*FT(X)/K )
ArBasPos(X) = ( X <= 0 ? 0/0 : X + 0.20*dy )
ArBasNeg(X) = ( X >= 0 ? 0/0 : X - ArLen(X) + 0.20*dy )

set style line 2 lw 3 lt 1 lc rgb '#ff2200'
set style line 3 lw 3 lt 1 lc rgb '#0066ff'
set style line 4 lw 2 lt 1 lc rgb '#888888'

set style arrow 2 size char 1.5,30 ls 2
set style arrow 3 size char 1.5,30 ls 3
set style arrow 4 size char 1.5,30 ls 4

plot \
  "${tmp}.dat" using 1:(dy*(10-2*column(0))):(0):(+100*dy) with vectors lw 2 lt 3 lc rgb '#888888', \
  \
  (Spring(SL0,SR0,SN,SL0, x-X0) + 9*dy) with lines lt 1 lw 8 lc rgb '#cccccc', \
  (Spring(SL0,SR0,SN,SL0, x-X0) + 7*dy) with lines lt 1 lw 8 lc rgb '#cccccc', \
  (Spring(SL0,SR0,SN,SL0, x-X0) + 5*dy) with lines lt 1 lw 8 lc rgb '#cccccc', \
  (Spring(SL0,SR0,SN,SL0, x-X0) + 3*dy) with lines lt 1 lw 8 lc rgb '#cccccc', \
  (Spring(SL0,SR0,SN,SL0, x-X0) + 1*dy) with lines lt 1 lw 8 lc rgb '#cccccc', \
  \
  "${tmp}.dat" using (ArBasPos(column(1))):(dy*(10-2*column(0)-1)):(ArLen(column(1))):(0) \
    with vectors as 2, \
  "${tmp}.dat" using (ArBasNeg(column(1))):(dy*(10-2*column(0)-1)):(ArLen(column(1))):(0) \
    with vectors as 3, \
  \
  (Spring(SL0,SR0,SN,(${XA[0]}-X0), x-X0) + 9*dy) with lines lt 1 lw 8 lc rgb '#000000', \
  (Spring(SL0,SR0,SN,(${XA[1]}-X0), x-X0) + 7*dy) with lines lt 1 lw 8 lc rgb '#000000', \
  (Spring(SL0,SR0,SN,(${XA[2]}-X0), x-X0) + 5*dy) with lines lt 1 lw 8 lc rgb '#000000', \
  (Spring(SL0,SR0,SN,(${XA[3]}-X0), x-X0) + 3*dy) with lines lt 1 lw 8 lc rgb '#000000', \
  (Spring(SL0,SR0,SN,(${XA[4]}-X0), x-X0) + 1*dy) with lines lt 1 lw 8 lc rgb '#000000'


#   (SR0 + 11*dy) with lines lw 1 lt 1 lc rgb '#888800', \
#   (SprR(SL0,SR0,SN,SLM, x-X0) + 11*dy) with lines lw 1 lt 1 lc rgb '#880088', \
#   \
#   (Spring(SL0,SR0,SN,SL0, x-X0) + 11*dy) with lines lt 1 lw 8 lc rgb '#cccccc', \
#   (Spring(SL0,SR0,SN,SLM, x-X0) + 11*dy) with lines lt 1 lw 8 lc rgb '#000000', \
# 

# ----------------------------------------------------------------------
set bmargin at screen (bfrac)
set tmargin at screen 0.95
set yrange [(Fmin):(Fmax)]

ExLabX = 0.027
ExLabY = -0.65
ExLabel = "Elongation"
set obj rect at first (ExLabX-0.001),(ExLabY+0.01) size char strlen(ExLabel)+1,1.1 lw 0 front fc rgb '#ffffff'
set label (ExLabel) at first (ExLabX),(ExLabY) center front

CmLabX = -0.027
CmLabY = ExLabY
CmLabel = "Compression"
set obj rect at first (CmLabX-0.001),(CmLabY+0.01) size char strlen(CmLabel)+1,1.1 lw 0 front fc rgb '#ffffff'
set label (CmLabel) at first (CmLabX),(CmLabY) center front

set label "Force"  at first -0.0025,2.8 rotate by 90

set label "{/arialbi X}" at first (0.8*W),+0.48 center front enhanced
set label "{/arialbi F}" at first 0.002,(0.8*Fmax) enhanced

unset ytics
unset xtics
unset yzeroaxis
unset xzeroaxis

set arrow from first 0,0 to first 0,0.98*Fmax    as 4
set arrow from first 0,0 to first -0.98*W        as 4
set arrow from first 0,0 to first +0.98*W        as 4

#set xtics out scale 2 ( '' (${XD//:/), '' (}100) ) nomirror

plot \
  "${tmp}.dat" using 1:(FT(column(1))):(0):(-1000) with vectors lw 2 lt 3 lc rgb '#888888', \
  (FT(x)) with lines lw 2 lt 2 lc rgb '#774433', \
  "${tmp}.dat" using 1:(FT(column(1))) with points lw 3 pt 7 ps 0.45 lc rgb '#553322', \
  (x >= 0 ? FL(x) : 0/0) with lines ls 2, \
  (x <= 0 ? FL(x) : 0/0) with lines ls 3

quit

EOF