fly/flightgear
1
0
Fork 0
Code Activity

Merge branch 'next' of gitorious.org:fg/flightgear into next

Browse source
This commit is contained in:
Torsten Dreyer 2010-12-05 13:06:56 +01:00
commit c467c91c07
4 changed files with 402 additions and 197 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

459
src/Instrumentation/dclgps.cxx
View file

@ -26,8 +26,10 @@
#include "dclgps.hxx" #include "dclgps.hxx"
#include <simgear/sg_inlines.h> #include <simgear/sg_inlines.h>
#include <simgear/misc/sg_path.hxx>
#include <simgear/timing/sg_time.hxx> #include <simgear/timing/sg_time.hxx>
#include <simgear/magvar/magvar.hxx> #include <simgear/magvar/magvar.hxx>
#include <simgear/structure/exception.hxx>
#include <Main/fg_props.hxx> #include <Main/fg_props.hxx>
#include <Navaids/fix.hxx> #include <Navaids/fix.hxx>
@ -35,6 +37,7 @@
#include <Airports/simple.hxx> #include <Airports/simple.hxx>
#include <Airports/runways.hxx> #include <Airports/runways.hxx>
#include <fstream>
#include <iostream> #include <iostream>
using namespace std; using namespace std;
@ -240,95 +243,7 @@ void DCLGPS::init() {
// Not sure if this should be here, but OK for now. // Not sure if this should be here, but OK for now.
CreateDefaultFlightPlans(); CreateDefaultFlightPlans();
// Hack - hardwire some instrument approaches for development. LoadApproachData();
// These will shortly be replaced by a routine to read ARINC data from file instead.
FGNPIAP* iap;
GPSWaypoint* wp;
GPSFlightPlan* fp;
const GPSWaypoint* cwp;
iap = new FGNPIAP;
iap->_aptIdent = "KHAF";
iap->_ident = "R12-Y";
iap->_name = ExpandSIAPIdent(iap->_ident);
iap->_rwyStr = "12";
iap->_approachRoutes.clear();
iap->_IAP.clear();
// -------
wp = new GPSWaypoint;
wp->id = "GOBBS";
// Nasty using the find any function here, but it saves converting data from FGFix etc.
cwp = FindFirstByExactId(wp->id);
if(cwp) {
*wp = *cwp;
wp->appType = GPS_IAF;
fp = new GPSFlightPlan;
fp->waypoints.push_back(wp);
} else {
//cout << "Unable to find waypoint " << wp->id << '\n';
}
// -------
wp = new GPSWaypoint;
wp->id = "FUJCE";
cwp = FindFirstByExactId(wp->id);
if(cwp) {
*wp = *cwp;
wp->appType = GPS_IAP;
fp->waypoints.push_back(wp);
iap->_approachRoutes.push_back(fp);
iap->_IAP.push_back(wp);
} else {
//cout << "Unable to find waypoint " << wp->id << '\n';
}
// -------
wp = new GPSWaypoint;
wp->id = "JEVXY";
cwp = FindFirstByExactId(wp->id);
if(cwp) {
*wp = *cwp;
wp->appType = GPS_FAF;
iap->_IAP.push_back(wp);
} else {
//cout << "Unable to find waypoint " << wp->id << '\n';
}
// -------
wp = new GPSWaypoint;
wp->id = "RW12";
wp->appType = GPS_MAP;
if(wp->id.substr(0, 2) == "RW" && wp->appType == GPS_MAP) {
// Assume that this is a missed-approach point based on the runway number, which appears to be standard for most approaches.
const FGAirport* apt = fgFindAirportID(iap->_aptIdent);
if(apt) {
// TODO - sanity check the waypoint ID to ensure we have a double digit number
FGRunway* rwy = apt->getRunwayByIdent(wp->id.substr(2, 2));
if(rwy) {
wp->lat = rwy->begin().getLatitudeRad();
wp->lon = rwy->begin().getLongitudeRad();
}
}
} else {
cwp = FindFirstByExactId(wp->id);
if(cwp) {
*wp = *cwp;
wp->appType = GPS_MAP;
} else {
//cout << "Unable to find waypoint " << wp->id << '\n';
}
}
iap->_IAP.push_back(wp);
// -------
wp = new GPSWaypoint;
wp->id = "SEEMS";
cwp = FindFirstByExactId(wp->id);
if(cwp) {
*wp = *cwp;
wp->appType = GPS_MAHP;
iap->_IAP.push_back(wp);
} else {
//cout << "Unable to find waypoint " << wp->id << '\n';
}
// -------
_np_iap[iap->_aptIdent].push_back(iap);
} }
void DCLGPS::bind() { void DCLGPS::bind() {
@ -685,6 +600,371 @@ string DCLGPS::ExpandSIAPIdent(const string& ident) {
return(name); return(name);
} }
/*
Load instrument approaches from an ARINC 424 file.
Tested on ARINC 424-18.
Known / current best guess at the format:
Col 1: Always 'S'. If it isn't, ditch it.
Col 2-4: "Customer area" code, eg "USA", "CAN". I think that CAN is used for Alaska.
Col 5: Section code. Used in conjunction with sub-section code. Definitions are with sub-section code.
Col 6: Always blank.
Col 7-10: ICAO (or FAA) airport ident. Left justified if < 4 chars.
Col 11-12: Based on ICAO geographical region.
Col 13: Sub-section code. Used in conjunction with section code.
"HD/E/F" => Helicopter record.
"HS" => Helicopter minimum safe altitude.
"PA" => Airport record.
"PF" => Approach segment.
"PG" => Runway record.
"PP" => Path point record. ???
"PS" => MSA record (minimum safe altitude).
------ The following is for "PF", approach segment -------
Col 14-19: SIAP ident for this approach (left justified). This is a standardised abbreviated approach name.
e.g. "R10LZ" expands to "RNAV (GPS) Z RWY 10 L". See the comment block for ExpandSIAPIdent for full details.
Col 20: Route type. This is tricky - I don't have full documentation and am having to guess a bit.
'A' => Arrival route? This seems to be used to encode arrival routes from the IAF to the approach proper.
Note that the final fix of the arrival route is duplicated in the approach proper.
'D' => VOR/DME or GPS
'N' => NDB or GPS
'P' => GPS (ARINC 424-18), GPS and RNAV (GPS) (ARINC 424-15 and before).
'R' => RNAV (GPS) (ARINC 424-18).
'S' => VOR or GPS
Col 21-25: Transition identifier. AFAICT, this is the ident of the IAF for this initial approach route, and left blank for the final approach course. See col 30-34 for the actual fix ident.
Col 26: BLANK
Col 27-29: Sequence number - position within the route segment. Rule: 10-20-30 etc.
Col 30-34: Fix identifer. The ident of the waypoint.
Col 35-36: ICAO geographical region code. I think we can ignore this for now.
Col 37: Section code - ??? I don't know what this means
Col 38 Subsection code - ??? ditto - no idea!
Col 40: Waypoint type.
'A' => Airport as waypoint
'E' => Essential waypoint (e.g. change of heading at this waypoint, etc).
'G' => Runway or helipad as waypoint
'H' => Heliport as waypoint
'N' => NDB as waypoint
'P' => Phantom waypoint (not sure if this is used in rev 18?)
'V' => VOR as waypoint
Col 41: Waypoint type.
'B' => Flyover, approach transition, or final approach.
'E' => end of route segment (transition waypoint). (Actually "End of terminal procedure route type" in the docs).
'N' => ??? I've also seen 'N' in this column, but don't know what it indicates.
'Y' => Flyover.
Col 43: Waypoint type. May also be blank when none of the below.
'A' => Initial approach fix (IAF)
'F' => Final approach fix
'H' => Holding fix
'I' => Final approach course fix
'M' => Missed approach point
'P' => ??? This is present, but I don't know what this means and it wasn't in the FAA docs that I found the above in!
??? Possibly procedure turn?
'C' => ??? This is also present in the data, but missing from the docs. Is at airport 00R.
Col 107-111 MSA center fix. We can ignore this.
*/
void DCLGPS::LoadApproachData() {
FGNPIAP* iap;
GPSWaypoint* wp;
GPSFlightPlan* fp;
const GPSWaypoint* cwp;
std::ifstream fin;
SGPath path = globals->get_fg_root();
path.append("Navaids/rnav.dat");
fin.open(path.c_str(), ios::in);
if(!fin) {
//cout << "Unable to open input file " << path.c_str() << '\n';
return;
} else {
//cout << "Opened " << path.c_str() << " for reading\n";
}
char tmp[256];
string s;
string apt_ident; // This gets set to the ICAO code of the current airport being processed.
string iap_ident; // The abbreviated name of the current approach being processed.
string wp_ident; // The ident of the waypoint of the current line
string last_apt_ident;
string last_iap_ident;
string last_wp_ident;
// There is no need to save the full name - it can be generated on the fly from the abbreviated name as and when needed.
bool apt_in_progress = false; // Set true whilst loading all the approaches for a given airport.
bool iap_in_progress = false; // Set true whilst loading a given approach.
bool iap_error = false; // Set true if there is an error loading a given approach.
bool route_in_progress = false; // Set true when we are loading a "route" segment of the approach.
int last_sequence_number = 0; // Position within the route, rule (rev 18): 10, 20, 30 etc.
int sequence_number;
char last_route_type = 0;
char route_type;
char waypoint_fix_type; // This is the waypoint type from col 43, i.e. the type of fix. May be blank.
int j;
// Debugging info
unsigned int nLoaded = 0;
unsigned int nErrors = 0;
//for(i=0; i<64; ++i) {
while(!fin.eof()) {
fin.getline(tmp, 256);
//s = Fake_rnav_dat[i];
s = tmp;
if(s.size() < 132) continue;
if(s[0] == 'S') { // Valid line
string country_code = s.substr(1, 3);
if(country_code == "USA") { // For now we'll stick to US procedures in case there are unknown gotchas with others
if(s[4] == 'P') { // Includes approaches.
if(s[12] == 'A') { // Airport record
apt_ident = s.substr(6, 4);
// Trim any whitespace from the ident. The ident is left justified,
// so any space will be at the end.
if(apt_ident[3] == ' ') apt_ident = apt_ident.substr(0, 3);
// I think that all idents are either 3 or 4 chars - could check this though!
if(!apt_in_progress) {
last_apt_ident = apt_ident;
apt_in_progress = 1;
} else {
if(last_apt_ident != apt_ident) {
if(iap_in_progress) {
if(iap_error) {
//cout << "ERROR: Unable to load approach " << iap->_ident << " at " << iap->_aptIdent << '\n';
nErrors++;
} else {
_np_iap[iap->_aptIdent].push_back(iap);
//cout << "** Loaded " << iap->_aptIdent << "\t" << iap->_ident << '\n';
nLoaded++;
}
iap_in_progress = false;
}
}
last_apt_ident = apt_ident;
}
iap_in_progress = 0;
} else if(s[12] == 'F') { // Approach segment
if(apt_in_progress) {
iap_ident = s.substr(13, 6);
// Trim any whitespace from the RH end.
for(j=0; j<6; ++j) {
if(iap_ident[5-j] == ' ') {
iap_ident = iap_ident.substr(0, 5-j);
} else {
// It's important to break here, since earlier versions of ARINC 424 allowed spaces in the ident.
break;
}
}
if(iap_in_progress) {
if(iap_ident != last_iap_ident) {
// This is a new approach - store the last one and trigger
// starting afresh by setting the in progress flag to false.
if(iap_error) {
//cout << "ERROR: Unable to load approach " << iap->_ident << " at " << iap->_aptIdent << '\n';
nErrors++;
} else {
_np_iap[iap->_aptIdent].push_back(iap);
//cout << "Loaded " << iap->_aptIdent << "\t" << iap->_ident << '\n';
nLoaded++;
}
iap_in_progress = false;
}
}
if(!iap_in_progress) {
iap = new FGNPIAP;
iap->_aptIdent = apt_ident;
iap->_ident = iap_ident;
iap->_name = ExpandSIAPIdent(iap_ident); // I suspect that it's probably better to just store idents, and to expand the names as needed.
// Note, we haven't set iap->_rwyStr yet.
last_iap_ident = iap_ident;
iap_in_progress = true;
iap_error = false;
}
// Route type
route_type = s[19];
sequence_number = atoi(s.substr(26,3).c_str());
wp_ident = s.substr(29, 5);
waypoint_fix_type = s[42];
// Trim any whitespace from the RH end
for(j=0; j<5; ++j) {
if(wp_ident[4-j] == ' ') {
wp_ident = wp_ident.substr(0, 4-j);
} else {
break;
}
}
// Ignore lines with no waypoint ID for now - these are normally part of the
// missed approach procedure, and we don't use them in the KLN89.
if(!wp_ident.empty()) {
// Make a local copy of the waypoint for now, since we're not yet sure if we'll be using it
GPSWaypoint w;
w.id = wp_ident;
bool wp_error = false;
if(w.id.substr(0, 2) == "RW" && waypoint_fix_type == 'M') {
// Assume that this is a missed-approach point based on the runway number, which appears to be standard for most approaches.
// Note: Currently fgFindAirportID returns NULL on error, but getRunwayByIdent throws an exception.
const FGAirport* apt = fgFindAirportID(iap->_aptIdent);
if(apt) {
string rwystr;
try {
rwystr = w.id.substr(2, 2);
// TODO - sanity check the rwystr at this point to ensure we have a double digit number
if(w.id.size() > 4) {
if(w.id[4] == 'L' || w.id[4] == 'C' || w.id[4] == 'R') {
rwystr += w.id[4];
}
}
FGRunway* rwy = apt->getRunwayByIdent(rwystr);
w.lat = rwy->begin().getLatitudeRad();
w.lon = rwy->begin().getLongitudeRad();
} catch(const sg_exception&) {
SG_LOG(SG_GENERAL, SG_WARN, "Unable to find runway " << w.id.substr(2, 2) << " at airport " << iap->_aptIdent);
//cout << "Unable to find runway " << w.id.substr(2, 2) << " at airport " << iap->_aptIdent << " ( w.id = " << w.id << ", rwystr = " << rwystr << " )\n";
wp_error = true;
}
} else {
wp_error = true;
}
} else {
cwp = FindFirstByExactId(w.id);
if(cwp) {
w = *cwp;
} else {
wp_error = true;
}
}
switch(waypoint_fix_type) {
case 'A': w.appType = GPS_IAF; break;
case 'F': w.appType = GPS_FAF; break;
case 'H': w.appType = GPS_MAHP; break;
case 'I': w.appType = GPS_IAP; break;
case 'M': w.appType = GPS_MAP; break;
case ' ': w.appType = GPS_APP_NONE; break;
//default: cout << "Unknown waypoint_fix_type: \'" << waypoint_fix_type << "\' [" << apt_ident << ", " << iap_ident << "]\n";
}
if(wp_error) {
//cout << "Unable to find waypoint " << w.id << " [" << apt_ident << ", " << iap_ident << "]\n";
iap_error = true;
}
if(!wp_error) {
if(route_in_progress) {
if(sequence_number > last_sequence_number) {
// TODO - add a check for runway numbers
// Check for the waypoint ID being the same as the previous line.
// This is often the case for the missed approach holding point.
if(wp_ident == last_wp_ident) {
if(waypoint_fix_type == 'H') {
if(!iap->_IAP.empty()) {
if(iap->_IAP[iap->_IAP.size() - 1]->appType == GPS_APP_NONE) {
iap->_IAP[iap->_IAP.size() - 1]->appType = GPS_MAHP;
} else {
//cout << "Waypoint is MAHP and another type! " << w.id << " [" << apt_ident << ", " << iap_ident << "]\n";
}
}
}
} else {
// Create a new waypoint on the heap, copy the local copy into it, and push it onto the approach.
wp = new GPSWaypoint;
*wp = w;
if(route_type == 'A') {
fp->waypoints.push_back(wp);
} else {
iap->_IAP.push_back(wp);
}
}
} else if(sequence_number == last_sequence_number) {
// This seems to happen once per final approach route - one of the waypoints
// is duplicated with the same sequence number - I'm not sure what information
// the second line give yet so ignore it for now.
// TODO - figure this out!
} else {
// Finalise the current route and start a new one
//
// Finalise the current route
if(last_route_type == 'A') {
// Push the flightplan onto the approach
iap->_approachRoutes.push_back(fp);
} else {
// All the waypoints get pushed individually - don't need to do it.
}
// Start a new one
// There are basically 2 possibilities here - either it's one of the arrival transitions,
// or it's the core final approach course.
wp = new GPSWaypoint;
*wp = w;
if(route_type == 'A') { // It's one of the arrival transition(s)
fp = new GPSFlightPlan;
fp->waypoints.push_back(wp);
} else {
iap->_IAP.push_back(wp);
}
route_in_progress = true;
}
} else {
// Start a new route.
// There are basically 2 possibilities here - either it's one of the arrival transitions,
// or it's the core final approach course.
wp = new GPSWaypoint;
*wp = w;
if(route_type == 'A') { // It's one of the arrival transition(s)
fp = new GPSFlightPlan;
fp->waypoints.push_back(wp);
} else {
iap->_IAP.push_back(wp);
}
route_in_progress = true;
}
last_route_type = route_type;
last_wp_ident = wp_ident;
last_sequence_number = sequence_number;
}
}
} else {
// ERROR - no airport record read.
}
}
} else {
// Check and finalise any approaches in progress
// TODO - sanity check that the approach has all the required elements
if(iap_in_progress) {
// This is a new approach - store the last one and trigger
// starting afresh by setting the in progress flag to false.
if(iap_error) {
//cout << "ERROR: Unable to load approach " << iap->_ident << " at " << iap->_aptIdent << '\n';
nErrors++;
} else {
_np_iap[iap->_aptIdent].push_back(iap);
//cout << "* Loaded " << iap->_aptIdent << "\t" << iap->_ident << '\n';
nLoaded++;
}
iap_in_progress = false;
}
}
}
}
}
// If we get to the end of the file, load any approach that is still in progress
// TODO - sanity check that the approach has all the required elements
if(iap_in_progress) {
if(iap_error) {
//cout << "ERROR: Unable to load approach " << iap->_ident << " at " << iap->_aptIdent << '\n';
nErrors++;
} else {
_np_iap[iap->_aptIdent].push_back(iap);
//cout << "*** Loaded " << iap->_aptIdent << "\t" << iap->_ident << '\n';
nLoaded++;
}
}
//cout << "Done loading approach database\n";
//cout << "Loaded: " << nLoaded << '\n';
//cout << "Failed: " << nErrors << '\n';
fin.close();
}
GPSWaypoint* DCLGPS::GetActiveWaypoint() { GPSWaypoint* DCLGPS::GetActiveWaypoint() {
return &_activeWaypoint; return &_activeWaypoint;
} }
@ -735,6 +1015,7 @@ void DCLGPS::DtoInitiate(const string& s) {
_fromWaypoint.id = "DTOWP"; _fromWaypoint.id = "DTOWP";
delete wp; delete wp;
} else { } else {
// TODO - Should bring up the user waypoint page.
_dto = false; _dto = false;
} }
} }

