sun_Angle.R

gistfile1.txt

# a <- 6371000
lonlat2xy <- function(lon, lat, a = 6371000) {
    A <- pi/2*a*(1 - lat/90)
    rot <- 0
    x <- cos((lon + rot)*pi/180)*A
    y <- sin((lon + rot)*pi/180)*A
    return(list(x = x, y = y))
}

xy2lonlat <- function(x, y, a = 6371000) {
    A <- sqrt(x^2 + y^2)
    lat <- -90*(A*2/(a*pi) - 1)
    lon <- ConvertLongitude(180/pi*atan2(y, x), from = 180)
    return(list(lon = lon, lat = lat))
}

# ha = altura de la atmósfera
# h = altura del sol
# D = distancia (horizontal) al sol 
# n = índice de refracción de la amtósfera
optical_path <- function(D, h.sun = 4800*1000, n = 1.000277, ha = 10*1000) {
    function(alpha) {
        alpha <- alpha*pi/180
        d1 <- ha/sin(alpha)
        l <- cos(alpha)*d1
        
        # distancia hasta el sol
        d2 <- sqrt((D-l)^2 + (h.sun-ha)^2)
        
        d2 + d1*n
    }
}

sun_angle <- function(lon = 0, lat = 0, lon.sun = 0, lat.sun = 0, 
                      h.sun = 4800*1000, n = 1.000277, ha = 10*1000) {
    # distancia hasta el tope de la atmósfera
    pos <- lonlat2xy(lon, lat)
    dsqr <- pos$x^2 + pos$y^2
    
    pos.sun <- lonlat2xy(lon.sun, lat.sun)
    d.sunsqr <- pos.sun$x^2 + pos.sun$y^2
    
    londif <- abs(lon.sun - lon)
    D <- sqrt(dsqr + d.sunsqr - 2*sqrt(dsqr*d.sunsqr)*cos(londif*pi/180))
    
    alpha1 <- optimise(optical_path(D, h.sun, n, ha),
                        c(0, 90))$minimum
    
    
    alpha2 <- optimise(optical_path(D, h.sun, n = 1, ha),
                       c(0, 90))$minimum
    list(refraction = alpha1, direct = alpha2, difference = alpha1 - alpha2)
}

s <- sun_angle(0, 0)