
Вараны острова Комодо, также называемые в литературе драконами, — самая крупная из живущих на земле ящериц. Длина его тела может достигать 3 метров, а масса 140 кг [1]. Это доминирующий хищник своего региона, который может добывать животных (свиньи, буйволы, олени), порой 10-ти кратно превосходящих его весу.
Важнейшим инструментом такой охотничьей эффективности являются зубы. У комодского варана их 60 штук [2], изогнутых как сабли и острых как бритва (край зуба усилен металлизированным слоем, образующим микро пилу [3]).
Этот комплект еще и регулярно, раз в 40 дней обновляется [4]. Не нужно ни стоматологов ни заточников — просто мечта. Однако фантастическая скорость роста зубов должна требовать и фантастических затрат «стройматериалов». Сколько, например, кальция и железа нужно варану в день для поддержания такого темпа?
Ниже мы оценим эти показатели, опираясь на «ангем», «матан» и python. Кто не испугался, welcome.
Содержание
Введение
Оценивать в миллиграммах содержание в ткани того или иного микроэлемента будем по формуле  ,  где  
 — плотность, V — объем.
Чтобы рассчитать объем, нужно знать форму объекта. А лучше — уравнение поверхности, его ограничивающей. С этого и начнем. Вообще задача моделирования поверхности интересна сама по себе, поэтому решим ее в нескольких вариантах, переходя от более простого к более точному.
Геометрия зуба
В качестве экспериментальных данных возьмем фотографии и описание из недавнего исследования, опубликованного в журнале Nature Ecology & Evolution [3]:

На изображении видно, что зубы варана имеют загнутую и слегка сплющенную с боков форму. Длина зубов 2.5–3 см (что также указано в [1]). Зафиксируем для дальнейших расчетов высоту зуба равной 25 мм, глубину 10 мм, ширину 4 мм.
Зуб варана v.1 (конус)
Объем конического зуба
В качестве первого приближения возьмем конус с эллиптическим основанием.
Объем конуса, как известно со школы,   
где S — площадь основания, h — высота.
В нашем случае h = 25 мм
В случае эллипса  ,
где a и b — большая и малая полуоси эллипса, лежащего в основании.
В нашем случае a = 5 мм, b = 2 мм.
Таким образом,   
Подставим численные значения:   
Итак объем зуба V = 261.7 мм3. Это наша первая грубая оценка.
Уравнение поверхности и визуализация
Способ 1 (быстрый)
Изобразим конус в логике аналитической геометрии, т.е. опираясь на уравнение его поверхности.
Уравнение конуса запишется как  
Или, после подстановки значений a, b и h:  
Чтобы задать конечную высоту зубу (сам по себе конус фигура бесконечная), наложим ограничение на диапазон изменения z:
  
Для быстрой визуализации можем воспользоваться, например, сервисом desmos.com. Вводим полученное уравнение как есть и наблюдаем изображение конуса.

Способ 2 (с полным контролем)
Изобразим конус при помощи Python и Matplotlib.
В данном случае будет удобнее взять уравнение конуса в параметрическом виде:
 
где параметры изменяются в следующих диапазонах:  , 
 
Код:
import numpy as np 
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D 
# Параметры конуса 
a, b, h = 2, 5, 25
# Координаты точек на поверхности конуса 
phi = np.linspace(0, 2 * np.pi, 30) 
tz = np.linspace(-h, 0, 20) 
phi, tz = np.meshgrid(phi, tz)
# Параметрическое уравнение конуса 
z = tz 
x = (a / h) * tz * np.cos(phi) 
y = (b / h) * tz * np.sin(phi)
# Создание графика 
ax = plt.figure().add_subplot(111, projection='3d') 
ax.view_init(elev=20, azim=-40, roll=0)
ax.plot_surface(x, y, z, color='w', edgecolor='royalblue',             	rstride=2, cstride=2, alpha=0.5) ax.set_aspect('equal')
# Диапазоны и настройки графика 
x_lims = [-2a, 2a] 
y_lims = [-1.5b, 1.5b] 
z_lims = [-h, 0]
ax.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
       xlabel='X', ylabel='Y', zlabel='Z',
       xticks=np.round(np.linspace(*x_lims, 5)))
plt.show()В результате получим следующий график:

Однако поверхность зуба все же изогнутая, а не линейная как у конуса, поэтому стоит добавить кривизны. Понизим степень в третьем слагаемом и получим уравнение параболоида.
Зуб варана v.2 (параболоид)
На сей раз начнем с уравнения поверхности.
Уравнение поверхности и визуализация
Уравнение параболоида имеет вид    
При подстановке заданных параметров получаем следующее уравнение боковой поверхности зуба:
  
Как и прежде, ограничим диапазон изменения z:  
Вводим полученное уравнение в сервисе desmos.com или аналогичном и наблюдаем соответствующий график параболоида.

Объем параболического зуба
Найдем объем зуба параболической формы как тройной интеграл по области пространства, ограниченной поверхностью параболоида и плоскостью .  
В общем случае в выражении через повторные интегралы это выглядит так
  
