# 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
# This corresponds to 23.439 deg rotation around X-axis
a <- matrix(0,nrow=3, ncol=3)
theta <- 23.439
a[1,1] <- 1
a[2,2] <-  cos(theta/180*PI)
a[2,3] <-  sin(theta/180*PI)
a[3,3] <-  cos(theta/180*PI)
a[3,2] <-  -sin(theta/180*PI)

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

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

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

print("lambda")
print(lambda)
print("beta")
print(beta)