This page gives more technical details on how the plate constant calculator is set up and works.
All the necessary galaxy information is in a simple text file, one line per galaxy. There are six pieces, in three pairs, at a minimum.
The first two fields are the catalog (NGC, IC, UGC, etc.) and the galaxy's number in it. Splitting the two up makes it easier to search by just the number.
The next two fields are the galaxy's astronomical coordinates, it's right ascension and declination, precessed to the time of observation. It would be possible to use a common reference epoch (like 2000.0) but then the program would need to know the time of observation so it could perform the precession. One could work in the reference epoch, but the rotation will be slightly different since precession causes a small amount of field rotation.
The final two fields are the galaxy's X and Y pixel coordinates measured in a frame.
I also use a flag field of a single byte, currently only to have the program include the galaxy in the fit or not. After the fit is based on two dozen or thirty galaxies it doesn't change very much with the addition of another galaxy, as the list is being built, even if the galaxy makes the fit worse for some reason. It's easier to flag it so it doesn't get included in the fit, while checking the numbers to see why it's discrepant. In my three fields, there was only one galaxy not near a frame edge or corner which is still off for unknown reasons. With the galaxies ordered by right ascension this means the flagged galaxies are all at the top or bottom of the list.
After reading in all the data, the calculator then does two linear least squares fits. One fits the pixel X values to the right ascension coordinates, while the other does the same for the Y values and declinations. The former needs to include a Cos(dec) term for its scale (slope of the fit) to be the same as that in the orthogonal direction. (Other than that, the slight convergence of the lines of constant right ascension towards the poles is not taken into account.) The two scales are then averaged.
The fit slope is in the units of pixels per radian (or the reciprocal, depending on how you do it), which can be converted into radians per mm using the manufacturer's spec (usually in the EXIF data) for the sensor's number of pixels per mm. The equation Scale = 1/F then yields the focal length (F) from the scale in angular rather than physical units.
The x- and y-intercepts that the two least squares fits produce gives the coordinates of the top and left sides of the frame (X and Y equal zero); knowing the scale and the number of pixels across and vertically, the coordinates of the bottom and right sides of the frame. It's convenient to denote these as the UL (upper left) and LR (lower right) corners of the frame.
The program can now internally generate predicted X and Y pixel coordinates from the astronomical coordinates using the two fits. By differencing these from the actual X and Y pixel coordinates, the RMS fit error can be calculated.
At this point the fit is likely to be mediocre because the camera is not perfectly level, causing the frame to be rotated slightly from the sky. The amount of this rotation is easily determined by converting both the predicted and actual pixel coordinates into polar coordinates. For an origin I thought it best to use the fit center rather than the frame center, because this makes the calculation independent of centering variations from frame to frame. The fit center is simply the mean X and Y values, which are calculated as a matter of course anyway during the linear regression part of the calculation; in other words, they are already variables that are present and don't have to be determined.
In polar coordinates rotation manifests itself as a systematic difference in the angles between the predicted and actual galaxy locations. The radial differences will be small. It is then a simple matter of averaging the angle differences to produce the amount of frame rotation. The actual X and Y pixel measures can then be rotated by this amount around the fit center to square the frame with the sky. After this is done, the least squares fit is re-done with the rotated pixel positions, and the slope, UL and LR corners, and RMS fit error are all re-calculated as before to give the final fit and the final values.
The only hitch in all this is that the average of the angle differences can be thrown off by a galaxy or two very close to the fit center. There is a sort of angle magnification effect as one gets closer and closer to the center, such that a couple of pixel difference (which is typical) between actual and predicted positions at a distance of only a few dozen pixels from the center is about a tenth of a radian (several degrees), whereas the same several pixels at many hundreds or more than a thousand pixels distance from the center is the fraction of a degree we want to average in. So, besides having a "dead zone" right near the center, I weighted the angle differences as each was added into the sum which is formed to calculate the average; the weighting is very simple (linear) such that the weight increases the farther a galaxy is from the fit center. Once the weight gets to 1 (the max) at some suitable distance, I then ramp it back down towards zero at larger distances, which include the slightly distorted images along the edges and in the frame corners where the errors are larger to begin with. With this little nuance addressed the rotation angle works as desired. In particular, different fields taken at the same time give the same result for the camera's tilt, even though they of course have entirely different galaxies and fit centers.