Чтобы удобно расставить пределы интегрирования выполним с помощью Matplotlib дополнительный чертеж.
Код для дополнительного чертежа
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import axes3d 
import numpy as np 
a, b, h = 2, 5, 25
# Создание поля для графика 
fig = plt.figure(constrained_layout=True, figsize=(7,5)) 
spec = fig.add_gridspec(ncols=2, nrows=1) 
spec_right = spec[1].subgridspec(ncols=1, nrows=5)
ax1 = fig.add_subplot(spec[0, 0], projection='3d') 
ax2 = fig.add_subplot(spec_right[1:4])
# ---------- График-1: Параболоид в 3D -----------
ax1.view_init(elev=20, azim=-40, roll=0)
# диапазоны по осям и сетка 
x_lims = [-1.5a, 1.5a] 
y_lims = [-1.5b, 1.5b] 
z_lims = [-1.0*h, 0]
X, Y = np.linspace(*x_lims, 50), np.linspace(*y_lims, 50) 
X, Y = np.meshgrid(X, Y)
def z_func(X, Y):
  return np.where((X2 / a2 + Y2 / b2) <= 1,
                  -h * (X2 / a2 + Y2 / b2) ,
                  -h)
Z = z_func(X,Y)
# Рисуем 3D поверхность 
ax1.plot_surface(X, Y, Z, edgecolor='royalblue', lw=0.5, 
                 rstride=3, cstride=3,
                 alpha=0.3) 
ax1.set_aspect('equal')
# Проекции на координатные плоскости 
ax1.contour(X, Y, Z, levels=[-h], zdir='z', offset=z_lims[0],
            colors=['black'])
# Подписи к графику 
ax1.text(x_lims[0]-1, y_lims[0], -1, r'', (0, 1, 0)) 
ax1.text(x_lims[0]-1, y_lims[0], -h, r'', (0, 1, 0))
ax1.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
        xlabel='X', ylabel='Y', zlabel='Z',
        xticks=np.linspace(*x_lims, 4))
# ------ График-2: Проекция на плоскость основания ------
# Эллипс 
phi = np.linspace(0, 2*np.pi, 50)
X = a * np.cos(phi) Y = b * np.sin(phi)
ax2.plot(X, Y, '--')
# Четвертинка 
X2 = np.linspace(0, a, 50) 
Y2 = b * np.sqrt(1 — X22 / a2)
ax2.plot(X2, Y2, color='orange') 
ax2.plot(X2, Y2*0, color='orange') 
ax2.fill_between(X2, Y2, 0, color='blue', alpha=.1)
# Аннотации 
ax2.text(0,-0.6, r'', rotation='vertical') 
ax2.text(a-0.2,-0.6, r'', rotation='vertical')
ax2.text(a + 0.2, 0, r'') ax2.text(a + 0.2, b, r'')
# Настройки отображения 
ax2.grid() 
ax2.set(xlabel=r'X', ylabel=r'Y')
plt.show()
Поскольку фигура симметричная, достаточно получить объем одной четверти, а потом умножить на 4. Это позволит выставить часть приделов равными 0, упростив тем самым выражение.
Таким образом для вычисления одной четверти объема:
- интеграл по dz берем в пределах от поверхности (плоскости) основания 
 до поверхности параболоида- , 
- в интеграле по dy смотрим на проекцию параболоида на плоскость основания (получается эллипс) и берем в пределах 
 от прямой
 до эллиптической границы основания- , 
- в интеграле по dx смотрим на проекцию основания на ось координат Ox берем в пределах от точки - до самой правой точки проекции 
 Объем всей фигуры тогда запишется как 
  
Вычислим его с помощью функции tplquad() библиотеки SciPy:
from scipy import integrate 
a, b, h = 2, 5, 25
# Подынтегральная фунция и пределы интегрирования 
f = lambda x,y,z: 1 
x_lims = [0, a] 
y_lims = [0, lambda x: b * (1 — x*2 / a**2)0.5] 
z_lims = [-h, lambda x, y: -h * (x2 / a2 + y2 / b*2)]
# Тройной интеграл 
V, err = integrate.tplquad(f, *x_lims, *ylims, *zlims)
print(f'Объем = {V * 4:.2f}')Объем зуба с параболическими обводами равен V=392.7 мм3 .
Продолжим совершенствование формы зуба и сделаем его более острым и изогнутым
Зуб варана v.3 (логарифм)
Уравнение поверхности и визуализация
Нам нужно «загнуть» центральную ось нашего параболического зуба. То есть, если посмотреть на рисунок ниже, — перейти от ситуации с левой верхней схемы к левой нижней.

Код рисунка
import sympy as sp
import numpy as np
from sympy.abc import x, y, z
from spb import *
from matplotlib.gridspec import GridSpec
from sympy.printing.latex import LatexPrinter
a, b, h = 2, 5, 25
la_print = LatexPrinter()  # рендеровщик LaTex формул
# Уравнение параболоида с логарифмическим сдвигом по "y"
left_part_eq = x**2/a**2 + (y + 2*sp.log(1-z))**2/b**2 + z/h
equations = [left_part_eq]
# Выразим "y" из уравнения
sol1, sol2 = sp.solve(equations, y, dict=True)
# Достанем программно логарифм из уравнения "загнутого" зуба
log_line = left_part_eq.args[2].args[1].args[0].args[1]
log_line
# Кривая смещения
eq_log = sp.Eq(y, -log_line)
# Уравнение параболической поверхности зуба
parab_left = x**2 / a**2 + y**2 / b**2 + z / h 
z_sol_parab = sp.solve([parab_left], z)
# Проекция параболоида
eq_parab = sp.Eq(z, z_sol_parab[z].subs({x:0}))
# Параметры отрисовки
params = {'aspect':(1,1), 'legend':False, 'show':False}
params_proj = {'aspect':(1,1), 'legend':False, 'show':False}
ranges_par = [(y, -4*b, 2*b), (z, -h, 0)]
ranges_log = [(y, -4*b, 2*b), (z, -h, 0)]
x_range = (x, -4*a, 4*a)
log_style = {"linestyles": "dashed", "linewidths": 1}
bacolor = 'dodgerblue'
# -------- График-1: Параболический зуб в профиль --------
# -- Аннотации
shift_arrows = []
for z_val in np.linspace(-0.2*h, -0.8*h, 3):
  y_val = -log_line.subs({z:z_val}).n()
  arrow = {'text': '', 'xy': (y_val, z_val), 
           'xytext': (0, z_val),
           'arrowprops': dict(arrowstyle="->")}
  shift_arrows.append(arrow)
