I'll suppose that the transformation from longitude/latitude $w=\begin{pmatrix}u\\v\end{pmatrix}$ to page coordinates $z=\begin{pmatrix}x\\y\end{pmatrix}$ is of the form $z=\rho R_\theta w + z_0$ where $R_\theta=\begin{pmatrix}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{pmatrix}$ and $z_0=\begin{pmatrix}x_0\\y_0\end{pmatrix}$, for some constants $\rho,\theta,x_0,y_0$. This seems to be what you want, since this transformation is "rotation anticlockwise by $\theta$, then scaling by $\rho$, then translation by $z_0$.
Now suppose for $j=1,2$, we are given $w_j,z_j$ (the two lat/long pairs and their corresponding points) so that $z_j=\rho R_\theta w_j+z_0$.
We wish to compute $\rho, \theta,z_0$ using this data. Since $\rho$ is the scaling factor, we have $\rho=\mathrm{dist}(z_1,z_2)/\mathrm{dist}(w_1,w_2)$, where $\mathrm{dist}$ is the usual Euclidean distance. Similarly, we can compute the angle $\theta$ by calculating the angle that $w_2-w_1$ makes with $z_2-z_1$, and then $z_0=z_1-\rho R_\theta w_1$.
This gives you all of the constants in the general formula, which you can use to find the coordinates corresponding to other lat/long pairs.
Here is some Haskell code which hopefully implements the above scheme and calculates the answer for your example.
neg z = [-x,-y] where [x,y]=z vplus a b = [a1+b1,a2+b2] where [a1,a2]=a [b1,b2]=b vdiff a b = a `vplus` neg b dotprod a b = a1*b1+a2*b2 where [a1,a2]=a [b1,b2]=b crossprod a b = a1*b2-a2*b1 where [a1,a2]=a [b1,b2]=b scalmult s w = [s*u,s*v] where [u,v]=w mag z = sqrt $ x*x+y*y where [x,y]=z dist a b = mag $ a `vdiff` b unitvec z = (1/(mag z)) `scalmult` z cos_sin a b = [dotprod a' b', crossprod a' b'] where [a',b']=[unitvec a,unitvec b] anglefromcs cs = sign * acos c where [c,s]=cs sign | s >= 0 = 1 | otherwise = -1 angle v w = anglefromcs cos_sin v w rho w1 w2 z1 z2 = (dist z1 z2) / (dist w1 w2) theta w1 w2 z1 z2 = angle (w2 `vdiff` w1) (z2 `vdiff` z1) rot theta w = [cos theta * u - sin theta * v, sin theta * u + cos theta * v] where [u,v]=w z0 w1 w2 z1 z2 = z1 `vdiff` ( rho' `scalmult` (rot theta' w1 ) ) where theta' = theta w1 w2 z1 z2 rho' = rho w1 w2 z1 z2 z w1 w2 z1 z2 w = rho' `scalmult` (rot theta' w) `vplus` z0' where rho' = rho w1 w2 z1 z2 theta' = theta w1 w2 z1 z2 z0' = z0 w1 w2 z1 z2 main = print z [37.725918,-122.223587] [37.619473,-122.373886] [307,267] [155,400] [37.7,-122.3]
The output is [226.55468797299545,303.8562915341063]
.