# PI value
PI <- 3.141592

# Input eqatorial coordinate
alpha <- 299.590316
delta <- 35.201606

print("alpha")
print(alpha)
print("delta")
print(delta)

# Input direction vector
xin <- cos(delta/180*PI)*cos(alpha/180*PI)
yin <- cos(delta/180*PI)*sin(alpha/180*PI)
zin <- sin (delta/180*PI)

Vin <- c(xin,yin,zin)

# Define the conversion matrix
# We need three rotaions around Z, Y and X
a <- matrix(0,nrow=3, ncol=3)
b <- matrix(0,nrow=3, ncol=3)
c <- matrix(0,nrow=3, ncol=3)

# rotation 266.405 around Z
phi = 266.405
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 28.93617 around Y
theta <- 28.93617
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 58.59866 around X
psi <- 58.59866
c[1,1] <- 1
c[2,2] <-  cos(psi/180*PI)
c[2,3] <-  sin(psi/180*PI)
c[3,3] <-  cos(psi/180*PI)
c[3,2] <-  -sin(psi/180*PI)


# Calculate the output direction vector
Vout <- c %*% b %*% a %*% Vin

xout <- Vout[1]
yout <- Vout[2]
zout <- Vout[3]

# Convert from the output direction vector
# to Galactic  coordinate
l <- atan(yout/xout)/PI*180.0
if(l <0){l <- l+360}
galb   <- atan(zout/sqrt(xout**2+yout**2))/PI*180.0

print("l")
print(l)
print("b")
print(galb)