# PI value
PI <- 3.141592

# four corners of CCD camera, on the satellite XY plane 
# Suzaku CCD camera has 18'x18' (0.3deg x 0.3 deg) square 
# field of view.  Namely, X,Y positions of the 
# four corners are as follows (in dgree).
X1 <-  0.15
Y1 <-  0.15
X2 <-  0.15
Y2 <- -0.15
X3 <- -0.15
Y3 <- -0.15
X4 <- -0.15
Y4 <-  0.15

print("X1")
print(X1)
print("Y1")
print(Y1)

print("X2")
print(X2)
print("Y2")
print(Y2)

print("X3")
print(X3)
print("Y3")
print(Y3)

print("X4")
print(X4)
print("Y4")
print(Y4)


vect1_x <- X1/180*PI
vect1_y <- Y1/180*PI
vect1_z <- sqrt(1.0-vect1_x**2-vect1_y**2)

vect2_x <- X2/180*PI
vect2_y <- Y2/180*PI
vect2_z <- sqrt(1.0-vect2_x**2-vect2_y**2)

vect3_x <- X3/180*PI
vect3_y <- Y3/180*PI
vect3_z <- sqrt(1.0-vect3_x**2-vect3_y**2)

vect4_x <- X4/180*PI
vect4_y <- Y4/180*PI
vect4_z <- sqrt(1.0-vect4_x**2-vect4_y**2)

Vin1 <- c(vect1_x,vect1_y,vect1_z)
Vin2 <- c(vect2_x,vect2_y,vect2_z)
Vin3 <- c(vect3_x,vect3_y,vect3_z)
Vin4 <- c(vect3_x,vect4_y,vect4_z)

# We need three rotation matrix
a <- matrix(0,nrow=3, ncol=3)
b <- matrix(0,nrow=3, ncol=3)
c <- matrix(0,nrow=3, ncol=3)

# rotation around Z
phi = 281.00
a[1,1] <- cos(phi/180*PI)
a[1,2] <- -sin(phi/180*PI)
a[2,1] <- sin(phi/180*PI)
a[2,2] <- cos(phi/180*PI)
a[3,3] <- 1

# rotation around Y
theta <- 94.08
b[1,1] <-  cos(theta/180*PI)
b[1,3] <-  sin(theta/180*PI)
b[3,1] <-  -sin(theta/180*PI)
b[3,3] <-  cos(theta/180*PI)
b[2,2] <- 1

# rotation around Z
psi = 184.47
c[1,1] <- cos(psi/180*PI)
c[1,2] <- -sin(psi/180*PI)
c[2,1] <- sin(psi/180*PI)
c[2,2] <- cos(psi/180*PI)
c[3,3] <- 1


# Calculate the output direction vector
Vout1 <- a %*% b %*% c %*% Vin1
Vout2 <- a %*% b %*% c %*% Vin2
Vout3 <- a %*% b %*% c %*% Vin3
Vout4 <- a %*% b %*% c %*% Vin4

xout_1 <- Vout1[1]
yout_1 <- Vout1[2]
zout_1 <- Vout1[3]

xout_2 <- Vout2[1]
yout_2 <- Vout2[2]
zout_2 <- Vout2[3]

xout_3 <- Vout3[1]
yout_3 <- Vout3[2]
zout_3 <- Vout3[3]

xout_4 <- Vout4[1]
yout_4 <- Vout4[2]
zout_4 <- Vout4[3]

# Convert from the output direction vector
# to equatorial coordinate
alpha1 <- atan(yout_1/xout_1)/PI*180.0
if (alpha1 <0 ) {alpha1 <- alpha1+360}
delta1 <- atan(zout_1/sqrt(xout_1**2+yout_1**2))/PI*180.0

alpha2 <- atan(yout_2/xout_2)/PI*180.0
if (alpha2 <0 ) {alpha2 <- alpha2+360}
delta2 <- atan(zout_2/sqrt(xout_2**2+yout_2**2))/PI*180.0

alpha3 <- atan(yout_3/xout_3)/PI*180.0
if (alpha3 <0 ) {alpha3 <- alpha3+360}
delta3 <- atan(zout_3/sqrt(xout_3**2+yout_3**2))/PI*180.0

alpha4 <- atan(yout_4/xout_4)/PI*180.0
if (alpha4 <0 ) {alpha4 <- alpha4+360}
delta4 <- atan(zout_4/sqrt(xout_4**2+yout_4**2))/PI*180.0

print("alpha1")
print(alpha1)
print("delta1")
print(delta1)

print("alpha2")
print(alpha2)
print("delta2")
print(delta2)

print("alpha3")
print(alpha3)
print("delta3")
print(delta3)

print("alpha4")
print(alpha4)
print("delta4")
print(delta4)