# -- Парабола, вертикальное сечение зуба
p1 = plot_implicit(z < z_sol_parab[z].subs({x:0}),
                   *ranges_par, **params, 
                   xlabel='y', ylabel='z')
# -- Кривая куда смещаем, искривляя, ось
p1_logos = plot_implicit(eq_log, *ranges_par, color='red',
                         **params, annotations=shift_arrows, 
                         rendering_kw=log_style,)
# -- Кастомная легенда
legend_xy = (ranges_par[0][1], -2)
pretty_log_leg = la_print.doprint(sp.Eq(y,-log_line)).replace('log', 'ln')
text_legend = {'text': f'\n${pretty_log_leg}$',
               'xy': (0, 0), 'xytext': legend_xy, 
               'fontsize': 8,
               'bbox': dict(boxstyle="square", fc="w")}
# -- Текущая прямая ось зуба-парабалода
p1_os = plot_implicit(sp.Eq(y, 0), *ranges_par, **params, 
                      color="white",
                      rendering_kw={"linestyles": "-.", 
                                    "linewidths": 4})
# -- График-2: Проекции на основание параболического зуба --
# -- Эллипс, проекция на основание
p1_base = plot_implicit(parab_left.subs({z:0}) <= 1,
                        (y, ranges_par[0][1], ranges_par[0][2]), x_range,
                        markers=[{'args': [0, 0, 'x'], 
                                  'color':'white'}],
                        #annotations=[text_para],
                        **params_proj, color=bacolor)
# -- Стрелка направление смещения основания в
new_cent = -log_line.subs({z:-h}).n()  # центр нового основания
shift_arrow_base = {'text': '', 'xy': (new_cent, 0), 'xytext': (0, 0),
                    'arrowprops': dict(arrowstyle="->")}
size_arrow_base = {'text': '', 'xy': (new_cent, 3), 'xytext': (0, 3),
                   'fontsize': 7,
                    'arrowprops': dict(arrowstyle="|-|")}
pretty_log = la_print.doprint(log_line).replace('log', 'ln')
text_log = {'text': f'${pretty_log}$',
             'xy': (0, 0), 'xytext': (new_cent, 5)}
# -- Эллипс, смещенная в новое положение проекция
p1_base_sh = plot_implicit(sp.Eq(left_part_eq.subs({z:-h}),0),
                           (y, ranges_par[0][1], ranges_par[0][2]), x_range,
                           annotations=[shift_arrow_base, size_arrow_base, text_log],
                           rendering_kw={"linestyles": "--", "linewidths": 1},
                           markers=[{'args': [new_cent, 0, 'x'], 'color':'red'}],
                           color='gray', **params_proj,
                           xlabel='y', ylabel='x')
# ------- График-3: Логарифмический зуб в профиль --------
# -- Вертикальное сечение изогнутого зуба
p2 = plot_implicit(sp.And(y >= sol1[y].subs({x:0}), y <= sol2[y].subs({x:0})),
                       *ranges_log, **params)
# -- Центральная ось изогнутого зуба
p2_os = plot_implicit(eq_log, *ranges_par, **params, color="white",
                      rendering_kw={"linestyles": "-.", "linewidths": 4})
# -- Кривая смещения (куда таки сместили, изогнув, параболоид)
p2_logos = plot_implicit(eq_log, *ranges_log, color='red',
                         **params, rendering_kw=log_style)
# -- График-4: Проекции на основание логарифмического зуба --
# -- Эллипс, проекция на основание
p2_base = plot_implicit(left_part_eq.subs({z:-h}) <= 0,
                        (y, ranges_log[0][1], ranges_log[0][2]), x_range,
                        markers=[{'args': [0, 0, 'x'], 'color':'gray'}],
                        #annotations=[text_loga],
                        **params_proj, color=bacolor)
# -- Эллипс, старая проекция, откуда смещались
p2_base_sh = plot_implicit(sp.Eq(parab_left.subs({z:-h}),0),
                           (y, ranges_log[0][1], ranges_log[0][2]), x_range,
                           annotations=[shift_arrow_base],
                           rendering_kw={"linestyles": "--", "linewidths": 1},
                           markers=[{'args': [new_cent, 0, 'x'], 'color':'red'}],
                           color='gray', **params_proj,
                           xlabel='y', ylabel='x')
