Browse Source

Didn't notice the original R port had lat/lon mixed up

master
boB Rudis 6 years ago
parent
commit
f84581f012
No known key found for this signature in database GPG Key ID: 1D7529BE14E2BBA9
  1. 2
      README.Rmd
  2. 50
      src/terminator-main.cpp

2
README.Rmd

@ -67,7 +67,7 @@ term_seq <- terminator_lat_lon()
chart <- ggplot(term_seq, aes(frame = frame)) + chart <- ggplot(term_seq, aes(frame = frame)) +
borders("world", colour = "gray90", fill = "gray85") + borders("world", colour = "gray90", fill = "gray85") +
geom_ribbon(aes(lat, ymax = lon), ymin = 90, alpha = 0.2) + geom_ribbon(aes(lon, ymax = lat), ymin = 90, alpha = 0.2) +
coord_equal(xlim = c(-180, 190), ylim = c(-58, 85), expand = 0) + coord_equal(xlim = c(-180, 190), ylim = c(-58, 85), expand = 0) +
ggthemes::theme_map() ggthemes::theme_map()

50
src/terminator-main.cpp

@ -44,11 +44,11 @@ inline double ecliptic_obliquity(double j_day) {
// compute epsilon // compute epsilon
return(23.43929111 - return(23.43929111 -
T * (46.836769 / 3600.0 T * (46.836769 / 3600.0
- T * (0.0001831 / 3600.0 - T * (0.0001831 / 3600.0
+ T * (0.00200340 / 3600.0 + T * (0.00200340 / 3600.0
- T * (0.576e-6 / 3600.0 - T * (0.576e-6 / 3600.0
- T * 4.34e-8 / 3600.0))))); - T * 4.34e-8 / 3600.0)))));
} }
@ -97,40 +97,40 @@ DataFrame terminator(int time, double from = -180, double to = 180, double by =
// calculate latitude and longitude of terminator within specified range using time (in POSIXct format, e.g. `Sys.time()`) // calculate latitude and longitude of terminator within specified range using time (in POSIXct format, e.g. `Sys.time()`)
double j_day = get_julian(time); double j_day = get_julian(time);
double gst = get_gmst(j_day); double gst = get_gmst(j_day);
std::vector< double > sunEclPos = sun_ecliptic_position(j_day); std::vector< double > sunEclPos = sun_ecliptic_position(j_day);
double eclObliq = ecliptic_obliquity(j_day); double eclObliq = ecliptic_obliquity(j_day);
std::vector< double > sunEqPos = sun_equatorial_position(sunEclPos[0], eclObliq); std::vector< double > sunEqPos = sun_equatorial_position(sunEclPos[0], eclObliq);
std::vector< double > out_lat, out_lon; std::vector< double > out_lat, out_lon;
out_lat.reserve(4000); out_lat.reserve(4000);
out_lon.reserve(4000); out_lon.reserve(4000);
int n=0; int n=0;
for (double i=from; i<=to; i+=by) { for (double i=from; i<=to; i+=by) {
n += 1; n += 1;
out_lat.push_back(i); out_lon.push_back(i);
out_lon.push_back( out_lat.push_back(
longitude( longitude(
hour_angle(i, sunEqPos, gst), hour_angle(i, sunEqPos, gst),
sunEqPos sunEqPos
) )
); );
} }
out_lat.resize(n); out_lat.resize(n);
out_lon.resize(n); out_lon.resize(n);
return( return(
DataFrame::create( DataFrame::create(
Named("lat") = out_lat, Named("lat") = out_lat,
Named("lon") = out_lon Named("lon") = out_lon
) )
); );
} }
Loading…
Cancel
Save