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]
.