# Собираем изображения в сетку общего рисунка
gs = GridSpec(2, 2)
mapping = {
    gs[0, 0]: p1 + p1_logos + p1_os, # + p1_log_legend,
    gs[0, 1]: p1_base_sh + p1_base,
    gs[1, 0]: p2 + p2_os + p2_logos,
    gs[1, 1]: p2_base_sh + p2_base,
}
plotgrid(gs=mapping, size=(7,7), legend=False);
Это можно сделать, например, задав переменное смещение для y-координаты.
Способ-1 (канонический)
Уравнение параболоида имеет вид    . 
Такое мы использовали выше для задания формы зуба в разделе Зуб варана v.2 (параболоид).
Уравнение параболоида с смещением s по y имеет вид
 
И если s — это просто некоторое положительное число, то весь параболоид сдвинется вдоль оси Oy влево. А если величина s зависит от z, то есть принимает какое-то свое значение для каждого горизонтального сечения  параболоида, то мы имеем возможность сместить каждое сечение персонально на нужную нам величину. Схемы такого смещения показаны в правой части вышеприведенного рисунка. 
Опираясь на знание, как ведет себя логарифм, и чувство прекрасного, зададим s как функцию от z следующего вида:   
Тогда уравнение примет вид  
Подставим значения a, b и h и получим искомое уравнение боковой поверхности искривленного драконьего зуба:
   
Снизу, как обычно, объем зуба ограничен плоскостью .
Вводим полученное уравнение в сервисе desmos.com или аналогичном и наблюдаем изображение.

Способ 2 (параметрический)
Уравнение параболоида можно записать и в параметрическом виде, где аналогично способу-1 задать смещение s для y-координаты:
после подстановки найденной выше  
  
где параметры изменяются в следующих диапазонах:  , 
 
Изобразим драконий зуб программно:
import numpy as np
import sympy as sp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from sympy.abc import x, y, z
a, b, h = 2, 5, 25
# функция нелинейного сдвига
log_line = 2*sp.log(1-z)
# векторизуем функцию логарифмического смещения для расчета сетки ниже
v_log_line = np.vectorize(lambda tz: float(log_line.subs({z:tz}).n()))
# сетка для переменных-параметров
tz = np.linspace(0, -np.sqrt(h), 30)
phi = np.linspace(0, 2 * np.pi, 30)  # одна сторона зуба
tz, phi = np.meshgrid(tz, phi)
# подставим сетку в параметрическое уравнение параболоида со смещением
Z = -tz**2
X = (a / np.sqrt(h)) * tz * np.cos(phi)
Y = (b / np.sqrt(h)) * tz * np.sin(phi) - v_log_line(-tz**2)
# График
ax = plt.figure().add_subplot(projection='3d')
ax.view_init(elev=30, azim=10, roll=0)
ax.plot_surface(X, Y, Z, color='w', edgecolor='royalblue', 
                rstride=2, cstride=2, alpha=0.3)
# Диапазоны и настройки графика
x_lims = [-2*a, 2*a]
y_lims = [-3*b, 0]
z_lims = [-h, 0]
ax.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
       xlabel='X', ylabel='Y', zlabel='Z',
       xticks=np.round(np.linspace(*x_lims, 5),1),
       yticks=np.round(np.linspace(*y_lims, 6),1)
       )
ax.set_aspect('equal')
plt.show()
Объем искривленного зуба варана
Разрежем мысленно фигуру на элементарные, предельно тонкие, слои — с высотой dh и площадью основания S. Объем каждого такого слоя будет равен :
  
Из параметрических уравнений поверхности  . Найдем дифференциал  
 
Площадь основания слоя S будет изменяться в зависимости от высоты, на которой вырезаем слой. По форме это в любом случае будет эллипс, но с различными полуосями. Обозначим их за a’ и b’ (штрихи просто как штрихи, не производные) и вспомним формулу площади эллипса:
Опять же из параметрических уравнений поверхности явно видны выражения для полуосей (множители перед синусом и косинусом):
   ,  
   
Подставим их в формулу площади эллипса:
А ее, в свою очередь, подставим в выражение для элементарного объема
 
Берем интеграл от левой и правой частей
  
Интегрируем мы в данной статье все численно, не будем делать исключение и для этого интеграла. Подставим значения a,b и h и найдем объем с помощью функции quad() библиотеки SciPy
from scipy import integrate
import numpy as np
a, b, h = 2, 5, 25
# Подынтегральная функция
def f(t):
 return - 2 * np.pi * a * b * t**3 / h
# Интеграл
V, err = integrate.quad(f, -np.sqrt(h), 0)
print(f'V = {V}')V = 392.7 мм3.
Объем получился такой же, как для зуба в форме обычного параболоида. Что логично, поскольку мы лишь параллельно сдвинули слои в сторону, не изменив ни их форму, ни толщину. Аналогично тому, как если положить на стол колоду карт и сдвинуть верхние слои относительно нижних, объем колоды останется прежним (привет принципу Кавальери).
Объем искривленного зуба варана с полостью
В заключение учтем тот факт, что зуб варана имеет полость внутри [4]. Ориентируясь на продольный срез зуба, приведенный на Fig.3 из [4], опишем форму полости как аналогично параболическую с высотой 80% от высоты зуба и размерами основания в 30% от соответствующих размеров основания зуба. То есть с h = 20 мм, a = 0.6 мм, b = 1.5 мм. Подставив эти значения в формулу/код из предыдущего раздела получим Vполости = 28.3 мм3
И значит объем костной ткани зуба V = Vобщий - Vполости = 392.7 - 28.3 = 364.4 мм3.
Кальций. Сколько вешать в граммах
Зная объем и плотность вещества можно найти его массу. Объем в наличии, осталось определится с веществом и плотностью.
Зуб варана состоит из дентина и покрывающей его эмали. Ввиду того, что слой эмали у комодского дракона чрезвычайно тонкий — 20 микрометров [3] - в рамках данной модели им пренебрегаем и считаем, что весь найденный выше объем занят дентином.
Каково экспериментальное значение плотности дентина у комодских варанов, мне, к сожалению, найти не удалось (если кто-то из читателей владеет такой информацией, буду признателен). Возьмем в качестве приближения параметры зубных тканей крокодила — рептилии схожей с вараном по размеру и хищному способу добывания пищи.
Минералом, придающим прочность зубам крокодила, является гидроксиапатит [5].
Найдем, каков процент его массы составляет чистый кальций.
Химическая формула гидроксиапатита    [6], плотность  
 =3.14 г/см3 [7] .
