Aitoff projection on R part 2

Aitoff projection に変換するコードと作図用の外枠を描くコードをメモ（補助線は省略）。

=============

d <- read.csv("nova_catal_2010.csv", header=T)

 pi <- 3.14159265
 rad <- pi/180
 deg <- 180/pi

#-----外枠
   L1 <- seq(-180, 180, length=100)
   B1 <- rep(-90:-90, length=100)
  
   B2 <- seq(-90, 90, length=100)
   L2 <- rep(180:180, length=100)
  
   L3 <- seq(180, -180, length=100)
   B3 <- rep(90:90, length=100)
  
   B4 <- seq(90, -90, length=100)
   L4 <- rep(-180:-180, length=100)
         
   outer_l <- c(L1,L2,L3,L4)
     outer_l <- outer_l*rad
   outer_b <- c(B1,B2,B3,B4)
     outer_b <- outer_b*rad

outer_a <- acos(cos(outer_b)*cos(outer_l/2))
sinc_outer_a <- sin(outer_a) / outer_a

outer_x <- (2*cos(outer_b)*sin(outer_l/2))/sinc_outer_a
outer_y <- sin(outer_b)/sinc_outer_a

outer_x <- outer_x*deg
outer_y <- outer_y*deg

plot(outer_x,outer_y, typ="l",
     ylim=c(-90, 90), xlim=c(-180,180),
     xaxt="n", yaxt="n", ann=F)

#-----読み込みデータの座標変換
 l <- (d[[2]] *-1) *rad
 b <- d[[3]] *rad

a <- acos(cos(b)*cos(l/2))
sinc_a <- sin(a) / a

x <- (2*cos(b)*sin(l/2))/sinc_a
y <- sin(b)/sinc_a

aitoff <- data.frame(x *deg, y *deg)
points(aitoff, col="red")