|
|
"tomcee" <tomcees_pc@xxxxxxxxx> wrote in message
news:24b3d18e-ac7f-46e9-a0a8-3ba9057ff035@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
> Given a table of sunset times, such as:
<snip>
> Has anyone determined the basic functions contained within this
> 'Sunset function'? Given the basic functions, I can then calculate the
> constants.
> Thanks in advance for your help,
> TomCee
My reply is going to infuriate everyone because of its horrendous
complexity, but I thought I would brush up on my vector algebra and solve
this problem by...
(1) considering the earth's orbit as an ellipse
(2) using Kepler's "equal area for equal time" rule of orbit speed
(3) general mapping of 3D coordinates from one basis (Earth's lat/lon grid)
to another (Earth's orbital plane 23.5 degrees oblique)
(4) the astronomer's publishing of GMT for perihelion and summer solstice
(5) looking iteratively for when the vector from the sun to the observer is
basically perpendicular to the vector from the earth's center to the
observer
Well, after working almost a month I've got a 'C' program that seems to get
within 6 minutes of the U.S. Navy's published times. I actually learned a
lot about vector operations. For now I'm going to put this thing to rest,
but here's the code...
=============================================================================================================
/*
SUNSET.C a program to calculate the time of sunset for any given latitude
and longitude, year, and day
of the year. It uses vector algebra as much as possible.
Sunset occurs when a line from the sun (considered a point source) runs
tangent to the earth and touches
earth exactly at the point of the observer's latitude and longitude.
Consider this to basically be a problem
in 3D coordinate space using rectangular coordinates (although polar
coordinates are used initially to compute
the earth's orbit).
The 3D coordinate space has the sun at the origin. The earth's orbit is an
ellipse in the X-Y plane with the
sun at one of the foci. When the earth is at perihelion, or closest to the
sun, it is on the positive X-axis
with Y=0.
We need two equations. One will determine the earth's center as a vector
from the sun as a function of time.
Call it the 'sun-earth' vector. The second will determine the distance and
direction from the center of the
earth to the observer's position on the surface as a function of time.
Call it the 'earth-observer' vector.
The sum of these two vectors will be a vector that goes from the sun to
the observer's position on the surface
of the earth. This third, computed, vector is called the 'sun-observer'
vector. When the 'sun-observer' vector
is perpendicular to the 'earth-observer' vector then the sun is sitting
on the horizon and it's either sunrise
or sunset. Perpendicular is tested for by taking the dot-product of the
two vectors 'earth-observer' and
'sun-observer' and testing for zero.
CALCULATING THE sun-earth VECTOR AS A FUNCTION OF TIME t
Using polar coordinates with the sun at the origin and zero degrees being
perihelion Kepler devised the
following two equations...
(1) M = 2 * PI * t / T
where 'M' is the "mean" anomaly, 't' is the time since perihelion, and
'T' is the orbital period
(2) M = E - e * sin (E)
where 'E' is the "eccentric" anomaly, and 'e' is the eccentricity of
the ellipse
The goal here is to first calculate 'M' from 't'. 'T' is 365 days 5 hours
49 minutes, or 31556940 seconds. Then
we use the second equation to solve for 'E'. Unfortunately, equation (2)
cannot be solved via algebra. But a
formula involving an infinite series will get it accurately enough with
about three terms...
E = M + (e - 1/8 * e^3) * sin (M) + 1/2 * e^2 * sin (2 * M) + 3/8 * e^3
* sin(3 * M) + ...
Once we get 'E' we substitute this value into the following equation to get
'v' (angle of the earth from the
sun)...
(3) v = arctan ( sqrt ((1 + e) / (1 - e)) * tan ( E / 2) ) * 2
where 'v' is the polar coordinate angle (from the sun) for the earth at
time 't'
Then we simply plug 'v' into the equation for an ellipse (in polar
coordinates) to get 'r' (the distance from
the sun)...
(4) r = a * (1 - e^2) / (1 + e * cos (v) )
where 'a' is the semi-major axis
A final conversion of 'v' and 'r' into cartesian coordinates is made by
x = r * cos ( v )
y = r * sin ( v )
VISUALIZING THE EARTH'S SPHERICAL SHAPE, AXIS OF ROTATION, AND
SEASON-CREATING TILT OF THE AXIS
Next we have to compute the vector 'earth_observer' based on any time
't'. This vector represents the direction
and distance of travel from the center of the earth to the observer's
location on the surface. This
'earth_observer' vector, expressed in the basis unit vectors of our
earth-orbit coordinate system, is computed
based on three new factors...
(1) Latitude and Longitude
The observer is on the surface of a sphere with an approximately 4000
mile radius. Using the observer's
latitude and longitude will allow us to compute the surface
coordinates and, thus, a basic 'observer'
vector from the earth's center to the observer. To simplify this
calculation we will invoke a new,
secondary, coordinate system aligned with the latitude and longitude
grid. The X-axis of this new system
starts from the earth's center and leaves the earth where the prime
meridian intersects with the equator.
The Z-axis of this new system is the axis of rotation of the earth and
starts at the earth's center and
leaves the earth through the north pole. The Y-axis of this new system is
now determined by the other two
axes. It starts at the earth's center and leaves the earth where the 90
degree east longitude meridian
intersects the equator. Using this secondary coordinate system we compute
the observer's location with
the three equations...
Z= 4000*sin(lat)
Y=(4000*cos(lat))*sin(lon)
X=(4000*cos(lat))*cos(lon)
Notice that the familiar notation of east and west longitude designation
has been replaced with the
mathematical requirement that east longitude is used as-is (a positive
number from 0 to 180) but west
longitude is negated in sign (so 0 to 180 west longitude becomes 0
to -180 degrees).
(2) Elapsed Time
The earth rotates on its own axis approximately every 23 hours 56
minutes 4 seconds (a sidereal day). So
even though the observer is fixed in geographic latitude and longitude
his effective longitude increases
by 360 degrees every 23h56m4s. Therefore, an incremented longitude
computed from the time 't' is added to
the basic longitude specified in factor (1) to give a correct computation
of the 'earth_observer' vector
including both factors (1) and (2). At this point the 'earth_observer'
vector is still relative to the
secondary coordinate system. The equations now become...
Z= 4000*sin(lat)
Y=(4000*cos(lat))*sin(lon+t/SIDEREAL_DAY*360)
X=(4000*cos(lat))*cos(lon+t/SIDEREAL_DAY*360)
(3) Coordinate Conversion from Secondary to Primary Coordinate System for
Earth's Axial Tilt
Because the earth's axis of rotation is "tilted" at an angle of 23.5
degrees with respect to the orbital
plane the basic calculation of the 'earth_observer' vector done by
using factors (1) and (2) with respect
to a secondary coordinate system based on latitude and longitude will
have to be converted to a vector based
on the original coordinate system that was used to describe the earth's
elliptical orbit before we can
use it in an algebraic solution for the time of sunset.
Returning to the original coordinate system the earth's north pole is
above the X-Y plane, always in
positive Z space and "roughly" near the top of the sphere that is the
earth. Due to the earth's tilt a
line from the earth's center through the north pole will always make an
angle of 23.45 degrees with the
Z-axis of the original coordinate system. The direction of this line in
the original coordinate system is
called the 'tilt' vector. It is a constant vector and does not change
over time. It never changes its
direction. Its specific X, Y, and Z components can only be determined by
researching astronomers' analysis
of earth's orbit. This 'tilt' vector the astronomers give us is the
Z-axis of our secondary coordinate
system. What are the X and Y axes? Astronomers will also have to give us
those since they determine the
orientation of the prime meridian to the original coordinate system at
time t=0. Once we have all three
of these secondary coordinate system axes defined *as components of the
original coordinate system* then we
use the transformation that follows for any coordinate system
conversion. Thus, the astronomers give us
three vectors aligned with the axes of our secondary coordinate system
but using the original coordinate
system's basis vectors. We call them 'tilt_x', 'tilt_y', and 'tilt_z'.
'tilt_z' is calulated by examining the earth's axis of rotation at the
time of summer solstice. 'tilt_x' is
also computed by noting where the prime meridian is positioned at the
time of summer solstice and then
"unwinding" it around the earth's axis of rotation back to the time t=0.
'tilt_y' is simply the cross
product of 'tilt_z' and 'tilt_x'. The rotation of the initial 'tilt_x'
vector constructed from the geometry
at summer solstice back to time t=0 is done using Rodrigues' solution for
rotating vector X around vector U
by angle theta...
_ _ _ _ _ _ _
X' = X * cos(theta) + U * (U . X)(1 - cos(theta)) + (U x X) *
sin(theta) //dot-prod, cross-prod, respectively
CONVERSION OF COORDINATES BETWEEN TWO COORDINATE SYSTEMS
The coordinate conversion from the secondary system to the original
system is accomplished via the
general-purpose linear transformation algebraic solution found on the top
of page 534 of my VNR book.
This solution handles the general case of both a translation of the
origin and a rotation of the axes for
two distinct coordinate systems. It requires the calculation of
"direction cosines" for each of nine angles.
Given any vectors aligned with the three axes for the two coordinate
systems these nine direction cosines
can be computed by using the dot-product relationship of two vectors.
Because the dot-product of two
vectors is the product of the magnitudes of both vectors multiplied by
the cosine of the angle between them
the cosine itself is the dot-product divided by the product of the
magnitudes of both vectors. The
calculation of the nine direction cosines follows:
a11 = dot-product(unit_x,tilt_x) / ( |unit_x|*|tilt_x| )
a12 = dot-product(unit_x,tilt_y) / ( |unit_x|*|tilt_y| )
a13 = dot-product(unit_x,tilt_z) / ( |unit_x|*|tilt_z| )
a21 = dot-product(unit_y,tilt_x) / ( |unit_y|*|tilt_x| )
a22 = dot-product(unit_y,tilt_y) / ( |unit_y|*|tilt_y| )
a23 = dot-product(unit_y,tilt_z) / ( |unit_y|*|tilt_z| )
a31 = dot-product(unit_z,tilt_x) / ( |unit_z|*|tilt_x| )
a32 = dot-product(unit_z,tilt_y) / ( |unit_z|*|tilt_y| )
a33 = dot-product(unit_z,tilt_z) / ( |unit_z|*|tilt_z| )
The direction cosines are the 'Amn' coefficients in the following
equations. The translational distance
between the origins of the two coordinate systems are the 'An'
constants.
The equations for converting a coordinate in the 'prime' coordinate
system to a coordinate in the 'non-prime'
coordinate system are on the left-hand side. The equations for
converting a coordinate in the 'non-prime'
coordinate system to a coordinate in the 'prime' coordinate system are
on the right-hand side.
x = a1 + a11 * x' + a12 * y' + a13 * z' | x' = a11 * (x - a1) +
a21 * (y - a2) + a31 * (z - a3)
y = a2 + a21 * x' + a22 * y' + a23 * z' | y' = a12 * (x - a1) +
a22 * (y - a3) + a32 * (z - a3)
z = a3 + a31 * x' + a32 * y' + a33 * z' | z' = a13 * (x - a1) +
a23 * (y - a3) + a33 * (z - a3)
In both sets of equations a1, a2, and a3 are the coordinates (x, y, and
z) of the origin of the 'prime'
system with respect to the 'non-prime' system. That is, the origin point
{0', 0', 0'} of the 'prime' system
is at location {a1, a2, a3} in the 'non-prime' system.
In both sets of equations the 'Amn' direction cosines refer to the
cosine of the angle between axis 'm' of
the 'non-prime' system and axis 'n' of the 'prime' system (where m=1 is
the X-axis of the 'non-prime' system,
m=2 is the Y-axis of the 'non-prime' system, and m=3 is the Z-axis of
the 'non-prime' system; where n=1 is
the X-axis of the 'prime' system, n=2 is the Y-axis of the 'prime'
system, and n=3 is the Z-axis of the
'prime' system). Notice that because of the symmetry of the cosine
function about the zero angle it makes no
difference whether the angle is rotated through according to a
right-handed orientation or a left-handed
orientation since cos(theta)=cos(360-theta). However, each of the nine
angles rotated through must all be
determined by choosing representative vectors of the axes of each
coordinate system that have the same
orientation to those axes with regard to positive (or negative)
direction.
Since we are dealing with vectors and not actual coordinates there is no
translation of the origin of our
two systems to be considered. a1, a2, and a3 are all zero. We will have
coordinates in the secondary system
that need to be converted to the primary system and the equations on the
left will be used.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "VCTRSUPPORT.h"
//SOME USEFUL CONSTANTS
#define PI 3.141592653589793
#define SIDEREAL_DAY (23*3600.0+56*60.0+4) //Number of seconds in
a sidereal day
#define SUN_RISE_SET_ANGLE 89.167 //Angle between
'sun_earth' & 'earth_observer' at sunrise/sunset
#define SEMI_MAJOR_AXIS 92955807.0 //Semi-major axis of
earth's orbit ellipse
double secs_per_orbit=365*86400+6*3600+13*60+53; //Number of seconds in
a true orbit (anomalistic) of the sun
double e=0.0167; //Earth's elliptical
orbit eccentricity
#define ECCENTRICITY 0.0167 //Earth's elliptical
orbit eccentricity
double a=92955807.0; //Semi-major axis of
earth's orbit ellipse
double obs_lat=40.267; //Observer's latitude (positive is northern
hemispere, negative is southern hemisphere)
double obs_lon=-80.117; //Observer's longitude (positive is east of
Greenwich, negative is west of Greenwich)
//double obs_lat=0.0; //Observer's latitude (positive is northern
hemispere, negative is southern hemisphere)
//double obs_lon=0.0; //Observer's longitude (positive is east of
Greenwich, negative is west of Greenwich)
typedef struct {
int yr; //2008, 2009, 2010, etc
int mo; //1=JAN, 2=FEB, etc
int da; //1-31
int hh; //0-23
int mm; //0-59
int ss; //0-59
int tzone; //0=GMT -5=EST, etc
} datim;
// YYYY MM DD HH MM SS Z
datim perihelion={2008, 1, 3, 0, 0, 0,0}; //Date/time (GMT) of perihelion
(earth closest to sun, t=0)
datim sp_equinox={2008, 3,20, 5,48, 0,0}; //Date/time (GMT) of spring
equinox (sun_earth perpend to tilt_z)
datim s_solstice={2008, 6,20,23,59, 0,0}; //Date/time (GMT) of summer
solstice (tilt_z leans toward sun)
datim aphelion ={2008, 7, 4, 8, 0, 0,0}; //Date/time (GMT) of aphelion
(earth farthest from sun)
datim fa_equinox={2008, 9,22,15,44, 0,0}; //Date/time (GMT) of fall
equinox (sun_earth perpend to tilt_z again)
datim w_solstice={2008,12,21,12, 4, 0,0}; //Date/time (GMT) of winter
solstice (tilt_z leans away from sun)
datim current_date;
POINT sun={0,0,0}; //Fixed position of the sun def'd in *primary
coords*. It never moves.
//The three vectors for the secondary coord system (described in primary
coord system components)
VECTOR tilt_z; //Z-axis of secondary coord syst def'd in *primary
coords* (earth's rotation axis)
VECTOR tilt_x; //X-axis of secondary coord syst def'd in *primary
coords* (prime meridian at t=0)
VECTOR tilt_y; //Y-axis of secondary coord syst def'd in *primary
coords* (90 degrees east longitude)
//Three unit vectors described in the primary coord system used for
computing direction cosines
VECTOR unit_x={1,0,0};
VECTOR unit_y={0,1,0};
VECTOR unit_z={0,0,1};
//The nine direction cosines
double a11,a12,a13,a21,a22,a23,a31,a32,a33;
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - -
//Function prototypes
int time_to_date(int t,datim *d); //Not coded yet. Not needed.
int date_to_time(datim *d,long *t);
int adjust_date (datim *dt,long time);
int calc_angle (datim *dt,double *angle);
int print_vec (VECTOR *v);
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - -
int main() {
int it; //Time 't' as an integer
double t; //Time in seconds from initial position of earth
double M; //The "mean" anomaly angle
double E; //The "eccentric" anomaly angle
double v; //Angle, in radians, of earth in its orbit from
its initial position
double r; //Polar coordinate 'radius' of earth
VECTOR sun_earth; //Vector from sun to earth's center
VECTOR normal; //Computed vector normal to orbital plane
through the sun
VECTOR earth_observer; //Vector from earth's center to the observer at
given latitude, longitude, and time
POINT obs; //Coordinates of observer's point (after time
rotations) in *secondary coordinates*
VECTOR obs_vec; //Vector from earth's center to observer's point
on surface in *secondary coordinates*
double lat,lon;
VECTOR sun_observer; //Vector from sun to the observer's tangent
point
double dp,angle;
double dangle;
long delta_t;
double temp;
VECTOR tempu,tempx,tempv1,tempv2;
double theta;
/*
OK, let's get started! First, calculate the constant 'tilt_z' vector (the
earth's axis of rotation) in the
primary coordinate system based on the fact that at summer solstice the
axis of the earth (which always forms
an angle of 23.45 degrees with the normal to the orbital plane) is, at that
moment, pointing *directly* at the
sun. Another way to put it is that the normal vector to the orbital
plane, passing through the sun, will
intersect the earth's extended axis at some point in the positive Z
direction. Thus, a right triangle is formed
by the earth, the sun, and that point of intersection in positive Z
space. (This condition also happens at
winter solstice, but then the intersection is in negative Z space.) The
'sun_earth' vector at time
t=summer_solstice is the short leg of this right triangle. The normal
vector to the earth's orbital plane and
passing through the sun is the long leg of this triangle. The hypotenuse
is the 'tilt_z' vector for the earth.
The tangent of 23.45 degrees (0.4337751161) is the ratio of the
'sun_earth' vector's length to the length of the
normal vector. So the length of the normal vector is computed by dividing
the length of the 'sun_earth' vector
by the tangent of 23.45 degrees.
*/
if (date_to_time(&s_solstice,&it)) { //Convert summer solstice
date/time to seconds since perihelion
printf("s_solstice conversion error to basic time\n");
goto bad_exit;
}
t=it; //Make it floating point
// printf("time (t) of summer solstice is %.1f\n",t);
//Now solve for the X and Y coordinates of the earth at summer solstice to
get 'sun_earth' vector
M=2*PI*t/secs_per_orbit; //Compute Kepler's 'Mean' anomaly
E=M //Compute Kepler's 'Eccentric'
anomaly (approximately)
+(e-1.0/8.0*pow(e,3.0))*sin(1.0*M) //using the first three terms of an
infinite series
+ 1.0/2.0*pow(e,2.0) *sin(2.0*M) // "
+ 3.0/8.0*pow(e,3.0) *sin(3.0*M); // "
v=atan(sqrt((1.0+e)/(1.0-e))*tan(E/2.0))*2.0; //Use 'Eccentric' anomaly
to solve for polar coord angle from sun
r=a*(1.0-pow(e,2.0))/(1.0+e*cos(v)); //Equation for ellipse
gives radius from polar coord angle
sun_earth.x=r*cos(v); //Convert polar coord angle
and radius to cartesian coords
sun_earth.y=r*sin(v); // "
sun_earth.z=0.0; // "
// printf("sun_earth="); print_vec(&sun_earth); printf("\n");
//OK, now have 'sun_earth' vector. Solve for normal vector's length from
right triangle relationship.
normal.x=0.0; //Normal vector has no X
component
normal.y=0.0; //Normal vector has no Y
component
normal.z=vector_length(&sun_earth)/tan(23.45*PI/180); //Z length is
'sun_earth' length divided by tan(23.45)
// printf("normal(+Z)="); print_vec(&normal); printf("\n");
//The 'tilt_z' vector is the sum of the 'sun_earth' vector (negated) and the
'normal' vector
tilt_z.x=-sun_earth.x+normal.x; //Add 'em
tilt_z.y=-sun_earth.y+normal.y; // "
tilt_z.z=-sun_earth.z+normal.z; // "
// printf("tilt_z="); print_vec(&tilt_z); printf("\n");
/*
OK, now have 'tilt_z' vector described in primary coordinates. We now
need the 'tilt_x' and 'tilt_y' vectors.
These are both perpendicular to the 'tilt_z' vector and the 3 of them are
the axes for the secondary coordinate
system. 'tilt_x' points from the earth's center through the earth's
equator at the longitude of the prime
meridian at time t=0. It will also be constant, not changing as the earth
orbits or rotates. 'tilt_x' will tell
us the exact orientation of the earth in its rotation about its axis at
t=0. 'tilt_y' is nothing more than the
vector from the earth's center through the equator at longitude east 90
degrees. It will be calculated via a
cross product of 'tilt_z' and 'tilt_x' once 'tilt_x' has been
ascertained.
Working with the summer solstice again, this time let's make a right
triangle using the 'sun_earth' vector as
the long leg. The normal vector to the earth's orbital plane going
through the sun and going in the *negative*
Z direction is the short leg. The hypotenuse is the 'tilt_x' vector at
12:00AM at the beginning of the date of
summer solstice. If we then back out all the rotations of this 'tilt_x'
vector that have happened since the
date/time of perihelion then we will have the correct 'tilt_x' vector value
we want. Now the tangent of 23.45
degrees is the ratio of the length of this negative-going normal vector to
the 'sun_earth' vector's length. So
the length of this negative-going normal vector is the product of the
'sun_earth' vector's length and the
tangent of 23.45 degrees.
*/
normal.x=0.0;
normal.y=0.0;
normal.z=-vector_length(&sun_earth)*tan(23.45*PI/180);
// printf("normal(-Z)="); print_vec(&normal); printf("\n");
//The 'tilt_x' vector (as a starting point) is the sum of the 'sun_earth'
vector and the 'normal' vector (negated)
tilt_x.x=sun_earth.x-normal.x;
tilt_x.y=sun_earth.y-normal.y;
tilt_x.z=sun_earth.z-normal.z;
// printf("tilt_x(12:00AM sum_solstice)="); print_vec(&tilt_x);
printf("\n");
//Unwind tilt_x vector to its position at perihelion (find the angle of
clockwise rotation needed)
current_date.yr=s_solstice.yr;
current_date.mo=s_solstice.mo;
current_date.da=s_solstice.da;
current_date.hh=0;
current_date.mm=0;
current_date.ss=0;
current_date.tzone=0;
if (date_to_time(¤t_date,&it)) { //Convert 12:00AM on summer
solstice to 't' secs since perihelion
printf("error converting date to time\n");
goto bad_exit;
}
t=it; //Convert time to float
// printf("time (t) of *midnight* on day of summer solstice is %.1f\n",t);
theta=-t/SIDEREAL_DAY*2*PI; //Get angle to unwind
// printf("angle to unwind=%.1f\n",theta);
//Use rotation of vector around a vector to get the unwound 'tilt_x' vector
at t=0
/*
The rotation of vector X around *unit* vector U by an angle 'theta' is
given by...
_ _ _ _ _ _ _
X' = X * cos(theta) + U * (U . X)(1 - cos(theta)) + (U x X) * sin(theta)
//dot-prod & cross-prod, respectively
*/
tempu=tilt_z; //Vector U is our tilt_z (earth's
axis of rotation). Make a copy of it.
vector_resize(&tempu,1.0); //Make it a unit vector for the
following calculations
tempx=tilt_x; //Make copy of X to begin with
vector_mpy(&tempx,cos(theta)); //Now have X*cos(theta)
temp=dot_product(&tempu,&tilt_x)*(1.0-cos(theta)); //Now have
(U.X)(1-cos(theta))
tempv1=tempu; //Prepare to multiply U by result in
previous line
vector_mpy(&tempv1,temp); //Now have U*(U.X)(1-cos(theta))
cross_product(&tempu,&tilt_x,&tempv2); //Now have (UxX)
vector_mpy(&tempv2,sin(theta)); //Now have (UxX)*sin(theta)
vector_add(&tempv1,&tempv2,&tilt_x); //Now have
U*(U.X)(1-cos(theta))+(UxX)*sin(theta)
vector_add(&tilt_x,&tempx,&tilt_x); //Now have
X*cos(theta)+U*(U.X)(1-cos(theta))+(UxX)*sin(theta)
// printf("tilt_x(final, at t=0)="); print_vec(&tilt_x); printf("\n");
//Take cross product of 'tilt_z' with 'tilt_x' to generate our 'tilt_y'
perpendicular to both
cross_product(&tilt_z,&tilt_x,&tilt_y);
//Compute the nine direction cosines for secondary to primary axes
conversion
a11=dot_product(&unit_x,&tilt_x)/(vector_length(&unit_x)*vector_length(&tilt_x));
a12=dot_product(&unit_x,&tilt_y)/(vector_length(&unit_x)*vector_length(&tilt_y));
a13=dot_product(&unit_x,&tilt_z)/(vector_length(&unit_x)*vector_length(&tilt_z));
a21=dot_product(&unit_y,&tilt_x)/(vector_length(&unit_y)*vector_length(&tilt_x));
a22=dot_product(&unit_y,&tilt_y)/(vector_length(&unit_y)*vector_length(&tilt_y));
a23=dot_product(&unit_y,&tilt_z)/(vector_length(&unit_y)*vector_length(&tilt_z));
a31=dot_product(&unit_z,&tilt_x)/(vector_length(&unit_z)*vector_length(&tilt_x));
a32=dot_product(&unit_z,&tilt_y)/(vector_length(&unit_z)*vector_length(&tilt_y));
a33=dot_product(&unit_z,&tilt_z)/(vector_length(&unit_z)*vector_length(&tilt_z));/*
OK, now test all our logic by choosing some date/time values and
thencalculating the 'sun_earth' vector, the 'earth_observer' vector, the
'sun_observer' vector, and the angle betweenthe 'sun_observer' and
'earth_observer' vectors. When the angle is close to 90 degrees then wehave a
sunrise or sunset.*//* current_date.yr=2008; current_date.mo=4;
current_date.da=3; current_date.hh=0; current_date.mm=0;
current_date.ss=0; current_date.tzone=-5; for
(;current_date.hh<24;current_date.hh++) { if (date_to_time(¤t_date,&it))
{ //Convert current_date to a time't' in seconds since perihelion
printf("error converting date to time\n"); goto bad_exit; } t=it;
//Convert time to float//First calculate sun_earth vector
M=2*PI*t/secs_per_orbit; //Mean anomaly
E=M+(e-1.0/8.0*pow(e,3.0))*sin(M)+1.0/2.0*pow(e,2.0)*sin(2.0*M)+3.0/8.0*pow(e,3.0)*sin(3.0*M);
//Approx v=atan(sqrt((1.0+e)/(1.0-e))*tan(E/2.0))*2.0; //Solved for polar
anglefrom sun r=a*(1.0-pow(e,2.0))/(1.0+e*cos(v)); //Equation for ellipse
gives radiusas funtion of polar angle sun_earth.x=r*cos(v); //Convert
polar to cartesian sun_earth.y=r*sin(v); // " sun_earth.z=0.0; // "
lon=obs_lon+t/SIDEREAL_DAY*360.0; //Compute longitude with timerotation done
lat=obs_lat ; //Latitude always stays the same obs.z=
4000.0*sin(lat*PI/180.0); //Observer'scoords in secondary
system at time 't' obs.y=(4000.0*cos(lat*PI/180.0))*sin(lon*PI/180.0);
//Observer'scoords in secondary system at time 't'
obs.x=(4000.0*cos(lat*PI/180.0))*cos(lon*PI/180.0); //Observer'scoords in
secondary system at time 't'//Change coords in secondary system to vector in
secondary system (trivial) obs_vec.x=obs.x; obs_vec.y=obs.y;
obs_vec.z=obs.z; earth_observer.x = a11 * obs_vec.x + a12 * obs_vec.y + a13
* obs_vec.z;//Now transform to primary system earth_observer.y = a21 *
obs_vec.x + a22 * obs_vec.y + a23 * obs_vec.z;// " earth_observer.z = a31 *
obs_vec.x + a32 * obs_vec.y + a33 * obs_vec.z;// "//Vector from sun to observer
is sum of vector from sun to earth plusobserver's vector from earth center
sun_observer.x=sun_earth.x+earth_observer.x;
sun_observer.y=sun_earth.y+earth_observer.y;
sun_observer.z=sun_earth.z+earth_observer.z;//The dot-product of the 'observer'
vector and 'sun_tangent' vector must bezero since they are at right angles
dp=dot_product(&sun_observer,&earth_observer);
angle=acos(dp/(vector_length(&sun_observer)*vector_length(&earth_observer)));
printf("%8.1f %4d/%02d/%02d %02dh %02dm
%02ds%5.1f\n",t,current_date.yr,current_date.mo,current_date.da,
current_date.hh,current_date.mm,current_date.ss,
angle*180.0/PI); }*//* Given a date/time determine the time of sunrise and
sunset for that date.Do this by iterating times during the date and check that
the angle of the sun is close enough toSUN_RISE_SET_ANGLE.*/ for
(current_date.da=1;current_date.da<=30;current_date.da++) {
current_date.yr=2008; current_date.mo=4;// current_date.da=3;
current_date.hh=12; //First use noon - sun isdefinitely
up! current_date.mm=0; current_date.ss=0; current_date.tzone=-5;//First, work
earlier in the day looking for snrise if (calc_angle(¤t_date,&dangle)) {
//Get the angle for noon printf("bad return from 'calc_angle'\n"); goto
bad_exit; } if (dangle<SUN_RISE_SET_ANGLE) { //If sun is below the
horizonat noon (!) printf("error: sun below horizon at noon!\n"); goto
bad_exit; } for (delta_t=-3600;abs(delta_t)>0;) { //Iterate to earlier
timestarting with one hour movement adjust_date(¤t_date,delta_t);
//Use 'delta_t' to getanother date/time to look at if
(calc_angle(¤t_date,&dangle)) { //Get the angle for this new time
printf("bad return from 'calc_angle'\n"); goto bad_exit; } if
(dangle<SUN_RISE_SET_ANGLE) { //If we're earlier than sunrise if
(delta_t<0) //If we were moving toward earliertimes
delta_t/=-2; //Halve the time increment andreverse its
direction continue; //Stay in the iteration }
if (dangle>SUN_RISE_SET_ANGLE) { //If we're later than sunrise if
(delta_t>0) //If we were moving toward latertimes
delta_t/=-2; //Halve the time increment andreverse its
direction continue; //Stay in the iteration }
break; //If we're right (exactly) atsunrise,
then stop the iteration } printf("%4d/%02d/%02d - %02d:%02d:%02d ",
current_date.yr,current_date.mo,current_date.da,
current_date.hh,current_date.mm,current_date.ss); current_date.hh=12;
//Reset starting date/time tonoon again current_date.mm=0;
current_date.ss=0; for (delta_t=3600;abs(delta_t)>0;) { //Iterate to
later timestarting with one hour movement
adjust_date(¤t_date,delta_t); //Use 'delta_t' to getanother
date/time if (calc_angle(¤t_date,&dangle)) { //Get the angle for this
new time printf("bad return from 'calc_angle'\n"); goto bad_exit; } if
(dangle>SUN_RISE_SET_ANGLE) { //If we're earlier than sunset if
(delta_t<0) //If we were moving toward earliertimes
delta_t/=-2; //Halve the time increment andreverse its
direction continue; //Stay in the iteration }
if (dangle<SUN_RISE_SET_ANGLE) { //If we're later than sunset if
(delta_t>0) //If we were moving toward latertimes
delta_t/=-2; //Halve the time increment andreverse its
direction continue; //Stay in the iteration }
break; //If we're right (exactly) atsunset,
then stop the iteration }
printf("%02d:%02d:%02d\n",current_date.hh,current_date.mm,current_date.ss); }
return(0);bad_exit: return(1);}//- - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/* Convert
a date and time-of-day value to the time elapsed since t=0*/int
days_in_month[12]={31,28,31,30, 31, 30, 31, 31, 30, 31, 30, 31};int
cum_days[12] ={ 0,31,59,90,120,151,181,212,242,273,303,334};int
date_to_time(datim *d,long *t) {int
max_days,cum_d,cum_sec_perihelion,cum_sec_input,cum_sec_intervening,y; if
(d->yr < 1901 || d->yr > 2099) goto bad_datim; //Century years notleap unless
divisible by 400 if (d->mo < 1 || d->mo > 12) goto bad_datim; //Month in
range? max_days=days_in_month[d->mo-1]+((!(d->yr%4) && d->mo==2) ? 1 :
0);//Account for Feb in leap year if (d->da < 1 || d->da > max_days)
goto bad_datim; //Day in range? if (d->hh < 0 || d->hh > 23) goto
bad_datim; //Hour in range? if (d->mm < 0 || d->mm > 59) goto
bad_datim; //Minute inrange? if (d->ss < 0 || d->ss > 59) goto
bad_datim; //Second inrange? if (d->tzone < -12 || d->tzone > 11) goto
bad_datim; //Time zone inrange?//Calc # seconds perihelion date/time is into
its own year cum_d=cum_days[perihelion.mo-1]+((!(perihelion.yr%4) &&
perihelion.mo>2)? 1 : 0); //Add 1 beyond Feb in leap yr
cum_d+=perihelion.da-1;//Add # days into mont cum_sec_perihelion=cum_d*86400
//Cum secs totalstarts w/secs in prior days
+(perihelion.hh-perihelion.tzone)*3600 //Plus secs inprior hrs of day (adj
for tzone) +perihelion.mm*60 //Plus secs in prior
mins +perihelion.ss; //Plus secs into current
min//Calc # seconds input date/time is into its own year
cum_d=cum_days[d->mo-1]+((!(d->yr%4) && d->mo>2) ? 1 : 0); //Add 1 beyondFeb in
leap yr cum_d+=d->da-1; //Add # days
intomonth cum_sec_input=cum_d*86400 //Cum secs
totalstarts w/secs in prior days +(d->hh-d->tzone)*3600
//Plus secs inprior hrs of day (adj for tzone) +d->mm*60
//Plus secs in prior mins +d->ss;
//Plus secs into currentmin//Calc # seconds in intervening
years starting at beginning of perihelionyear for
(cum_sec_intervening=0,y=perihelion.yr;y<d->yr;y++) //For allintervening years
cum_sec_intervening+=(365*86400)+(!(y%4) ? 86400 : 0); //Add # secs inyear,
accounting for leap *t=cum_sec_intervening-cum_sec_perihelion+cum_sec_input;
//Subtperihelion start fract, add input start fract return(0);bad_datim:
return(1);}//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - -int adjust_date(datim *dt,long
t_sec) {/* Take a date/time in 'dt' and adjust it forward (t_sec positive)
orbackward (t_sec negative) and construct a new value for 'dt'. Do not move
beyond midnight if movingbackward or 23:59:59 if moving forward. Return 0 as
a normal status. Return 1 if a limit for themovement was reached, but otherwise
set the value to the limit.*/long adj_hh,adj_mm,adj_ss; adj_hh=t_sec/3600;
adj_mm=(t_sec-adj_hh*3600)/60; adj_ss=t_sec-adj_hh*3600-adj_mm*60;
dt->ss+=adj_ss; dt->mm+=adj_mm; dt->hh+=adj_hh; if (dt->ss > 59) { dt->ss-=60;
dt->mm++; } if (dt->ss < 0) { dt->ss+=60; dt->mm--; } if (dt->mm > 59) {
dt->mm-=60; dt->hh++; } if (dt->mm < 0) { dt->mm+=60; dt->hh--; } if (dt->hh
> 23) { dt->hh=23; dt->mm=59; dt->ss=59; goto bad_exit; } if (dt->hh <
0) { dt->hh=0; dt->mm=0; dt->ss=0; goto bad_exit; } return(0);bad_exit:
return(1);}//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - -int calc_angle(datim *dt,double
*angle) {/* Given a date/time in 'dt' and using the observer's position
('obs_lat'and 'obs_lon') compute the angle between the 'earth_observer' vector
and the 'sun_observer' vector.Return its value in degrees. Return 0 if no
error. Return 1 if an error is found in the input date/time'dt'.*/long
it;double t,M,E,v,r,lon,lat,dp;VECTOR
sun_earth,obs,earth_observer,sun_observer; if (date_to_time(dt,&it)) {
//Convert date/time structure to a time't' in seconds since perihelion
printf("error converting date to time\n"); goto bad_exit; } t=it;
//Convert time to float//First calculate sun_earth vector based
on the time 't' M=2*PI*t/secs_per_orbit;
//CalculateKepler's "Mean" anomaly E=M
//ComputeKepler's "Eccentric" anomaly (approximately)
+(ECCENTRICITY-1.0/8.0*pow(ECCENTRICITY,3.0))*sin(1.0*M) //using thefirst
three terms of an infinite series +
1.0/2.0*pow(ECCENTRICITY,2.0) *sin(2.0*M) // " +
3.0/8.0*pow(ECCENTRICITY,3.0) *sin(3.0*M); // "
v=atan(sqrt((1.0+ECCENTRICITY)/(1.0-ECCENTRICITY))*tan(E/2.0))*2.0;//"Eccentric"
anomaly gives polar coord angle
r=SEMI_MAJOR_AXIS*(1.0-pow(ECCENTRICITY,2.0))/(1.0+ECCENTRICITY*cos(v));//Eq
for ellipse gives radius from angle sun_earth.x=r*cos(v); //Convert polar to
cartesian sun_earth.y=r*sin(v); // " sun_earth.z=0.0; // "//Second
calculate observer's latitude/longitude in secondary coords withelapsed time
increasing longitude lon=obs_lon+t/SIDEREAL_DAY*360.0; //Compute logitude
with time rotationdone lat=obs_lat ; //Latitude always
stays the same obs.z= 4000.0*sin(lat*PI/180.0);
//Observer's vectorfrom earth's center (secondary coords)
obs.y=(4000.0*cos(lat*PI/180.0))*sin(lon*PI/180.0); // "
obs.x=(4000.0*cos(lat*PI/180.0))*cos(lon*PI/180.0); // "//Change vector in
secondary system to vector in primary system coords earth_observer.x = a11 *
obs.x + a12 * obs.y + a13 * obs.z; //Nowtransform to primary system
earth_observer.y = a21 * obs.x + a22 * obs.y + a23 * obs.z; // "
earth_observer.z = a31 * obs.x + a32 * obs.y + a33 * obs.z; // "//Vector from
sun to observer is sum of vector from sun to earth plusobserver's vector from
earth center sun_observer.x=sun_earth.x+earth_observer.x;
sun_observer.y=sun_earth.y+earth_observer.y;
sun_observer.z=sun_earth.z+earth_observer.z;//The dot-product of the
'earth_observer' vector and 'sun_observer' vectorallows us to calculate the
angle dp=dot_product(&sun_observer,&earth_observer);
*angle=acos(dp/(vector_length(&sun_observer)*vector_length(&earth_observer)));
*angle*=180.0/PI; //Convert radians to degrees return(0);bad_exit:
return(1);}//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - -int print_vec(VECTOR *v) {
printf("%.1fx %.1fy %.1fz",v->x,v->y,v->z); return(0);}//- - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - -
-======================================================================================================
|
|