First I shall begin by laying out some notation (I shall be using the conventions that are used by both DLMF and Mathematica which occasionally differ from the standard literature):
Let $\Lambda:=\mathbb{Z}\omega_1+\mathbb{Z}\omega_3$ (note that, as per DLMF and Mathematica, we use $\omega_3$ in place of the standard $\omega_2$, but it should be considered to be defined the same) be the period lattice generated by $\omega_1,\omega_3\in\mathbb{C}$. Furthemore, let $$\theta_1\left(z,q\right):= 2\sum_{n=0}^{\infty}(-1)^nq^{\left(n+\frac{1}{2}\right)^2}\sin\left((2n+1)z\right),\quad\eta(\tau):=q^{1/24}\prod_{n=1}^{\infty}\left(1-q^n\right),$$ be the Jacobi theta function of order $1$, with elliptic nome $q=e^{\pi i \tau}$, and Dedekind eta function respectively, with $\tau:=\frac{\omega_3}{\omega_1}$ and $q:=e^{2\pi i\tau}$ the elliptic nome for $\eta$. Let $\wp$ be the standard Weierstrass elliptic function satisfying the differential equation $$\wp'(z)^2=4\wp^3(z)-g_2\wp(z)-g_3.$$ It should be worth noting (for anyone intending to computationally verify the proceeding formulae) that Mathematica uses $\wp$ with parameters $g_2,g_3$ as opposed to $\omega_1,\omega_3$.
Given some elliptic curve $E$ of conductor $N$, with positive rank, the canonical height of some generator point $P=(x,y)\in E\left(\mathbb{Q}\right)$ is given by \begin{equation}\tag{1}\hat{h}\left(P\right)=2\log\left|\frac{\pi \cdot\eta(\tau)^3}{\omega_1 \theta_1\left(\frac{\pi z_0}{2\omega_1},q\right)}\right|+\frac{2\pi\cdot\Im\left(\frac{z_0}{2\omega_1}\right)^2}{\Im(\tau)}-\sum_{\substack{p\in\mathbb{P},\\ p\mid N}}\frac{k_p}{c_p}\log(p),\end{equation} where $c_p$ is the Tamagawa number associated to each bad prime $p$, $z_0:=\wp^{-1}\left(\tilde{x}\right)$, and $k_p\in\mathbb{Z}_{\geq 0}$ (conjecturally). We note that an elliptic curve $y^2+a_1xy+a_3y=x^3+a_2x^2+a_4x+a_6$ should be transformed to short Weierstrass form in the following particular way:
After multiplying by $4$, one should add $(a_1x+a_3)^2$ on either side, transforming the LHS into a perfect square, which should now be relabelled as $y^2$. One should then proceed to transform $x\mapsto x-\frac{\left(a_1^2+4a_2\right)}{12}$ as to make the $x^2$ term vanish. We then take $g_2,g_3$ to be the coefficients of the $x$ and constant term (accounting for the negatives of course). When inputting our $x$-coordinate into $\wp^{-1}$ for $z_0$, we should input the transformed value, namely given a starting generator point $P=(x,y)$, we use $z_0=\wp^{-1}\left(x+\frac{\left(a_1^2+4a_2\right)}{12}\right)=:\wp^{-1}(\tilde{x})$.
To demonstrate, take the elliptic curve given by LMFDB label 102.a1: $$y^2+xy=x^3+x^2-2x$$ which has generator point $P=(2,2)$ with canonical height $\hat{h}\left(P\right)\approx 0.143253\dots$. We transform it to the elliptic curve $y^2=4x^3-\frac{121}{12}x+\frac{845}{216}$. Taking $g_2=\frac{121}{12},g_3=-\frac{845}{216},z_0=\wp^{-1}\left(2+\frac{5}{12}\right)$ we have:
i.e. $k_2=k_3=1$ and $k_{17}=0$.
The goal of this question is to determine why the values $k_p$ arise, and ideally, determine a general formula for them.
To arrive at equation $(1)$ we essentially set up four recurrence relations:
Assume our elliptic curve is in short Weierstass form. Take $P=(1,0)$, then define the sequence of rational points $2^i P:=(X(i),Y(i))$. Now, $X(i)$ satisfies some recurrence relation in terms of $X(i-1)$ and $Y(i-1)$ that depends on the given elliptic curve. We then write $X(i)=\frac{x_1(i)}{x_2(i)}$ and $Y(i)=\frac{y_1(i)}{y_2(i)}$ with $x_1,x_2,y_1,y_2\in\mathbb{Z}$ as reduced as possible (this is the first place that could explain where $k_p$ arise as in general our recurrence relations need not produce reduced values, but there is likely something from the theory of elliptic divisibility sequences that can help here). Plugging this into the recurrence relation for $X(i)$ and comparing numerators and denominators, we can obtain some rather horrible recurrence relations for $x_1$ and $x_2$, and similarly one for $y_1$ and another $y_2$ by using the fact that the points are all on the same line. Now assuming $x_1(i)>x_2(i)$, the height is then $$\lim_{n\to\infty}\frac{\log(x_2(n))}{4^{n+1}}.$$ Now using the fact that with $X(i)=\wp\left(2^i z_0\right)$, we have $Y(i)^2=\wp'\left(2^i z_0\right)^2$ for some $z_0$ we can solve a recurrence relation to arrive at $$x_2(n)=x_2(1)^{c^{n-a}}\prod_{i=1}^{n-1}\left(Y(i)^2\right)^{4^{n-1-i}}$$ where $a$ is some integer and $c$ is the constant given by whatever the limiting value of $\frac{x_2(i)}{x_2(i-1)^4 Y(i-1)^2}$ is. This constant is intrinsically related to the $k_p$ values. Through all this, we essentially arrive at the fact that the height is given by the sum $$\tag{2}\frac{1}{2}\sum_{i=1}^{\infty}\frac{\log\left|\wp'\left(2^i z_0\right)\right|}{4^i}$$ plus some extra $\log$ terms that are (conjecturally) rational combinations of logarithms of the initial $x$ and $y$ coordinates of the generating points. The equation $(2)$ was determined by Tate in a letter to Serre (see here, here and here) albeit in a slightly different, seemingly less precise, form. Now using the fact that $$\wp'(z,\tau)=-\frac{\pi^3}{\omega_1^3}\eta^9(\tau)\cdot\frac{\theta_1\left(\frac{\pi z}{\omega_1},q\right)}{\theta_1\left(\frac{\pi z}{2\omega_1},q\right)^4},$$ it is then possible to evaluate $(2)$ exactly to arrive at equation $(1)$, which I believe is novel as I have been unable to locate it in the literature.
For the sake of computational verification, I have computed $k_p$ for some rank $1$ and $2$ elliptic curves which can be found in this spreadsheet. I conjecture that $k_p=0$ if $c_p=\mathrm{ord}(N)=\mathrm{ord}(\Delta)=\mathrm{ord}(j)=1$. For most cases, it appears that $k_p = \mathrm{ord}(\Delta)-\mathrm{ord}(N)$, but cases like $117.a4$ break this significantly.