Массы атомов, входящих в формулу:   ,  
, 
,  
 (где u — атомная единица массы).
Откуда общая масса рассматриваемой молекулы
  
= 40.078 ∙ 5 + (30.974 + 15.999 ∙ 4) ∙ 3 + (15.999 + 1.008) = 502.307 u
И значит на долю кальция по массе приходится    = (40.078 ∙ 5) / 502.307 = 0.3989 = 39.89 %
Содержание самого гидроксиапатита в дентине составляет   = 0.713 = 71.3 % . Остальное — органические компоненты [5].
Таким образом, чистого кальция в дентине содержится   = 0.713 ∙ 0.3989 = 0.2844 = 28.44 % .
Оставшиеся 42.9% — это суммарная доля фосфора, кислорода и водорода.
Массовые доли рассматриваемых компонентов вещества зуба показаны на диаграмме ниже. Минеральная часть выделена сплошной заливкой, органическая — заполнением точками.

Найдем массу дентина одного драконьего зуба.
   = 364.4 мм3 ∙ 3.14 ∙ 10-3 г/мм3 = 1.144 г
Значит масса чистого кальция в одном зубе комодского варана составит
  = 0.325 г
Поскольку всего зубов у дракона с острова Комодо 60 штук, для их полной замены необходимо 60 ∙  = 19.5 г кальция.
Если поделить это на 40 дней (период обновления зуба согласо [4]), то получим 0.488 гр./сутки = 488 мг/сутки.
Именно столько кальция требуется варану для обеспечения динамики регулярного обновления комплекса зубов.
Норма потребления по кальцию в целом будет выше, т.к. кроме зубов он содержится в костной ткани, необходим для сокращения скелетных мышц, передачи нервных импульсов и ряда других процессов в организме. К примеру, норма потребления кальция для человека в возрасте от 16 до 50 лет, когда роста зубов уже не происходит, составляет 800-1000 мг/сутки [8].
Режущая кромка
Режущая кромкаКромка v.1 (ровная)
Зуб комодского дракона замечателен не только своей формой, но и наличием укрепленной железом кромки [3]. На рисунке ниже она выделяется оранжевым цветом и волнистой формой.

Для начала опишем геометрию ровной кромки, без учета идущей по ней «волны».
Уравнение кривой и визуализация
Поскольку кромка зуба это линия, принадлежащая его поверхности, воспользуемся полученным ранее параметрическим представлением поверхности.  
Кромку будет образовывать множество точек со значениями углового параметра     и  
 (соответствующие противоположным краям), при том, что параметр t, отвечающий за положение по высоте, пробегает прежний диапазон значений  
.
Такое представление пригодится нам для трехмерной визуализации.
Однако в данном случае (и для будущего расчета длины кромки) можно обойтись и  двумерным графиком в центральной вертикальной плоскости (), содержащей оба противоположных режущих края. 
При фиксированном x = 0 в системе остается два уравнения:
Которые можно свести в одно, выразив t через z из первого уравнения и подставив во второе. Итого
  
Подставим известные значения b и h.
Для переднего края (обозначен на рисунке ниже L1) с    приходим к следующему уравнению в декартовых координатах:
   
Для заднего края (обозначен на рисунке ниже L2) с    — к уравнению: 
   
Построим графики:

Код для построения графиков
import numpy as np
import sympy as sp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from sympy.abc import x, y, z
a, b, h = 2, 5, 25
# Кривая смещения
log_line = 2*sp.log(1-z)
# Поле рисунка для построения графиков
fig = plt.figure(constrained_layout=True, figsize=(8,5))
spec = fig.add_gridspec(ncols=2, nrows=1)
spec_right = spec[1].subgridspec(ncols=1, nrows=5)
ax1 = fig.add_subplot(spec[0, 0], projection='3d')
ax2 = fig.add_subplot(spec_right[1:4])
# ------- График-1: зуб с выделенной режущей кромкой --------
# векторизуем функцию логарифмического смещения для расчета сетки ниже
v_log_line = np.vectorize(lambda tz: float(log_line.subs({z:tz}).n()))
# сетка для переменных-параметров
tz = np.linspace(0, np.sqrt(h), 30)
phi_0 = np.linspace(0, 2 * np.pi, 30)  # поверхность зуба
phi_1 = np.linspace((3/2)*np.pi, (3/2)*np.pi, 30)  # режущая кромка зуба
phi_2 = np.linspace(np.pi / 2, np.pi/ 2, 2)  # режущая кромка зуба
tz_0, phi_0 = np.meshgrid(tz, phi_0)
tz_1, phi_1 = np.meshgrid(tz, phi_1)
tz_2, phi_2 = np.meshgrid(tz, phi_2)
# График
ax1.view_init(elev=30, azim=-35, roll=0)
for tz, phi, col, lw in ([tz_0, phi_0, 'royalblue', 1], [tz_1, phi_1, 'r', 2], [tz_2, phi_2, 'r', 2]):
  # подставим сетку в параметрическое уравнение параболоида со смещением
  Z = -tz**2
  X = (a / np.sqrt(h)) * tz * np.cos(phi)
  Y = (b / np.sqrt(h)) * tz * np.sin(phi) - v_log_line(-tz**2)
  ax1.plot_surface(X, Y, Z, color='w', edgecolor=col, lw=lw,
                  rstride=5, cstride=5, alpha=0.3)
