# R/solarPos.R In solarPos: Solar Position Algorithm for Solar Radiation Applications

```##  Julian day
julianDay <- function(year, month, day, hour=12, min=0, sec=0, tz=0, dut1=0){
day_decimal = day + (hour - tz + (min + (sec + dut1)/60.0)/60.0)/24.0
julian_day <- julian_day+(julian_day > 2299160.0)*B
return(julian_day)
}

##  Solar zenith angle and azimuth angle
solarPosition <- function(jd, lon, lat, delta_t=32.184, elev=0, temp=16, pres=1013.25){
limit_degrees <- function(degrees){
limited <- degrees/360.0
return(360.0*(limited-floor(limited)))
}
d2r <- function(degrees){
return(degrees*pi/180)
}
}
julian_ephemeris_day <- function(jd, delta_t=32.184){
jde <- jd+delta_t/86400.0
return(jde)
}
julian_century <- function(jd){
jc <- (jd-2451545.0)/36525.0
return(jc)
}
julian_ephemeris_century <- function(jde){
jce <- (jde - 2451545.0)/36525.0
return(jce)
}
julian_ephemeris_millennium <- function(jce){
jme <- jce/10.0
return(jme)
}

sunGeocentricPosition <- function(jme){
EPT_L <- list(
c(175347046.0,0,0,
3341656.0,4.6692568,6283.07585,
34894.0,4.6261,12566.1517,
3497.0,2.7441,5753.3849,
3418.0,2.8289,3.5231,
3136.0,3.6277,77713.7715,
2676.0,4.4181,7860.4194,
2343.0,6.1352,3930.2097,
1324.0,0.7425,11506.7698,
1273.0,2.0371,529.691,
1199.0,1.1096,1577.3435,
990,5.233,5884.927,
902,2.045,26.298,
857,3.508,398.149,
780,1.179,5223.694,
753,2.533,5507.553,
505,4.583,18849.228,
492,4.205,775.523,
357,2.92,0.067,
317,5.849,11790.629,
284,1.899,796.298,
271,0.315,10977.079,
243,0.345,5486.778,
206,4.806,2544.314,
205,1.869,5573.143,
202,2.458,6069.777,
156,0.833,213.299,
132,3.411,2942.463,
126,1.083,20.775,
115,0.645,0.98,
103,0.636,4694.003,
102,0.976,15720.839,
102,4.267,7.114,
99,6.21,2146.17,
98,0.68,155.42,
86,5.98,161000.69,
85,1.3,6275.96,
85,3.67,71430.7,
80,1.81,17260.15,
79,3.04,12036.46,
75,1.76,5088.63,
74,3.5,3154.69,
74,4.68,801.82,
70,0.83,9437.76,
62,3.98,8827.39,
61,1.82,7084.9,
57,2.78,6286.6,
56,4.39,14143.5,
56,3.47,6279.55,
52,0.19,12139.55,
52,1.33,1748.02,
51,0.28,5856.48,
49,0.49,1194.45,
41,5.37,8429.24,
41,2.4,19651.05,
39,6.17,10447.39,
37,6.04,10213.29,
37,2.57,1059.38,
36,1.71,2352.87,
36,1.78,6812.77,
33,0.59,17789.85,
30,0.44,83996.85,
30,2.74,1349.87,
25,3.16,4690.48),
c(628331966747.0,0,0,
206059.0,2.678235,6283.07585,
4303.0,2.6351,12566.1517,
425.0,1.59,3.523,
119.0,5.796,26.298,
109.0,2.966,1577.344,
93,2.59,18849.23,
72,1.14,529.69,
68,1.87,398.15,
67,4.41,5507.55,
59,2.89,5223.69,
56,2.17,155.42,
45,0.4,796.3,
36,0.47,775.52,
29,2.65,7.11,
21,5.34,0.98,
19,1.85,5486.78,
19,4.97,213.3,
17,2.99,6275.96,
16,0.03,2544.31,
16,1.43,2146.17,
15,1.21,10977.08,
12,2.83,1748.02,
12,3.26,5088.63,
12,5.27,1194.45,
12,2.08,4694,
11,0.77,553.57,
10,1.3,6286.6,
10,4.24,1349.87,
9,2.7,242.73,
9,5.64,951.72,
8,5.3,2352.87,
6,2.65,9437.76,
6,4.67,4690.48),
c(52919.0,0,0,
8720.0,1.0721,6283.0758,
309.0,0.867,12566.152,
27,0.05,3.52,
16,5.19,26.3,
16,3.68,155.42,
10,0.76,18849.23,
9,2.06,77713.77,
7,0.83,775.52,
5,4.66,1577.34,
4,1.03,7.11,
4,3.44,5573.14,
3,5.14,796.3,
3,6.05,5507.55,
3,1.19,242.73,
3,6.12,529.69,
3,0.31,398.15,
3,2.28,553.57,
2,4.38,5223.69,
2,3.75,0.98),
c(289.0,5.844,6283.076,
35,0,0,
17,5.49,12566.15,
3,5.2,155.42,
1,4.72,3.52,
1,5.3,18849.23,
1,5.97,242.73),
c(114.0,3.142,0,
8,4.13,6283.08,
1,3.84,12566.15),
c(1,3.14,0)
)
EPT_L <- lapply(EPT_L, function(x) matrix(data=x, ncol=3, byrow=TRUE))

EPT_B <- list(
c(280.0,3.199,84334.662,
102.0,5.422,5507.553,
80,3.88,5223.69,
44,3.7,2352.87,
32,4,1577.34),
c(9,3.9,5507.55,
6,1.73,5223.69)
)
EPT_B <- lapply(EPT_B, function(x) matrix(data=x, ncol=3, byrow=TRUE))

EPT_R <- list(
c(100013989.0,0,0,
1670700.0,3.0984635,6283.07585,
13956.0,3.05525,12566.1517,
3084.0,5.1985,77713.7715,
1628.0,1.1739,5753.3849,
1576.0,2.8469,7860.4194,
925.0,5.453,11506.77,
542.0,4.564,3930.21,
472.0,3.661,5884.927,
346.0,0.964,5507.553,
329.0,5.9,5223.694,
307.0,0.299,5573.143,
243.0,4.273,11790.629,
212.0,5.847,1577.344,
186.0,5.022,10977.079,
175.0,3.012,18849.228,
110.0,5.055,5486.778,
98,0.89,6069.78,
86,5.69,15720.84,
86,1.27,161000.69,
65,0.27,17260.15,
63,0.92,529.69,
57,2.01,83996.85,
56,5.24,71430.7,
49,3.25,2544.31,
47,2.58,775.52,
45,5.54,9437.76,
43,6.01,6275.96,
39,5.36,4694,
38,2.39,8827.39,
37,0.83,19651.05,
37,4.9,12139.55,
36,1.67,12036.46,
35,1.84,2942.46,
33,0.24,7084.9,
32,0.18,5088.63,
32,1.78,398.15,
28,1.21,6286.6,
28,1.9,6279.55,
26,4.59,10447.39),
c(103019.0,1.10749,6283.07585,
1721.0,1.0644,12566.1517,
702.0,3.142,0,
32,1.02,18849.23,
31,2.84,5507.55,
25,1.32,5223.69,
18,1.42,1577.34,
10,5.91,10977.08,
9,1.42,6275.96,
9,0.27,5486.78),
c(4359.0,5.7846,6283.0758,
124.0,5.579,12566.152,
12,3.14,0,
9,3.63,77713.77,
6,1.87,5573.14,
3,5.47,18849.23),
c(145.0,4.273,6283.076,
7,3.92,12566.15),
c(4,2.56,6283.08)
)
EPT_R <- lapply(EPT_R, function(x) matrix(data=x, ncol=3, byrow=TRUE))

earth_heliocentric_longitude <- function(jme){
LX <- do.call(rbind,lapply(EPT_L, function(y) rowSums(matrix(data=apply(y,1,function(x) x[1]*cos((x[2]+x[3]*jme))), nrow=length(jme)))))
pow <- apply(matrix(jme, ncol=1), 1, function(x) x^(0:(nrow(LX)-1)))
L <- r2d(colSums(LX*pow)/10^8)
return(limit_degrees(L))
}
earth_heliocentric_latitude <- function(jme){
BX <- do.call(rbind,lapply(EPT_B, function(y) rowSums(matrix(data=apply(y,1,function(x) x[1]*cos((x[2]+x[3]*jme))), nrow=length(jme)))))
pow <- apply(matrix(jme, ncol=1), 1, function(x) x^(0:(nrow(BX)-1)))
B <- r2d(colSums(BX*pow)/10^8)
return(B)
}
RX <- do.call(rbind,lapply(EPT_R, function(y) rowSums(matrix(data=apply(y,1,function(x) x[1]*cos((x[2]+x[3]*jme))), nrow=length(jme)))))
pow <- apply(matrix(jme, ncol=1), 1, function(x) x^(0:(nrow(RX)-1)))
R <- colSums(RX*pow)/10^8
return(R)
}
sun_geocentric_longitude <- function(L){
Theta <- L+180
return(limit_degrees(Theta))
}
sun_geocentric_latitude <- function(B){
beta <- -B
return(limit_degrees(beta))
}

L <- earth_heliocentric_longitude(jme)
B <- earth_heliocentric_latitude(jme)

Theta <- sun_geocentric_longitude(L)
beta <- sun_geocentric_latitude(B)

return(matrix(data=c(Theta, beta, R), ncol=3))
}
nutation <- function(jce){
Y_TERMS <- c(
0,0,0,0,1,
-2,0,0,2,2,
0,0,0,2,2,
0,0,0,0,2,
0,1,0,0,0,
0,0,1,0,0,
-2,1,0,2,2,
0,0,0,2,1,
0,0,1,2,2,
-2,-1,0,2,2,
-2,0,1,0,0,
-2,0,0,2,1,
0,0,-1,2,2,
2,0,0,0,0,
0,0,1,0,1,
2,0,-1,2,2,
0,0,-1,0,1,
0,0,1,2,1,
-2,0,2,0,0,
0,0,-2,2,1,
2,0,0,2,2,
0,0,2,2,2,
0,0,2,0,0,
-2,0,1,2,2,
0,0,0,2,0,
-2,0,0,2,0,
0,0,-1,2,1,
0,2,0,0,0,
2,0,-1,0,1,
-2,2,0,2,2,
0,1,0,0,1,
-2,0,1,0,1,
0,-1,0,0,1,
0,0,2,-2,0,
2,0,-1,2,1,
2,0,1,2,2,
0,1,0,2,2,
-2,1,1,0,0,
0,-1,0,2,2,
2,0,0,2,1,
2,0,1,0,0,
-2,0,2,2,2,
-2,0,1,2,1,
2,0,-2,0,1,
2,0,0,0,1,
0,-1,1,0,0,
-2,-1,0,2,1,
-2,0,0,0,1,
0,0,2,2,1,
-2,0,2,0,1,
-2,1,0,2,1,
0,0,1,-2,0,
-1,0,1,0,0,
-2,1,0,0,0,
1,0,0,0,0,
0,0,1,2,0,
0,0,-2,2,2,
-1,-1,1,0,0,
0,1,1,0,0,
0,-1,1,2,2,
2,-1,-1,2,2,
0,0,3,2,2,
2,-1,0,2,2)
Y_TERMS <- matrix(data=Y_TERMS, ncol=5, byrow=TRUE)

PE_TERMS <- c(
-171996,-174.2,92025,8.9,
-13187,-1.6,5736,-3.1,
-2274,-0.2,977,-0.5,
2062,0.2,-895,0.5,
1426,-3.4,54,-0.1,
712,0.1,-7,0,
-517,1.2,224,-0.6,
-386,-0.4,200,0,
-301,0,129,-0.1,
217,-0.5,-95,0.3,
-158,0,0,0,
129,0.1,-70,0,
123,0,-53,0,
63,0,0,0,
63,0.1,-33,0,
-59,0,26,0,
-58,-0.1,32,0,
-51,0,27,0,
48,0,0,0,
46,0,-24,0,
-38,0,16,0,
-31,0,13,0,
29,0,0,0,
29,0,-12,0,
26,0,0,0,
-22,0,0,0,
21,0,-10,0,
17,-0.1,0,0,
16,0,-8,0,
-16,0.1,7,0,
-15,0,9,0,
-13,0,7,0,
-12,0,6,0,
11,0,0,0,
-10,0,5,0,
-8,0,3,0,
7,0,-3,0,
-7,0,0,0,
-7,0,3,0,
-7,0,3,0,
6,0,0,0,
6,0,-3,0,
6,0,-3,0,
-6,0,3,0,
-6,0,3,0,
5,0,0,0,
-5,0,3,0,
-5,0,3,0,
-5,0,3,0,
4,0,0,0,
4,0,0,0,
4,0,0,0,
-4,0,0,0,
-4,0,0,0,
-4,0,0,0,
3,0,0,0,
-3,0,0,0,
-3,0,0,0,
-3,0,0,0,
-3,0,0,0,
-3,0,0,0,
-3,0,0,0,
-3,0,0,0)
PE_TERMS <- matrix(data=PE_TERMS, ncol=4, byrow=TRUE)

mean_elongation_moon_sun <- function(jce){
X0 <- 297.85036 + 445267.11148*jce - 0.0019142*jce^2 + (1.0/189474.0)*jce^3
return(X0)
}
mean_anomaly_sun <- function(jce){
X1 <- 357.52772 + 35999.05034*jce - 0.0001603*jce^2 - (1.0/300000.0)*jce^3
return(X1)
}
mean_anomaly_moon <- function(jce){
X2 <- 134.96298 + 477198.867398*jce + 0.0086972*jce^2 + (1.0/56250.0)*jce^3
return(X2)
}
argument_latitude_moon <- function(jce){
X3 <- 93.27191 + 483202.017538*jce - 0.0036825*jce^2 +(1.0/327270.0)*jce^3
return(X3)
}
ascending_longitude_moon <- function(jce){
X4 <- 125.04452 - 1934.136261*jce + 0.0020708*jce^2 + (1.0/450000.0)*jce^3
return(X4)
}
nutation_longitude <- function(jce, X0, X1, X2, X3, X4){
X <- c(X0, X1, X2, X3, X4)
f1 <- apply(PE_TERMS[,1:2], 1, function(x) x[1]+x[2]*jce)
f2 <- apply(Y_TERMS, 1, function(Y) sin(sum(X*Y)*pi/180))
delta_psi <- matrix(data=f1, nrow=length(jce)) %*% matrix(data=f2,ncol=1)/36000000
return(as.numeric(delta_psi))
}
nutation_obliquity <- function(jce, X0, X1, X2, X3, X4){
X <- c(X0, X1, X2, X3, X4)
f1 <- apply(PE_TERMS[,3:4], 1, function(x) x[1]+x[2]*jce)
f2 <- apply(Y_TERMS, 1, function(Y) cos(sum(X*Y)*pi/180))
delta_epsilon <- matrix(data=f1, nrow=length(jce)) %*% matrix(data=f2,ncol=1)/36000000
return(as.numeric(delta_epsilon))
}

X0 <- mean_elongation_moon_sun(jce)
X1 <- mean_anomaly_sun(jce)
X2 <- mean_anomaly_moon(jce)
X3 <- argument_latitude_moon(jce)
X4 <- ascending_longitude_moon(jce)
delta_psi <- nutation_longitude(jce, X0, X1, X2, X3, X4)
delta_epsilon <- nutation_obliquity(jce, X0, X1, X2, X3, X4)

return(matrix(data=c(delta_psi,delta_epsilon), ncol=2))
}

ecliptic_mean_obliquity <- function(jme, delta_epsilon){
u <- jme/10.0
epsilon0 <- 84381.448 + -4680.93*u + -1.55*u^2 + 1999.25*u^3 -51.38*u^4 -249.67*u^5 -
39.05*u^6 + 7.12*u^7 + 27.87*u^8 + 5.79*u^9 + 2.45*u^10
epsilon <- epsilon0/3600 + delta_epsilon
return(epsilon)
}
aberration_correction <- function(R){
delta_tau <- -20.4898 / (3600.0*R)
return (delta_tau)
}
apparent_sun_longitude <- function(Theta, delta_psi, delta_tau){
lambda <- Theta + delta_psi + delta_tau
return(lambda)
}
greenwich_sidereal_time <- function(jd, jc, delta_psi, epsilon){
nu0 <- 280.46061837 + 360.98564736629 * (jd - 2451545.0) + 0.000387933*jc^2 - jc^3/38710000.0
nu0 <- limit_degrees(nu0)
nu <- nu0 + delta_psi*cos(d2r(epsilon))
return(nu)
}
geocentric_right_ascension <- function(lamda, epsilon, beta){
return(limit_degrees(alpha))
}
geocentric_declination <- function(beta, epsilon, lamda){
return(delta)
}
observer_hour_angle <- function(nu, longitude, alpha){
H <- nu + longitude - alpha
H <- limit_degrees(H)
return(H)
}
sun_equatorial_horizontal_parallax <- function(R){
xi <- 8.794/(3600*R)
return(xi)
}
parallax_right_ascension <- function(latitude, elevation, xi, H, delta){
y <- 0.99664719 * sin(u) + elevation*sin(lat_rad)/6378140.0
return(list(r2d(delta_alpha),x,y))
}
topocentric_right_ascension <- function(alpha, delta_alpha){
alpha_prime <- alpha + delta_alpha
return(alpha_prime)
}
topocentric_declination <- function(xi, H, delta, delta_alpha, x, y){
return(r2d(delta_prime))
}
topocentric_local_hour_angle <- function(H, delta_alpha){
H_prime <- H - delta_alpha
return(H_prime)
}
topocentric_elevation_angle <- function(latitude, delta_prime, H_prime){
return(r2d(e0))
}
atmospheric_refraction_correction <- function(pressure, temperature, e0){
delta_e <- (pressure/1010.0)*(283.0/(273.0 + temperature))*
1.02/(60.0*tan(d2r(e0 + 10.3/(e0 + 5.11))))
#if (e0 >= -1*(SUN_RADIUS + atmos_refract))
return(delta_e)
}
topocentric_elevation_angle_corrected <- function(e0, delta_e){
e <- e0 + delta_e
return(e)
}
topocentric_zenith_angle <- function(e){
theta <- 90-e
return(theta)
}
topocentric_azimuth_angle <- function(H_prime, latitude, delta_prime){
return(limit_degrees(Phi))
}

jde <- julian_ephemeris_day(jd,delta_t=delta_t)
jc <- julian_century(jd)
jce <- julian_ephemeris_century(jde)
jme <- julian_ephemeris_millennium(jce)
sgp <- sunGeocentricPosition(jme)
nut <- nutation(jce)
epsilon <- ecliptic_mean_obliquity(jme, nut[,2])
delta_tau <- aberration_correction(sgp[,3])
lambda <- apparent_sun_longitude(sgp[,1], nut[,1], delta_tau)
nu <- greenwich_sidereal_time(jd, jc, nut[,1], epsilon)
alpha <- geocentric_right_ascension(lambda, epsilon, sgp[,2])
delta <- geocentric_declination(sgp[,2], epsilon, lambda)
H <- observer_hour_angle(nu, lon, alpha)
xi <- sun_equatorial_horizontal_parallax(sgp[,3])
PRA <- parallax_right_ascension(lat, elev, xi, H, delta)
delta_alpha <- PRA[[1]]
x <- PRA[[2]]
y <- PRA[[3]]
alhpa_prime <- topocentric_right_ascension(alpha, delta_alpha)
delta_prime <- topocentric_declination(xi, H, delta, delta_alpha,x,y)
H_prime <- topocentric_local_hour_angle(H, delta_alpha)
e0 <- topocentric_elevation_angle(lat, delta_prime, H_prime)
delta_e <- atmospheric_refraction_correction(pres, temp, e0)
e <- topocentric_elevation_angle_corrected(e0, delta_e)
theta <- topocentric_zenith_angle(e)
Phi <- topocentric_azimuth_angle(H_prime, lat, delta_prime)

return(matrix(data=c(theta,Phi),ncol=2, dimnames=list(NULL, c("zenith","azimuth"))))
}
```

## Try the solarPos package in your browser

Any scripts or data that you put into this service are public.

solarPos documentation built on May 2, 2019, 9:59 a.m.