2
src/Instrumentation/dclgps.hxx
View file

@ -342,6 +342,8 @@ protected:
protected: protected:
void LoadApproachData();
// Find first of any type of waypoint by id. (TODO - Possibly we should return multiple waypoints here). // Find first of any type of waypoint by id. (TODO - Possibly we should return multiple waypoints here).
GPSWaypoint* FindFirstById(const string& id) const; GPSWaypoint* FindFirstById(const string& id) const;
GPSWaypoint* FindFirstByExactId(const string& id) const; GPSWaypoint* FindFirstByExactId(const string& id) const;

115
src/Time/light.cxx
View file

@ -48,38 +48,6 @@
#include "light.hxx" #include "light.hxx"
#include "sunsolver.hxx" #include "sunsolver.hxx"
/**
* Map i.e. project a vector onto a plane.
* @param normal (in) normal vector for the plane
* @param v0 (in) a point on the plane
* @param vec (in) the vector to map onto the plane
*/
static SGVec3f map_vec_onto_cur_surface_plane(const SGVec3f& normal,
const SGVec3f& v0,
const SGVec3f& vec)
{
// calculate a vector "u1" representing the shortest distance from
// the plane specified by normal and v0 to a point specified by
// "vec". "u1" represents both the direction and magnitude of
// this desired distance.
// u1 = ( (normal <dot> vec) / (normal <dot> normal) ) * normal
SGVec3f u1 = (dot(normal, vec) / dot(normal, normal)) * normal;
// calculate the vector "v" which is the vector "vec" mapped onto
// the plane specified by "normal" and "v0".
// v = v0 + vec - u1
SGVec3f v = v0 + vec - u1;
// Calculate the vector "result" which is "v" - "v0" which is a
// directional vector pointing from v0 towards v
// result = v - v0
return v - v0;
}
// Constructor // Constructor
FGLight::FGLight () FGLight::FGLight ()
: _ambient_tbl( NULL ), : _ambient_tbl( NULL ),
@ -432,16 +400,22 @@ void FGLight::updateSunPos()
SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t->getGst() ); SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t->getGst() );
double sun_l; double sun_l;
double sun_gd_lat; double sun_gc_lat;
fgSunPositionGST(t->getGst(), &sun_l, &sun_gd_lat); fgSunPositionGST(t->getGst(), &sun_l, &sun_gc_lat);
set_sun_lon(sun_l); set_sun_lon(sun_l);
set_sun_lat(sun_gd_lat); // It might seem that sun_gc_lat needs to be converted to geodetic
SGVec3d sunpos(SGVec3d::fromGeod(SGGeod::fromRad(sun_l, sun_gd_lat))); // latitude here, but it doesn't. The sun latitude is the latitude
// of the point on the earth where the up vector has the same
// angle from geocentric Z as the sun direction. But geodetic
// latitude is defined as 90 - angle of up vector from Z!
set_sun_lat(sun_gc_lat);
SGVec3d sunpos(SGVec3d::fromGeoc(SGGeoc::fromRadM(sun_l, sun_gc_lat,
SGGeodesy::EQURAD)));
SG_LOG( SG_EVENT, SG_DEBUG, " t->cur_time = " << t->get_cur_time() ); SG_LOG( SG_EVENT, SG_DEBUG, " t->cur_time = " << t->get_cur_time() );
SG_LOG( SG_EVENT, SG_DEBUG, SG_LOG( SG_EVENT, SG_DEBUG,
" Sun Geodetic lat = " << sun_gd_lat " Sun Geocentric lat = " << sun_gc_lat
<< " Geodetic lat = " << sun_gd_lat ); << " Geodcentric lat = " << sun_gc_lat );
// update the sun light vector // update the sun light vector
sun_vec() = SGVec4f(toVec3f(normalize(sunpos)), 0); sun_vec() = SGVec4f(toVec3f(normalize(sunpos)), 0);
@ -450,8 +424,8 @@ void FGLight::updateSunPos()
// calculate the sun's relative angle to local up // calculate the sun's relative angle to local up
SGVec3d viewPos = v->get_view_pos(); SGVec3d viewPos = v->get_view_pos();
SGQuatd hlOr = SGQuatd::fromLonLat(SGGeod::fromCart(viewPos)); SGQuatd hlOr = SGQuatd::fromLonLat(SGGeod::fromCart(viewPos));
SGVec3f world_up = toVec3f(hlOr.backTransform(-SGVec3d::e3())); SGVec3d world_up = hlOr.backTransform(-SGVec3d::e3());
SGVec3f nsun = toVec3f(normalize(sunpos)); SGVec3d nsun = normalize(sunpos);
// cout << "nup = " << nup[0] << "," << nup[1] << "," // cout << "nup = " << nup[0] << "," << nup[1] << ","
// << nup[2] << endl; // << nup[2] << endl;
// cout << "nsun = " << nsun[0] << "," << nsun[1] << "," // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
@ -461,62 +435,11 @@ void FGLight::updateSunPos()
SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = " SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
<< get_sun_angle() ); << get_sun_angle() );
// calculate vector to sun's position on the earth's surface // Get direction to the sun in the local frame.
SGVec3d rel_sunpos = sunpos - v->get_view_pos(); SGVec3d local_sun_vec = hlOr.transform(nsun);
// vector in cartesian coordinates from current position to the // Angle from south. XXX Is this correct in the southern hemisphere?
// postion on the earth's surface the sun is directly over double angle = atan2(local_sun_vec.x(), -local_sun_vec.y());
SGVec3f to_sun = toVec3f(rel_sunpos); set_sun_rotation(angle);
// printf( "Vector to sun = %.2f %.2f %.2f\n",
// v->to_sun[0], v->to_sun[1], v->to_sun[2]);
// Given a vector from the view position to the point on the
// earth's surface the sun is directly over, map into onto the
// local plane representing "horizontal".
// surface direction to go to head towards sun
SGVec3f surface_to_sun;
SGVec3f view_pos = toVec3f(v->get_view_pos());
surface_to_sun = map_vec_onto_cur_surface_plane(world_up, view_pos, to_sun);
surface_to_sun = normalize(surface_to_sun);
// cout << "(sg) Surface direction to sun is "
// << surface_to_sun[0] << ","
// << surface_to_sun[1] << ","
// << surface_to_sun[2] << endl;
// cout << "Should be close to zero = "
// << sgScalarProductVec3(nup, surface_to_sun) << endl;
// calculate the angle between surface_to_sun and
// v->get_surface_east(). We do this so we can sort out the
// acos() ambiguity. I wish I could think of a more efficient
// way. :-(
SGVec3f surface_east(toVec3f(hlOr.backTransform(SGVec3d::e2())));
float east_dot = dot( surface_to_sun, surface_east );
// cout << " East dot product = " << east_dot << endl;
// calculate the angle between v->surface_to_sun and
// v->surface_south. this is how much we have to rotate the sky
// for it to align with the sun
SGVec3f surface_south(toVec3f(hlOr.backTransform(-SGVec3d::e1())));
float dot_ = dot( surface_to_sun, surface_south );
// cout << " Dot product = " << dot << endl;
if (dot_ > 1.0) {
SG_LOG( SG_ASTRO, SG_INFO,
"Dot product = " << dot_ << " is greater than 1.0" );
dot_ = 1.0;
}
else if (dot_ < -1.0) {
SG_LOG( SG_ASTRO, SG_INFO,
"Dot product = " << dot_ << " is less than -1.0" );
dot_ = -1.0;
}
if ( east_dot >= 0 ) {
set_sun_rotation( acos(dot_) );
} else {
set_sun_rotation( -acos(dot_) );
}
// cout << " Sky needs to rotate = " << angle << " rads = " // cout << " Sky needs to rotate = " << angle << " rads = "
// << angle * SGD_RADIANS_TO_DEGREES << " degrees." << endl; // << angle * SGD_RADIANS_TO_DEGREES << " degrees." << endl;
} }

