230.sagews

var('R')
d2 = 2.265*2.
d1 = 1.
d12  = (d2 + d1)/2
beta = 0.85
du = 3.015
rho = 0.2
@interact
def _(u=(0.001, 2, 0.1)):
  f1 = exp(beta*u)*(-((R^2 + d12^2 -du^2)/(2*R))^2 + d12^2)
  f2 = exp(2*beta*u)*(-(R/2)^2 + ((R^2 + d12^2 -du^2)/(2*R))^2)
  f3 = (-(R/2)^2 + du^2)
  f4 = exp(beta*u)*(-((R^2 + du^2 - d12^2)/(2*R))^2 + (R/2)^2)
  f = - pi*rho/beta*(f1 + f2 + (1 - exp(beta*u))*(f3 + f4))
  g2 = exp(2*beta*u)*( - d12^2 + ((R^2 + d12^2 -du^2)/(2*R))^2)
  g = - pi*rho/beta*(f1 + g2 + (1 - exp(beta*u))*(f3 + f4))
  h4 = exp(beta*u)*(-du^2 + (R/2)^2)
  h = - pi*rho/beta*(1 - exp(beta*u))*(f3 + h4)
  j = 0
  pt1 = plot(f, (d2, d12*2), color='red', thickness=1)
  dot1 = point((2*d12,f(2*d12)),pointsize=10,rgbcolor=(0,0,0))
  pt2 = plot(g, (d12*2, d12 + du), color='red', thickness=1)
  dot2 = point((du + d12, g(du + d12)),pointsize=10,rgbcolor=(0,0,0))
  pt3 = plot(h, (d12 + du, 2*du), color='red', thickness=1)
  dot3 = point((2*du, h(2*du)),pointsize=10,rgbcolor=(0,0,0))
  pt4 = plot(j, (2*du, 2.5*du), color='red', thickness=1)
  show(pt1 + pt2 + pt3 + pt4 + dot1 + dot2 + dot3, xmin = d2, xmax = 2.5*du)