# colors
cc1 = "#5e81b5"
cc2 = "#e19c24"
cc2l = "#f9ebd3"
cc3 = "#648FFF"
cc4 = "#DC267F"
cc5 = "#000000"

Tc = 17.26

a1 = 71.13597475161484
b1 = 0.010425748328804992
c1 = 1665.064389018724
d1 = 6.283065722894796
e1 = 14.579664883988315

C10(T) = a1 - b1 * T
C1(T) = C10(T) / (1 + 1 / (1 / d1 + abs(T - Tc) / c1) / C10(T))
C11(T) = C10(T) \
  + e1**2 / (2 * C10(T)) \
  + e1**2 * ((T - Tc) / c1) / 12 \
  - e1**2 * (9 * C10(T) + d1) / (6 * C10(T) * (9 * C10(T) + d1 + C10(T) * d1 * (T - Tc) / c1)) \
  + d1 * (C10(T) + d1) * (12 * C10(T)**2 + 5 * e1**2) / (12 * C10(T) * ((C10(T) + d1) + C10(T) * d1 * (T - Tc) / c1)**2) \
  - (24 * C10(T)**2 * d1 + (4 * C10(T) + 9 * d1) * e1**2) / (12 * C10(T) * (C10(T) + d1 + C10(T) * d1 * (T - Tc) / c1))
C12(T) = T < Tc ? C11(T) : C1(T)

a2 = 144.345
b2 = 0.019492
c2 = 120.462 

C20(T) = a2 - b2 * T**2 / (c2 + T)

set terminal epslatex size 18.00cm, 11.12cm standalone header \
  '\usepackage{xcolor}\definecolor{mathc3}{HTML}{1a85ff}\definecolor{mathc4}{HTML}{d41159}'
set output "fig-stiffnesses.tex"

set multiplot  layout 2,2  margins 0.075, 0.915, 0.1, 0.99  spacing 0.01, 0.01

set nokey


set mxtics 2
set mx2tics 2
set mytics 2
set my2tics 2
set format y '\scriptsize $%1.f$'
set format y2 '\scriptsize $%.1f$'
set format x '\scriptsize $%.0f$'
#set format x2 '\scripsize $%.0f$'
set format x2 ''

set style rect back fc rgb cc2l
set object 1 rect from 0,0 to Tc,200
set border lw 2 front
set tics front

set arrow 1 from Tc,graph 0 to Tc,graph 1 nohead lw 5 lc rgb cc2
unset xtics
set x2tics 50, 50, 250 offset 0,-0.5 mirror
#set x2label '\tiny $T / \mathrm K$' offset 0,-1.0
set xrange [0:300]

set title '(a)' offset 12,-2.7
set ylabel '$C_{\mathrm{B_{2\mathrm g}}}$ (GPa)' offset 3
set yrange [140.2:144.7]
set ytics 139,1,146 offset 0.5 mirror
plot "data/c66.dat" using 1:(100 * $2) with lines lw 6 lc rgb cc3, \
     C20(x) dt 3 lw 4 lc rgb cc4

set title '(b)' offset 12,-2.7
unset ylabel
set y2label '$C_{\mathrm{B_{1\mathrm g}}}$ (GPa)' offset -4 rotate by -90
set yrange [65.05:65.68]
set y2tics 62.1,0.1,66 offset -0.5 mirror
set samples 1000
plot "data/c11mc12.dat" using 1:(100 * $2) with points pt 6 ps 0.4 lc rgb cc3, \
     x > Tc ? C12(x) : 1 / 0 lw 3 lc rgb cc5, \
     x <= Tc ? C12(x) : 1 / 0 dt (2,1) lw 3 lc rgb cc5, \
     C1(x)  dt (1,2) lw 3 lc rgb cc4



unset x2label
unset x2tics

set title '(c)' offset 12,-2.7
set xlabel '$T$ (K)\\~' offset 0,0.5
set xtics 50, 50, 250 offset 0,0.5 mirror
unset y2tics
set ytics 63,1,74
unset y2label
set ylabel '$\textcolor{mathc3}{C_{\mathrm{B_{1\mathrm g}}}},\textcolor{mathc4}{C_{\mathrm{B_{1\mathrm g}}}^0}$ (GPa)' offset 3

set yrange [64.5:71.5]
plot "data/c11mc12.dat" using 1:(100 * $2) with lines lw 6 lc rgb cc3, \
     C10(x) dt 3 lw 4 lc rgb cc4

set ylabel ''
set y2label '$\Big[C^0_{\mathrm{B_{1\mathrm g}}}(C^0_{\mathrm{B_{1\mathrm g}}}/C_{\mathrm{B_{1\mathrm g}}} - 1)\Big]^{-1}$ (GPa$^{-1}$)' offset -3 rotate by -90

set title '(d)' offset 12,-2.7
set format y2 '\scriptsize$%0.2f$'
set format y '\scriptsize$%0.2f$'
set yrange [0.14:0.38]
set y2tics 0.05,0.05,0.5 offset -0.7 mirror
plot "data/c11mc12.dat" using 1:(1 / (C10($1)*(C10($1) / (100 * $2) - 1))) with lines lw 6 lc rgb cc3, \
     x > Tc ? 1/(C10(x) * (C10(x) / C12(x) - 1)) : 1 / 0 lw 3 lc cc5, \
     x <= Tc ? 1/(C10(x) * (C10(x) / C12(x) - 1)) : 1 / 0 dt (2,1) lw 3 lc cc5, \
     1/(C10(x) * (C10(x) / C1(x) - 1)) dt (1,2) lw 3 lc rgb cc4

set size 0.37,0.3
set origin 0.55,0.52
unset margins
set xtics 0,5,35
set y2tics 65.1,0.03,65.3
set xrange [7:27]
set yrange [65.14:65.26]
unset title
unset xlabel
unset y2label

plot "data/c11mc12.dat" using 1:(100 * $2) with points pt 2 ps 0.15 lc rgb cc3, \
     x > Tc ? C12(x) : 1 / 0 lw 3 lc rgb cc5, \
     x <= Tc ? C12(x) : 1 / 0 dt (2,1) lw 3 lc rgb cc5, \
     C1(x)  dt (1,2) lw 3 lc rgb cc4

unset margins
set origin 0.5,0.26
set yrange [0.154:0.166]
set y2tics 0.155,0.005,0.167
set format y2 '\scriptsize$%0.3f$'

plot "data/c11mc12.dat" using 1:(1 / (C10($1)*(C10($1) / (100 * $2) - 1))) with points pt 2 ps 0.15 lc rgb cc3, \
     x > Tc ? 1/(C10(x) * (C10(x) / C12(x) - 1)) : 1 / 0 lw 3 lc cc5, \
     x <= Tc ? 1/(C10(x) * (C10(x) / C12(x) - 1)) : 1 / 0 dt (2,1) lw 3 lc cc5, \
     1/(C10(x) * (C10(x) / C1(x) - 1)) dt (1,2) lw 3 lc rgb cc4

unset multiplot