# Диапазоны и настройки графика
x_lims = [-2*a, 2*a]
y_lims = [-3*b, 0]
z_lims = [-h, 0]
ax1.set_aspect('equal')
ax1.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
        xlabel='X', ylabel='Y', zlabel='Z',
        xticks=np.round(np.linspace(*x_lims, 5),1),
        yticks=np.round(np.linspace(*y_lims, 6),1)
        )
# -------- График-2: режущая кромка на плоскости ---------
# Уравнение боковой поверхности изогнутого зуба
left_part_eq = x**2/a**2 + (y + log_line)**2/b**2 + z/h
equations = [left_part_eq]
sol1, sol2 = sp.solve(equations, y, dict=True)
print(sol1, sol2)
# Проекция вдоль Ox на плоскость Oyz
zy_section1 = sol1[y].subs({x:0})
zy_section2 = sol2[y].subs({x:0})
print(zy_section1, zy_section2)
# Диапазоны
zs = np.linspace(*z_lims, 30)
ys1 = [zy_section1.subs({z:val}).n() for val in zs]
ys2 = [zy_section2.subs({z:val}).n() for val in zs]
# График
ax2.plot(ys1, zs, label='$L_1$')
ax2.plot(ys2, zs, label='$L_2$')
# Настройки графика
ax2.set(xlim=y_lims, ylim=z_lims, xlabel='Y', ylabel='Z')
ax2.set_aspect('equal')
ax2.grid()
ax2.legend()
plt.show()Длина ровной режущей кромки
Длину кромки найдем при помощи криволинейного интеграла 1го рода.
Длина передней кромки:
 
Длина задней кромки вычисляется по той же схеме, отличаясь лишь знаками в выражении для :
  
И полная длина   
Возьмем данные интегралы численно и определим итоговую длину.
import sympy as sp 
from sympy.abc import x, y, z 
from scipy import integrate 
a, b, h = 2, 5, 25
# Кривая смещения оси 
log_line = 2*sp.log(1-z)
# Уравнение боковой поверхности изогнутого зуба 
left_part_eq = x**2/a**2 + (y + log_line)**2/b**2 + z/h 
equations = [left_part_eq] 
# разрешим относительно y
sol1, sol2 = sp.solve(equations, y, dict=True) 
print('Решения отностительно y: \n', sol1, sol2)
# Проекция вдоль Ox на плоскость Oyz 
zy_section1 = sol1[y].subs({x:0}) 
zy_section2 = sol2[y].subs({x:0}) 
print('Проекции: \n', zy_section1, zy_section2)
# Берем производные 
derivative_y1 = sp.diff(zy_section1, z) 
derivative_y2 = sp.diff(zy_section2, z) 
print('Их производные: \n', derivative_y1, derivative_y2)
# Подынтегральная функция для криволинейного интеграла 
def f(z_val, curve_name):
  derivative = {'L1': derivative_y1, 'L2': derivative_y2}
  return (1 + derivative[curve_name].subs({z:z_val}).n()**2) ** 0.5
# Интеграл 
length = integrate.quad(f, -h, 0, 'L1')[0] + integrate.quad(f, -h, 0, 'L2')[0] 
print(f'Длина кромки {length} мм')Таким образом длина ровного режущего края l=54.18 мм
Кромка v.2 (волнистая)
Циклоида путешественница
Для начала вспомним, как выглядит параметрическое уравнение циклоиды:
 
Если производящая окружность будет катиться не по горизонтальной прямой , а по кривой 
, то уравнение примет следующий вид: 
 
Как можно заметить, здесь использован тот же принцип с заданием функции смещения по y, что и при искривлении оси параболоида в разделе Зуб варана v.3 (логарифм).
Такая кривая является частным случаем квазитрохоиды (которая в общем случае описывает сочетание произвольного числа разновидностей поступательного и вращательного движений).
Демонстрационный пример
Рассмотрим демонстрационный пример — с опорной кривой вида , по которой отправим в путешествие производящую окружность радиусом R=0.5.
или
Что на графике будет выглядеть так:

Код для построения графика:
import matplotlib.pyplot as plt 
import numpy as np 
R = 0.5  # радиус производящей окружности
# Уравнение циклоиды, опирающейся на наклонную прямую 
x = lambda t: R * (t — np.sin(t))
line = lambda t: 0.1 * x(t)  # опорная прямая 
y = lambda t: R * (1 — np.cos(t)) + line(t)
# Подготовка данных для графика 
ts = np.linspace(0, 6*np.pi, 100) 
xs = [x(t) for t in ts] 
ys = [y(t) for t in ts] 
framework = [line(t) for t in ts]
# График 
fig, ax = plt.subplots()
ax.plot(xs, ys, label='L') 
ax.plot(xs, framework)
# Оформление графика
ax.grid() 
ax.legend()
plt.show()Уравнение волнистой режущей кромки и визуализация
Теперь направим производящую окружность квазитрохоиды вдоль режущей кромки зуба, уравнения для которой мы получили ранее в разделе Кромка v.1 (ровная).
Передняя кромка зуба описывается уравнением   
Значит параметрическое уравнение квазитрохоиды, опирающейся на эту кромку будет иметь вид
 