23
src/Time/sunsolver.cxx
View file

@ -58,16 +58,14 @@ void fgSunPositionGST(double gst, double *lon, double *lat) {
SGPropertyNode* sun = fgGetNode("/ephemeris/sun"); SGPropertyNode* sun = fgGetNode("/ephemeris/sun");
assert(sun); assert(sun);
double beta = sun->getDoubleValue("lat-deg");
// double r = globals->get_ephem()->get_sun()->getDistance();
double xs = sun->getDoubleValue("xs"); double xs = sun->getDoubleValue("xs");
double ys = sun->getDoubleValue("ys"); double ys = sun->getDoubleValue("ys");
double ye = sun->getDoubleValue("ye"); double ye = sun->getDoubleValue("ye");
double ze = sun->getDoubleValue("ze"); double ze = sun->getDoubleValue("ze");
alpha = atan2(ys - tan(beta)*ze/ys, xs); double ra = atan2(ye, xs);
delta = asin(sin(beta)*ye/ys + cos(beta)*ze); double dec = atan2(ze, sqrt(xs * xs + ye * ye));
tmp = alpha - (SGD_2PI/24)*gst; tmp = ra - (SGD_2PI/24)*gst;
if (tmp < -SGD_PI) { if (tmp < -SGD_PI) {
do tmp += SGD_2PI; do tmp += SGD_2PI;
while (tmp < -SGD_PI); while (tmp < -SGD_PI);
@ -77,7 +75,7 @@ void fgSunPositionGST(double gst, double *lon, double *lat) {
} }
*lon = tmp; *lon = tmp;
*lat = delta; *lat = dec;
} }
static double sun_angle( const SGTime &t, const SGVec3d& world_up, static double sun_angle( const SGTime &t, const SGVec3d& world_up,
@ -85,17 +83,18 @@ static double sun_angle( const SGTime &t, const SGVec3d& world_up,
SG_LOG( SG_EVENT, SG_DEBUG, " Updating Sun position" ); SG_LOG( SG_EVENT, SG_DEBUG, " Updating Sun position" );
SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t.getGst() ); SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t.getGst() );
double sun_lon, sun_gd_lat; double sun_lon, sun_gc_lat;
fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat ); fgSunPositionGST( t.getGst(), &sun_lon, &sun_gc_lat );
SGVec3d sunpos = SGVec3d::fromGeod(SGGeod::fromRad(sun_lon, sun_gd_lat)); SGVec3d sunpos = SGVec3d::fromGeoc(SGGeoc::fromRadM(sun_lon, sun_gc_lat,
SGGeodesy::EQURAD));
SG_LOG( SG_EVENT, SG_DEBUG, " t.cur_time = " << t.get_cur_time() ); SG_LOG( SG_EVENT, SG_DEBUG, " t.cur_time = " << t.get_cur_time() );
SG_LOG( SG_EVENT, SG_DEBUG, SG_LOG( SG_EVENT, SG_DEBUG,
" Sun Geodetic lat = " << sun_gd_lat ); " Sun Geocentric lat = " << sun_gc_lat );
// calculate the sun's relative angle to local up // calculate the sun's relative angle to local up
SGVec3f nup = normalize(toVec3f(world_up)); SGVec3d nup = normalize(world_up);
SGVec3f nsun = normalize(toVec3f(sunpos)); SGVec3d nsun = normalize(sunpos);
// cout << "nup = " << nup[0] << "," << nup[1] << "," // cout << "nup = " << nup[0] << "," << nup[1] << ","
// << nup[2] << endl; // << nup[2] << endl;
// cout << "nsun = " << nsun[0] << "," << nsun[1] << "," // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","