计算图像旋转角度

import astropy.io.fits as fits
import math

headname1= "TS-1015-500_18_1-offset_1024-w_4096-h_4096-Mono16.fits"
img1 = fits.open("/Users/zhouxu/zhaoyuanyuan/tmp/"+headname1)
h1 = img1[0].header
ra10  =  h1["CRVAL1"]
dec10  =  h1["CRVAL2"]
x10  =  h1["CRPIX1"]
y10  =  h1["CRPIX2"]
a1   =   h1["CD1_1"]
b1   =   h1["CD1_2"]
c1   =   h1["CD2_1"]
d1   =   h1["CD2_2"]

headname2 = "TS-990-500-D1_21_1-offset_1024-w_4096-h_4096-Mono16.fits"
img2 = fits.open("/Users/zhouxu/zhaoyuanyuan/tmp/"+headname2)
h2 = img2[0].header
ra20  =  h2["CRVAL1"]
dec20  =  h2["CRVAL2"]
x20  =  h2["CRPIX1"]
y20  =  h2["CRPIX2"]
a2   =   h2["CD1_1"]
b2   =   h2["CD1_2"]
c2   =   h2["CD2_1"]
d2   =   h2["CD2_2"]

'''
ra =ra10+a1*(x1-x10)+b1*(y1-y10)
dec=dec10+c1*(x1-x10)+d1*(y1-y10)

and 

ra =ra20+a2*(x2-x20)+b2*(y2-y20)
dec=dec20+c2*(x2-x20)+d2*(y2-y20)

===>

ra10+a1*(x1-x10)+b1*(y1-y10)=ra20+a2*(x2-x20)+b2*(y2-y20)
dec10+c1*(x1-x10)+d1*(y1-y10)=dec20+c2*(x2-x20)+d2*(y2-y20)

x1 = k1 + k2*(x2-x20) + k3*(y2-y20)

k1 = x10+( (ra20-ra10)/b1-(dec20-dec10)/d1 )/(a1/b1-c1/d1) 
k2 = (a2/b1-c2/d1)/(a1/b1-c1/d1)
k3 = (b2/b1-d2/d1)/(a1/b1-c1/d1)

y1 = k4 + k5*(x2-x20) + k6*(y2-y20)
k4 = y10+( (ra20-ra10)/a1-(dec20-dec10)/c1 )/(b1/a1-d1/c1) 
k5 = (a2/a1-c2/c1)/(b1/a1-d1/c1)
k6 = (b2/a1-d2/c1)/(b1/a1-d1/c1)

if:

x1 = a + cosA*x2 + sinA*y2
y1 = d - sinA*x2 + cosA*y2

then

sinA = k3  or -1*k5
cosA = k2  or k6

'''

k1 = x10+( (ra20-ra10)/b1-(dec20-dec10)/d1 )/(a1/b1-c1/d1)
k2 = (a2/b1-c2/d1)/(a1/b1-c1/d1)
k3 = (b2/b1-d2/d1)/(a1/b1-c1/d1)

k4 = y10+( (ra20-ra10)/a1-(dec20-dec10)/c1 )/(b1/a1-d1/c1)
k5 = (a2/a1-c2/c1)/(b1/a1-d1/c1)
k6 = (b2/a1-d2/c1)/(b1/a1-d1/c1)

deg=180.0/math.pi

A1 = math.asin(k3)*deg
A2 = math.asin(-1*k5)*deg
A3 = math.acos(k2)*deg
A4 = math.acos(k6)*deg


print("Rotation angle:",A1,A2,A3,A4)