Подставим выражение для y1
 
Заменим z во втором уравнении на ее выражение через параметр t из первого уравнения
 
Получили параметрическое уравнение волнистой передней кромки зуба при радиусе производящей окружности R=0.5.
Такое значение радиуса обеспечит нам наглядность чертежа. Далее при расчете длины возьмем значение по экспериментальным данным (меньше на порядок).
Аналогично и для задней кромки.  Уравнение опорной кривой имеет вид  
И, значит, параметрическое уравнение для задней волнистой кромки:
 
Параметр t, пробегает значения от    до 0, где  
  — решение уравнения 
. 
Воспроизведем это в программном коде и графической форме.
Первым этапом находим уравнения для ровных передней и задней кромок.
import matplotlib.pyplot as plt 
import numpy as np 
import sympy as sp 
from sympy.abc import x, y, z, u 
a, b, h = 2, 5, 25
R = 0.5  # радиус производящей окружности
# Кривая смещения оси
log_line = 2*sp.log(1-z)
# Уравнение боковой поверхности изогнутого зуба
left_part_eq = x2/a2 + (y + log_line)2/b2 + z/h
equations = [left_part_eq]
# Разрешим относительно y
sol1, sol2 = sp.solve(equations, y, dict=True)
print('Решения относительно y: \n', sol1, sol2)
# Проекция вдоль Ox на плоскость Oyz
zy_section1 = sol1[y].subs({x:0})
zy_section2 = sol2[y].subs({x:0})
print('Проекции: \n', zy_section1, zy_section2)Получили выражения для передней и задней ровных кромок (в коде выше zy_section1 и zy_section2). Применим их в качестве направляющих для развертывания волнообразного режущего края.
# Параметрическое уравнение квазитрохоиды на направляющих кромках
Z = lambda t: R * (t - np.sin(t))
# Направляющие
line_1 = lambda t: float(zy_section1.subs({z:Z(t)}).n())
line_2 = lambda t: float(zy_section2.subs({z:Z(t)}).n())
Y_1 = lambda t: -R * (1 - np.cos(t)) + line_1(t)
Y_2 = lambda t: R * (1 - np.cos(t)) + line_2(t)
Определим пределы изменения переменных
# Границы и диапазоны
y_lims = [-3*b, 1]
z_lims = [-h, 1]
# при каком значении параметра z = -h (т.е. на нижнем крае)
equation = sp.Eq(R * (u - sp.sin(u)), -h)
bottom_u = float(sp.nsolve(equation, u, -1))
# z = 0 очевидно при u = 0
top_u = 0
print(f'Пределы изменения параметра u: от {bottom_u} до {top_u}')
# генерируем данные для построения графика
ts = np.linspace(bottom_u, top_u, 100)
ys_1, ys_2 = [Y_1(t) for t in ts], [Y_2(t) for t in ts]
zs_1, zs_2 = [Z(t) for t in ts], [Z(t) for t in ts]
framework_y1 = [float(zy_section1.subs({z:z_val}).n()) for z_val in zs_1]
framework_y2 = [float(zy_section2.subs({z:z_val}).n()) for z_val in zs_2]
Построим график
# График
fig, ax = plt.subplots()
ax.plot(ys_1, zs_1, label='$C_1$')
ax.plot(framework_y1, zs_1, '--', label='$L_1$')
ax.plot(ys_2, zs_2, label='$C_2$')
ax.plot(framework_y2, zs_2, '--', label='$L_2$')
# параметры отображения графика
ax.set(xlim=y_lims, ylim=z_lims, xlabel='Y', ylabel='Z')
ax.set_aspect('equal')
ax.grid()
ax.legend()
plt.show()
Длина волнистой режущей кромки
Длина кривой L, заданной параметрически в зависимостями  и 
, где  
 , будет вычисляться по формуле
  
