SDK Question How to calculate Sun elevation

[cristiapi]

I need to calculate the Sun elevation as seen from any point on Mars surface (the input is lat, lon on Mars and the output is the Sun elevation).

I'm trying to use oapiGlobalToLocal() and oapiGetPlanetObliquityMatrix(), but I don't get any useful result. Please, could somebody help me?

 

[boogabooga]

What about time?

 

[cristiapi]

What do you mean?

 

[boogabooga]

Your inputs are lat and lon.

That is insufficient to determine the sun elevation, which changes with time. (Dawn, noon, dusk, etc.)

 

[Keithth G]

Step 1: use OAPI's GetGlobalPosition function to return the position of the Sun, your position (i.e., the relevant point on the surface of Mars), and the centre of Mars.

Step 2: calculate a vector from the Sun to your position. Call it SunPos. Calculate a vector from the centre of Mars to your position. Call it MarsPos.

Step 3. set up a local co-ordinate system at your position to identify unit vectors that x_hat point East; y_hat points north and z_hat points vertically away from the centre of Mars. (Start with z_hat, since this depends simply on MarsPos and your position on the surface of the planet)

Step 4. project SunPos on to those coordinates by taking the 'dot product' of SunPos with each of these three vectors - x_hat, y_hat and z_hat. Call these three numbers (x', y', z').

Step 5. calculate ArcTan[z' / x']. This is the Sun's altitude.

N.B. This calculation assumes that the light-time from Sun to your location is zero - i.e., that light travels very fast; it assumes that Mars is spherical; and it assumes that there is no bending of light as it passes through the Martian atmosphere.

 

[cristiapi]

Step 2: calculate a vector from the Sun to your position. Call it SunPos. Calculate a vector from the centre of Mars to your position. Call it MarsPos.

How do I calculate those vectors? My skill in matrix calculation is poor.

Step 3. set up a local co-ordinate system at your position to identify unit vectors that x_hat point East; y_hat points north and z_hat points vertically away from the centre of Mars. (Start with z_hat, since this depends simply on MarsPos and your position on the surface of the planet)

Double dutch for me.
I posted in the SDK forum because I hope that there is some API functions to do most of the work (like the ones I posted).

Step 4. project SunPos on to those coordinates by taking the 'dot product' of SunPos with each of these three vectors - x_hat, y_hat and z_hat. Call these three numbers (x', y', z').

Step 5. calculate ArcTan[z' / x']. This is the Sun's altitude.

Ok for those steps.

N.B. This calculation assumes that the light-time from Sun to your location is zero - i.e., that light travels very fast; it assumes that Mars is spherical; and it assumes that there is no bending of light as it passes through the Martian atmosphere.

In this case I don't need very accurate result, but afaik, all the celestial bodies are spherical in Orbiter.

 

[Keithth G]

I have to do a few chores for a few hours. When I finish, I'll respond in more detail.

---------- Post added at 07:26 AM ---------- Previous post was at 12:17 AM ----------

OK, here's the same thing in a bit more detail - using Lua script-like commands:

sObj     =   Object Handle or the Sun
mObj     =   Object Handle for Mars
shObj    =   Object Handle for the location on Mars

sun      =   oapi.get_globalpos( sObj  )
mars     =   oapi.get_globalpos( mObj  )
ship     =   oapi.get_globalpos( shObj )

sunpos   =   vec.sub( sun,  ship)  -- a vector from the ship to the Sun
marspos  =   vec.sub( mars, ship)  -- a vector from the Ship to the centre of Mars


-- calculate a unit vector that points in the direction perpendicular to the surface
-- away from the centre of Mars.
z_hat    = - vec.div( marspos , sqrt( vec.dotp( marspos, marspos ) ) )


-- calculate the component of 'sunpos' that is perpendicular to the surface
alpha    =   vec.dotp( z_hat, sunpos )
sunpos_z =   vec.mul ( z_hat, alpha  )



-- calculate the component of 'sunpos' that lies in the plane at the surface
sunpos_p =   vec.sub ( subpos, subpos_z )
rsq        =   vec.dotp( sunpos_p, sunpos_p)


-- calculate the elevation of the Sun (in degrees)
angle    =   (180 /pi ) ArcTan( alpha / sqrt(rsq) )

 

[indy91]

Of course, if you have no ship but only the latitude and longitude, you can use oapiEquToGlobal from the API to get the global ship position for that time.

 

[cristiapi]

It works! :thumbup: Thanks a lot.

But since the function will be used in my Mars atmospheric module, I don't have any particular ship for the position, but only lat e lon.

Here is the C++ code for a point on Mars surface, for future use...

double lat= 12 * PI/180, lon= 123 * PI/180; VECTOR3 pos;
	
OBJHANDLE sObj= oapiGetGbodyByName("sun");
VECTOR3 sun; oapiGetGlobalPos(sObj, &sun);

OBJHANDLE mObj= oapiGetGbodyByName("mars");
VECTOR3 mars; oapiGetGlobalPos(mObj, &mars);

oapiEquToGlobal(mObj, lon, lat, oapiGetSize(mObj), &pos);
VECTOR3 sunpos= sun - pos; // Vector from pos to Sun
VECTOR3 marspos= mars - pos; // Vector from pos to the centre of Mars
// calculate a unit vector that points in the direction perpendicular to the surface away from the centre of Mars.
VECTOR3 z_hat= - marspos / length(marspos);

// Calculate the component of 'sunpos' that is perpendicular to the surface
double alpha= dotp(z_hat, sunpos);
VECTOR3 sunpos_z= alpha * z_hat;

// Calculate the component of 'sunpos' that lies in the plane at the surface
VECTOR3 sunpos_p= sunpos - sunpos_z;
double r= length(sunpos_p),
	EL= atan(alpha / r) * 180/PI;

 

[cristiapi]

Hi Keithth G,

Is there any chance to calculate Ls (the solar longitude asked in the other thread) using a similar method?
I currently use the "long method" that starts from the state of the body to calculate the true anomaly; very CPU expensive.

 

[Keithth G]

Probably

Too late for me to do much this evening. But I'll have a think about it tomorrow.

 

[Urwumpe]

Why not simply calculate the angle between sun (~Global position) and local radius vector (Relative position) and subtract it from 90°?

 

[cristiapi]

Are you referring to the solar zenith angle? Or just zenith angle in this link:
http://education.gsfc.nasa.gov/experimental/July61999siteupdate/inv99Project.Site/Pages/solar.insolation.html

 

[Keithth G]

Why not simply calculate the angle between sun (~Global position) and local radius vector (Relative position) and subtract it from 90°?

Since I think the goal is to calculate Martian seasons, at some point in the calculation one has to take into account the orientation of Mars' equatorial plane with respect to its orbital plane about the Sun. In concept, and just as a sketch, I think the way to do calculate the (Martian) solar longitude is:

For the date in question:

1: calculate the orientation of Mars' orbital plane about the Sun (i.e., find a vector that points in a direction that is perpendicular to its orbital plane);

2: calculate the orientation of Mars' equatorial plane (i.e. find a vector that points along the Martian axis of rotation);

3. from these two vectors find a vector that lies in Mars' equatorial plane that points to the (Martian) vernal (Spring) equinox - i.e., the 'place' where the Sun ascends though the Martian equatorial plane. This will be one of the two vectors that points along the 'line of nodes' of Mars' orbital plane and Mars' equatorial plane. The other vector will point to the 'other' equinox - i.e. where the Sun descends through the Martian equatorial plane.

4. define a co-ordinate system in which 'x' points towards the vernal equinox; 'y' is at right-angles to the x-axis but lies in the orbital plane; and 'z' is perpendicular to the Martian equatorial plane.

5. find the 'x' and 'y' co-ordinates of the Sun in this reference frame.

6. from these 'x' and 'y' co-ordinates use an 'atan2' function to calculate the angle of the Sun from the vernal equinox. This should be the Sun's longitude, 'Ls'.

I'll write up some psuedo-oapi code to do this.

(P.S., if it works, this idea can be extended to any other body).

---------- Post added at 07:14 AM ---------- Previous post was at 04:58 AM ----------

Try this:

oapi.set_pause(true) 

-- get the object handles for the Sun and for Earth
sun        = oapi.get_objhandle('Sun')
mars       = oapi.get_objhandle('Mars')

-- find the state vectors for the Sun
q_sun      = oapi.get_globalpos(sun)
p_sun      = oapi.get_globalvel(sun)

-- find the state vectors for Mars
q_mar      = oapi.get_globalpos(mars)
p_mar      = oapi.get_globalvel(mars) 

-- get two points on Mars' equatorial plane - separated by about 90 degrees
pos1	   = {lng = 0.0000, lat = 0.0000, rad = 3000000.0}
pos2       = {lng = 1.5708, lat = 0.0000, rad = 3000000.0}
marspos1   = oapi.equ_to_global(mars, pos1)
marspos2   = oapi.equ_to_global(mars, pos2)

-- calculate the unit vectors normal to the orbital and equatorial planes
vec1       = vec.sub( q_sun, q_mar )
vec2       = vec.sub( p_sun, p_mar )
vec1       = vec.div( vec1, vec.length( vec1 ) )
vec2       = vec.div( vec2, vec.length( vec2 ) )
normal1    = vec.crossp( vec2, vec1  )

vec1       = vec.sub( marspos1, q_mar)
vec2       = vec.sub( marspos2, q_mar)
vec1       = vec.div( vec1, vec.length( vec1 ) )
vec2       = vec.div( vec2, vec.length( vec2 ) )
normal2    = vec.crossp( vec2, vec1  )

-- calculate the unit vector that points towards the vernal equinox
vernal     = vec.crossp( normal1, normal2)
vernal     = vec.div( vernal, vec.length( vernal ) )

-- set up a co-ordinate system in Mars' orbital plane with the x-axis pointing to the
-- vernal equinox
x_hat      = vernal
z_hat      = normal1
y_hat      = vec.crossp( x_hat, z_hat)

-- calculate the earth-centric position of the Sun in the co-ordinate system
sunpos     = vec.sub(  q_sun,  q_mar )
x_coord    = vec.dotp( sunpos, x_hat )
y_coord    = vec.dotp( sunpos, y_hat )
z_coord    = vec.dotp( sunpos, z_hat )

-- calculate the solar longitude in degrees
longitude  = 180 * math.atan2( y_coord, x_coord ) / 3.141592654
if longitude < 0 then longitude = longitude + 360 end

oapi.set_pause(false)

I've tested this for a few spot dates using the widget provided here (http://www-mars.lmd.jussieu.fr/mars/time/martian_time.html)and the above algorithm appears to be in good agreement

 

[Urwumpe]

Why not use the rotation axis of Mars instead of some point in the equatorial plane... and then calculate normal2 which is the same as the rotation axis?

 

[Keithth G]

Why not use the rotation axis of Mars instead of some point in the equatorial plane... and then calculate normal2 which is the same as the rotation axis?

Sure. Do you know the oapi command for returning the rotation axis of a planet? I couldn't find it.

---------- Post added at 08:31 AM ---------- Previous post was at 08:03 AM ----------

P.S. Just in case others didn't quite catch the point that Urwumpe is making, there is a block of code in the above piece of lua script code:

pos1	   = {lng = 0.0000, lat = 0.0000, rad = 3000000.0}
pos2       = {lng = 1.5708, lat = 0.0000, rad = 3000000.0}
marspos1   = oapi.equ_to_global(mars, pos1)
marspos2   = oapi.equ_to_global(mars, pos2)

vec1       = vec.sub( marspos1, q_mar)
vec2       = vec.sub( marspos2, q_mar)
vec1       = vec.div( vec1, vec.length( vec1 ) )
vec2       = vec.div( vec2, vec.length( vec2 ) )
normal2    = vec.crossp( vec2, vec1  )

the sole purpose of which is to return Mars' rotation axis as a variable, 'normal2'. If one has some other means of calculating Mars rotation axis, then one can replace the above (correct but slightly long-winded) code with the more straightforward:

normal2 = 'the rotation axis of Mars as a vector in Orbiter's global coordinate system coordinates'

 

[indy91]

You can find the rotation axis of a planet by using equation 7 in the precession.pdf under Technotes. R(t) is the rotation matrix (oapiGetRotationMatrix).

 

[Keithth G]

Actually, I think this works quite well:

oapi.set_pause(true) 

-- get the object handles for the Sun and for Earth
sun        = oapi.get_objhandle('Sun')
mars       = oapi.get_objhandle('Mars')

-- find the state vectors for the Sun
q_sun      = oapi.get_globalpos(sun)
p_sun      = oapi.get_globalvel(sun)

-- find the state vectors for Mars
q_mar      = oapi.get_globalpos(mars)
p_mar      = oapi.get_globalvel(mars) 

-- calculate the unit vectors normal to mars' orbital planes
vec1       = vec.sub( q_sun, q_mar )
vec2       = vec.sub( p_sun, p_mar )
vec1       = vec.div( vec1, vec.length( vec1 ) )
vec2       = vec.div( vec2, vec.length( vec2 ) )
normal1    = vec.crossp( vec2, vec1  )

-- calculate mars' rotation axis
normal2    = oapi.equ_to_global(mars, {lng = 0.0, lat = 1.5707963267949, rad = 3000000.0}) 
normal2    = vec.sub( normal2, q_mar)
normal2    = vec.div( normal2, vec.length( normal2 ) )

-- calculate the unit vector that points towards the vernal equinox
vernal     = vec.crossp( normal1, normal2)
vernal     = vec.div( vernal, vec.length( vernal ) )

-- set up a co-ordinate system in Mars' orbital plane with the x-axis pointing to the
-- vernal equinox
x_hat      = vernal
z_hat      = normal1
y_hat      = vec.crossp( x_hat, z_hat)

-- calculate the earth-centric position of the Sun in the co-ordinate system
sunpos     = vec.sub(  q_sun,  q_mar )
x_coord    = vec.dotp( sunpos, x_hat )
y_coord    = vec.dotp( sunpos, y_hat )

-- calculate the solar longitude in degrees
longitude  = 180 * math.atan2( y_coord, x_coord ) / 3.141592654
if longitude < 0 then longitude = longitude + 360 end

oapi.set_pause(false)

 

[cristiapi]

I compared the result of your code with the value returned by SPICE (with no correction for the aberration). In order to obtain from Orbiter the same body state as in SPICE, I used CeBoMo.
I calculated only few values for Mars, Titan and Pluto and the differences are less than 0.4 deg.
For my needs, your code is perfectly good. :thumbup: I'll use it in my next version of Mars atmosphere...
Thanks a lot #2.

I changed the name of the variables related to Mars with a generic 'bdy'.
Probably the function oapiGetPlanetCurrentRotation(OBJHANDLE hPlanet) can be used to calculate the body rotation axis, but I don't know how.

OBJHANDLE sun= oapiGetGbodyByName("sun");
OBJHANDLE bdy= oapiGetGbodyByName("mars");

// find the state vectors for the Sun
VECTOR3 q_sun; oapiGetGlobalPos(sun, &q_sun);
VECTOR3 p_sun; oapiGetGlobalVel(sun, &p_sun);

// find the state vectors for bdy
VECTOR3 q_bdy; oapiGetGlobalPos(bdy, &q_bdy);
VECTOR3 p_bdy; oapiGetGlobalVel(bdy, &p_bdy);

// calculate the unit vector normal to bdy orbital planes
VECTOR3 vec1= q_sun - q_bdy; vec1 /= length(vec1);
VECTOR3 vec2= p_sun - p_bdy; vec2 /= length(vec2);
VECTOR3 normal1= crossp(vec2, vec1);

// calculate bdy rotation axis
VECTOR3 normal2; oapiEquToGlobal(bdy, 0.0, PI/2, 3000000.0, &normal2);
normal2= (normal2 - q_bdy) / length(normal2);

// calculate the unit vector that points towards the vernal equinox
VECTOR3 vernal= crossp(normal1, normal2);
vernal /= length(vernal);

// set up a co-ordinate system in bdy orbital plane with the x-axis pointing to the vernal equinox
VECTOR3 x_hat= vernal;
VECTOR3 z_hat= normal1;
VECTOR3 y_hat= crossp(x_hat, z_hat);

// calculate the bodycentric position of the Sun in the co-ordinate system
VECTOR3 sunpos= q_sun - q_bdy;
double x_coord= dotp(sunpos, x_hat);
double y_coord= dotp(sunpos, y_hat);

// calculate the solar longitude in degrees
double Ls= atan2(y_coord, x_coord) * (180/PI);
if(Ls < 0) Ls += 360;

---------- Post added at 12:43 ---------- Previous post was at 12:17 ----------

Probably the differences with SPICE come from a different reference frame.
Here is what SPICE calculates:

Ls is the longitude of the body-sun vector in a right-handed frame whose basis vectors are defined as follows:

- The positive Z direction is given by the instantaneous angular velocity vector of the orbit of the body about the sun.

- The positive X direction is that of the cross product of the instantaneous north spin axis of the body with the positive Z direction.

- The positive Y direction is Z x X.

 

[Keithth G]

Looking carefully at those definitions, I'm pretty sure that my x-y-z co-ordinate system matches those definitions - at least, in the principle of what is going on. So it is a little disturbing that there is still a sizeable difference.

However, there may be some subtle differences in interpretation in calculating instantaneous angular momentum, e.g.:

1. When calculating the instantaneous angular momentum vector, should one subtract off the Sun's instantaneous position/velocity; or should one reference the Solar System barycentre - i.e., the global coordinate system's origin (0,0,0)?

2. Does SPICE use mean orbital elements to calculate the Sun's position - or does it use a full ephemeris solution?

2. Does Orbiter use the same spin axes for bodies as SPICE?

Anyway, as you say, for your calculations it doesn't make much difference.