Займемся ее вычислением.
Зададимся уравнением боковой поверхности зуба и найдем уравнения для передней и задней кромок. Эту часть кода мы уже использовали в предыдущих разделах. Единственным отличием здесь будет значение радиуса порождающей окружности. Установим его равным R=0.01 мм в соответствии с экспериментальными данными [3]
Уравнения волнистых режущих кромок (Sympy)
import sympy as sp 
from sympy.abc import x, y, z, u 
from scipy import integrate 
a, b, h, R = 2, 5, 25, 0.01
# Кривая смещения оси 
log_line = 2*sp.log(1-z)
# Уравнение боковой поверхности изогнутого зуба 
left_part_eq = x**2/a**2 + (y + log_line)**2/b**2 + z/h 
equations = [left_part_eq] 
# Разрешим относительно y
sol1, sol2 = sp.solve(equations, y, dict=True) 
print('Решения отностительно y: \n', sol1, sol2)
# Проекция вдоль Ox на плоскость Oyz 
zy_section1 = sol1[y].subs({x:0}) 
zy_section2 = sol2[y].subs({x:0}) 
print('Проекции: \n', zy_section1, zy_section2Используем полученные как проекции выражения для передней и задней ровных кромок в качестве направляющих для развертывания волнообразного режущего края
# Параметрические уравнения передней и задней кромки 
zu = R * (u — sp.sin(u)) 
# Опорные кривые
line_1 = zy_section1.subs({z:zu}) 
line_2 = zy_section2.subs({z:zu}) 
Y_1 = -R * (1 — sp.cos(u)) + line_1 
Y_2 = R * (1 — sp.cos(u)) + line_2После чего посчитаем необходимые производные и возьмем интеграл
# Их производные по параметру 
derivative_zu = sp.diff(zu, u) 
derivative_yu1 = sp.diff(Y_1, u) 
derivative_yu2 = sp.diff(Y_2, u) 
print('Их производные: \n', derivative_zu, derivative_yu1, derivative_yu2) 
# Подынтегральная функция для криволинейного интеграла 
def f(u_val, curve_name):
  derivative = {'L1': derivative_yu1, 'L2': derivative_yu2}
  return (derivative_zu.subs({u:u_val}).n()**2 +
          derivative[curve_name].subs({u:u_val}).n()**2) ** 0.5
  
# Пределы 
# при каком значении параметра z = -h (т.е. на нижнем крае) 
equation = sp.Eq(R * (u — sp.sin(u)), -h) 
bottom_u = float(sp.nsolve(equation, u, -1))
# z = 0 очевидно при u = 0 
top_u = 0
# Интеграл
length = integrate.quad(f, bottom_u, top_u, 'L1')[0] + integrate.quad(f, bottom_u, top_u, 'L2')[0]
print(f'Длина кромки {length} мм')Итого, длина волнообразной режущей кромки составляет 66.47 мм. Циклоидные дуги увеличили длину кромки на 23%.
Объем волнистой режущей кромки
Из [3] можно заключить, что характерный размер идущего по кромке оранжевого канта в поперечном сечении составляет 20 мкм. Представим его как цилиндр диаметром D = 20 мкм, уложенный по найденной выше линии волнистой кромки и рассчитаем объем.
Радиус сечения r = D / 2 = 10 мкм = 0.01 мм.
Площадь сечения   
Объем цилиндра   
Итак, объем железосодержащего "шнура", образующего режущую кромку, равен V = 0.021 мм3.
Железо оранжевой кромки
Согласно исследованию [3] оранжевая режущая кромка зуба варана состоит из ферригедрита. Формула ферригидрита:     [9] Плотность  
 =3.96 г/см3 = 3.96∙10-3 г/мм3  [10] 
Массы атомов, входящих в формулу:  ,  
,  
 (где u — атомная единица массы).
Откуда общая масса рассматриваемой молекулы    = 5 (55.845 ∙ 2 + 15.999 ∙ 3) + 9 * (1.008 ∙ 2 + 15.999) = 960.57 u
Массовая доля Fe составляет
    = (10 ∙ 55.845) / 960.57 = 0.581 = 58.1 %
Масса чистого железа в покрытии режущей кромки зуба комодского варана составит
  = 0.581 ∙ 3.96∙10-3 ∙ 0.021 = 0.048 ∙ 10-3 г = 0.048 мг
Масса железа на кромках полного комплекта из 60 зубов    = 60 ∙ 0.048 мг = 2.88 мг
Зная, что период замены зуба у варана составляет 40 суток [4], найдем, сколько железа нужно в день при условии его равномерного расходования на покрытие формирующихся режущих кромок.  / 40   = 2.88 мг/ 40 суток = 0.072 мг/сутки. 
Не так уж много. Если сравнить с человеком, то средняя суточная доза потребления железа для взрослых мужчин и пожилых людей обоих полов составляет 8 мг (для женщин в 2-3 раза больше) [11]. При этом оно расходуется в основном на не связанные напрямую с зубами нужды организма.
Варану на покрытие острых кромок своих зубов нужно порядка 1% от такого объема.
Заключение
В рамках рассмотренных моделей мы получили, что варану для обеспечения динамики обновления зубов необходимо порядка 488 мг кальция и 72 мкг железа в сутки. Для сравнения: это примерно половина суточной потребности человека в кальции и менее одного процента суточной потребности человека в железе соответственно. Цифры не астрономические. И для комодского дракона, с его охотничьими способностями и положением на вершине пищевой цепи своего региона, вполне достижимые.
Что касается зубов человека, биоинженерные технологии понемногу идут в сторону создания (включения) механизма выращивания новых зубов взамен поврежденных или отсутствующих ([12] или подробнее в [13]). Но, думаю, еще долгое время варану конкуренция в этом вопросе с нашей стороны не грозит.
Источники
- https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0295002 
- https://www.handbookofmineralogy.org/pdfs/hydroxylapatite.pdf 
- https://www.consultant.ru/document/cons_doc_LAW_383904/def867cf1af3e9abc3ff8fe0463d65565a4b6264/ 
- https://hepd.pnpi.spb.ru/hepd/events/abstract/2024/Getalov_16_04_24.pdf 
- https://www.sciencedirect.com/science/article/pii/S2352320423000044 
 
           
 
BrazilMaterializator
Поскольку всего зубов у дракона с острова Комодо 60 штук, для их полной замены необходимо 60 ∙ = 14.1 г кальция.
  = 14.1 г кальция.
У меня нет под рукой компиляторов. Но 0,3*60 это уже 18.
shbma Автор
Действительно. Благодарю за внимательный взгляд, исправил.