|
|
Hey, let's try this one more time now that I've set line-wrap to 132 columns...
/*
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 in solving the mathematics of elliptical orbits...
(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 6 hours 13
minutes and 53 seconds, the time
in an "anomalistic" orbit (perihelion to perihelion). This is 31558433
seconds. Then we use the second
equation to solve for 'E' since we now know 'M' and '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, where v=0 at perihelion)...
(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 )
The 'sun-earth' X coordinate and Y coordinate are now known as a function of
time 't'. The Z coordinate is
always zero since the orbit is in the X-Y plane.
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.5 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)
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
#define SECONDS_PER_ORBIT (365*86400+6*3600+13*60+53) //Number of seconds in a
true orbit (anomalistic) of the sun
#define ECCENTRICITY 0.0167 //Earth's elliptical
orbit eccentricity
#define AXIS_TILT 23.45 //Number of degrees of
tilt in the earth's axis
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=45.0; //Observer's latitude (positive is northern
hemispere, negative is southern hemisphere)
//double obs_lon=-120.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, our t=0)
datim s_solstice={2008, 6,20,23,59, 0,0}; //Date/time (GMT) of summer solstice
(tilt_z leans toward sun)
//Other interesting, but unneeded, dates
//datim sp_equinox={2008, 3,20, 5,48, 0,0}; //Date/time (GMT) of spring equinox
(sun_earth perpend to tilt_z)
//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)
//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 date_to_time (datim *d,long *t);
int adjust_date (datim *dt,long time);
int calc_angle (datim *dt,double *angle);
int calc_sun_earth_vector (datim *dt,VECTOR *se);
int print_vec (VECTOR *v); //Debug support
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - -
/*
Convert a date and time-of-day value to the time elapsed since perihelion
(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;
long cum_d,cum_sec_perihelion,cum_sec_input,cum_sec_intervening,y;
if (d->yr < 1901 || d->yr > 2099) goto bad_datim; //Century years
not leap 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 in range?
if (d->ss < 0 || d->ss > 59) goto bad_datim; //Second in range?
if (d->tzone < -12 || d->tzone > 11) goto bad_datim; //Time zone in
range?
//Calc # seconds perihelion date/time is into its own year
cum_d=cum_days[perihelion.mo-1] //Start with #
days in prior months
+((!(perihelion.yr%4) && perihelion.mo>2) ? 1 : 0); //Add 1 if beyond
Feb in leap yr
cum_d+=perihelion.da-1; //Add # days into
month
cum_sec_perihelion=cum_d*86400 //Cum secs total
starts w/secs in prior days
+(perihelion.hh-perihelion.tzone)*3600 //Plus secs in
prior 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] //Start with #
days in prior months
+((!(d->yr%4) && d->mo>2) ? 1 : 0); //Add 1 if beyond
Feb in leap yr
cum_d+=d->da-1; //Add # days into
month
cum_sec_input=cum_d*86400 //Cum secs total
starts w/secs in prior days
+(d->hh-d->tzone)*3600 //Plus secs in
prior hrs of day (adj for tzone)
+d->mm*60 //Plus secs in
prior mins
+d->ss; //Plus secs into
current min
//Calc # seconds in intervening years starting at beginning of perihelion year
for (cum_sec_intervening=0,y=perihelion.yr;y<d->yr;y++) //For all
intervening years
cum_sec_intervening+=(365*86400)+(!(y%4) ? 86400 : 0); //Add # secs in
year, accounting for leap
*t=cum_sec_intervening-cum_sec_perihelion+cum_sec_input; //Subt perihelion
start fract, add input start fract
return(0);
bad_datim:
printf("date_to_time: invalid datim structure (%d/%d/%d %d:%d:%d
tZone:%d\n)\n",
d->yr,d->mo,d->da,d->hh,d->mm,d->ss,d->tzone);
return(1);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - -
int main() {
int it; //Time 't' as an integer
double t; //Time in seconds from initial position of earth
VECTOR sun_earth; //Vector from sun to earth's center
VECTOR normal; //Computed vector normal to orbital plane through the sun
datim current_date; //Work area for constructing dates to be analyzed
double angle;
long delta_t;
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 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 of 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 (calc_sun_earth_vector(&s_solstice,&sun_earth)) {
printf("sunset: cannot use summer solstice date/time\n");
goto bad_exit;
}
// printf("sun_earth="); print_vec(&sun_earth); printf("\n"); //Debug
//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(AXIS_TILT*PI/180); //Z length is
'sun_earth' length divided by tan(23.45)
// printf("normal="); print_vec(&normal); printf("\n"); //Debug
//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"); //Debug
/*
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 related to the 'tilt_x'
vector since it goes from the earth's
center out the equator but at a longitude that is *not* the prime meridian.
*If* the prime meridian were
aligned exactly with this hypotenuse at the time of summer solstice then it
would be exactly noon GMT at the
summer solstice. But, instead, it is some other time (GMT) at summer
solstice (23:59 for 2008, not 12:00) and so
we know the prime meridian is rotated (on the 'tilt_z' axis) an angle of
(23:59-12:00)/SIDEREAL_DAY*360) degrees
from the orientation of the hypotenuse of our right triangle. But this is
all we need to know in order to "back
out" all the rotations of the prime meridian (which is the 'tilt_x' vector)
that have happened since the
date/time of perihelion. Then we will have the correct 'tilt_x' vector value
we want. Back to our second right
triangle. The tangent of 23.45 degrees is the ratio of the length of the
negative-Z-going normal vector to the
'sun_earth' vector's length. So the length of this negative-Z-going normal
vector is the product of the
'sun_earth' vector's length and the tangent of 23.45 degrees.
*/
normal.x=0.0; //There's no X component of any vector normal to the orbital
plane
normal.y=0.0; //Also, no Y component
normal.z=-vector_length(&sun_earth)*tan(AXIS_TILT*PI/180); //Z component
from right triangle relationship
// 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(start)="); print_vec(&tilt_x); printf("\n");
//The prime meridian is (Time at GMT)/SIDEREAL_DAY*360 degrees counterclockwise
from the "midnight" vector
it=s_solstice.hh*3600+s_solstice.mm*60+s_solstice.ss;
t=it; //Convert time to float
theta=t/SIDEREAL_DAY*2*PI;
vector_rotate(&tilt_x,&tilt_z,theta); //Rotate small amt to get prime
meridian at GMT of summer solstice
// printf("tilt_x(at GMT)="); print_vec(&tilt_x); printf("\n");
//Unwind tilt_x vector to its position at perihelion (find the angle of
clockwise rotation needed)
if (date_to_time(&s_solstice,&it)) { //Convert summer solstice date/time to
secs since perihelion
printf("sunset: unable to convert summer solstice date/time to
seconds\n");
goto bad_exit;
}
t=it; //Convert time to float
theta=-t/SIDEREAL_DAY*2*PI; //This is angle to unwind *for the specific
meridian* aligned with 'sun-earth' vector
vector_rotate(&tilt_x,&tilt_z,theta); //Unwind 'tilt_x' vector until t=0
// 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 then
calculating the 'sun_earth' vector, the
'earth_observer' vector, the 'sun_observer' vector, and the angle between
the 'sun_observer' and
'earth_observer' vectors. When the angle is close to 90 degrees then we have
a sunrise or sunset.
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 to
SUN_RISE_SET_ANGLE.
*/
for (current_date.yr=2008;current_date.yr<=2008;current_date.yr++) {
//Do several years, if wanted
printf("YEAR %d\n\n",current_date.yr);
printf(" JAN FEB MAR APR MAY JUN
JUL AUG SEP OCT
NOV DEC\n");
printf(" --------- --------- --------- --------- ---------
--------- --------- --------- --------- --------- ---------
---------\n");
for (current_date.da=1;current_date.da<=31;current_date.da++) {
//Consider 31 days in each month
printf("%2d ",current_date.da);
//Print line starts with the day #
for (current_date.mo=1;current_date.mo<=12;current_date.mo++) {
//Consider 12 months in the year
if (current_date.da>(days_in_month[current_date.mo-1]
//If day-of-month beyond max for the month
+((!(current_date.yr%4) && current_date.mo==2) ? 1 : 0))) {
//Account for Feb in leap year
printf(" ");
//Issue blanks for non-existent date
continue;
//But keep going
}
// current_date.yr=2008;
// current_date.mo=1;
// current_date.da=1;
current_date.hh=12; //Start iterative time search at noon - sun should
be definitely up!
current_date.mm=0;
current_date.ss=0;
current_date.tzone=-5;
//First, work earlier in the day looking for sunrise
if (calc_angle(¤t_date,&angle)) { //Get the angle for noon
printf("bad return from 'calc_angle'\n");
goto bad_exit;
}
if (angle<SUN_RISE_SET_ANGLE) { //If sun is below the horizon at
noon (!)
printf("error: sun below horizon at noon!\n");
goto bad_exit;
}
for (delta_t=-3600;abs(delta_t)>0;) { //Iterate to earlier time
starting with one hour movement
adjust_date(¤t_date,delta_t); //Use 'delta_t' to get another
date/time to look at
if (calc_angle(¤t_date,&angle)) { //Get the angle for this new time
printf("bad return from 'calc_angle'\n");
goto bad_exit;
}
// printf("trying %02d:%02d:%02d gives %.3f
degrees\n",current_date.hh,current_date.mm,current_date.ss,angle);
if (angle<SUN_RISE_SET_ANGLE) { //If we're earlier than sunrise
if (delta_t<0) //If we were moving toward
earlier times
delta_t/=-2; //Halve the time increment and
reverse its direction
continue; //Stay in the iteration
}
if (angle>SUN_RISE_SET_ANGLE) { //If we're later than sunrise
if (delta_t>0) //If we were moving toward later
times
delta_t/=-2; //Halve the time increment and
reverse its direction
continue; //Stay in the iteration
}
break; //If we're right (exactly) at
sunrise, then stop the iteration
}
printf("%02d%02d ",current_date.hh,current_date.mm);
current_date.hh=12; //Reset starting date/time to
noon again
current_date.mm=0;
current_date.ss=0;
for (delta_t=3600;abs(delta_t)>0;) { //Iterate to later time starting
with one hour movement
adjust_date(¤t_date,delta_t); //Use 'delta_t' to get another
date/time
if (calc_angle(¤t_date,&angle)) { //Get the angle for this new time
printf("bad return from 'calc_angle'\n");
goto bad_exit;
}
// printf("trying %02d:%02d:%02d gives %.3f
degrees\n",current_date.hh,current_date.mm,current_date.ss,angle);
if (angle>SUN_RISE_SET_ANGLE) { //If we're earlier than sunset
if (delta_t<0) //If we were moving toward
earlier times
delta_t/=-2; //Halve the time increment and
reverse its direction
continue; //Stay in the iteration
}
if (angle<SUN_RISE_SET_ANGLE) { //If we're later than sunset
if (delta_t>0) //If we were moving toward later
times
delta_t/=-2; //Halve the time increment and
reverse its direction
continue; //Stay in the iteration
}
break; //If we're right (exactly) at
sunset, then stop the iteration
}
printf("%02d%02d ",current_date.hh,current_date.mm);
}
printf("\n");
}
}
return(0);
bad_exit:
return(1);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - -
int adjust_date(datim *dt,long t_sec) {
/*
Take a date/time in 'dt' and adjust it forward (t_sec positive) or backward
(t_sec negative) and
construct a new value for 'dt'. Do not move beyond midnight if moving
backward or 23:59:59 if
moving forward. Return 0 as a normal status. Return 1 if a limit for the
movement 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,lon,lat,dp;
VECTOR sun_earth; //Vector from sun to earth's center
VECTOR earth_observer; //Vector from earth's center to the observer at given
latitude, longitude, and time
VECTOR obs; //Vector from earth's center to observer's point (after
time rotation) in *secondary coords*
VECTOR sun_observer; //Vector from sun to the observer's tangent point
if (calc_sun_earth_vector(dt,&sun_earth)) {
printf("calc_angle: unable to calculate sun_earth vector\n");
goto bad_exit;
}
if (date_to_time(dt,&it)) { //Convert date/time structure to a time 't'
in seconds since perihelion
printf("calc_angle: cannot use input date/time\n");
goto bad_exit;
}
t=it; //Convert time to float
//Second calculate observer's latitude/longitude in secondary coords with
elapsed time increasing longitude
lon=obs_lon+t/SIDEREAL_DAY*360.0; //Compute longitude with
time rotation done
lat=obs_lat ; //Latitude always stays
the same
obs.z= 4000.0*sin(lat*PI/180.0); //Observer's vector from
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; //Now transform 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 plus
observer'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' vector
allows 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 calc_sun_earth_vector(datim *dt,VECTOR *se) {
long it; //Time in seconds past perihelion (integer)
double t; //Time in seconds past perihelion (floating point)
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
int i;
if (date_to_time(dt,&it)) { //Convert date/time structure to a time 't'
in seconds since perihelion
printf("calc_sun_earth_vector: cannot use input date/time\n");
goto bad_exit;
}
t=it; //Convert time to
float
M=2*PI*t/SECONDS_PER_ORBIT; //Calculate
Kepler's "Mean" anomaly
/*
The following equation solves for E, given M, using the first 3 terms of an
infinite series. The equation
that it's solving is Kepler's relationship between the Mean anomaly M and
the Eccentric anomaly E...
E = M - e * sin E
*/
E=M //Compute
Kepler's "Eccentric" anomaly (approximately)
+(ECCENTRICITY-1.0/8.0*pow(ECCENTRICITY,3.0))*sin(1.0*M) //using the first
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); // "
/*
Another way to solve for E is iteratively using the conditions...
E[0] = M
E[i+1] = M + e * sin( E[i] )
E=M;
for (i=0;i<5;i++) { //Here's 5 iterations
E=M+ECCENTRICITY*sin(E);
}
But I've found the three terms of the infinite series works fine.
*/
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
se->x=r*cos(v); //Convert polar to cartesian
se->y=r*sin(v); // "
se->z=0.0; //There's never any Z component for earth's center
return(0);
bad_exit:
return(1);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - -
int print_vec(VECTOR *v) {
printf("%.1fx %.1fy %.1fz",v->x,v->y,v->z);
return(0);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - -
// VCTRSUPPORT.h
//= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
= = = = = = = = = = = = = = = = = =
//= = = = = = = = = = = = = = = 3D VECTOR AND GEOMETRIC PROCESSING ROUTINES = =
= = = = = = = = = = = = = = = = = =
//= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
= = = = = = = = = = = = = = = = = =
/*
Minimum angle to call lines parallel or perpendicular is arbitrarily set
to 0.1 degrees. The sine of 0.1 degrees is 0.001745.
*/
#define MIN_SINE 0.001745
#define MIN_COSINE 0.001745
//- - - - - - - - - DEFINE GEOMETRIC TYPES (POINT, VECTOR, LINE, PLANE,
TRIANGLE, etc) - - - - - - - - - - - - - -
typedef struct { //A POINT has a single definition and is defined
by...
double x,y,z; //Its 3 orthogonal (x,y,z) coords rel to the
origin
} POINT;
typedef struct { //A VECTOR has a single definition and is defined
by...
double x,y,z; //Its 3 orthogonal (x,y,z) composite vector
lengths
} VECTOR;
typedef struct { //A LINE's preferred definition
is...
POINT lp; //Any arbitrary POINT on the
line
VECTOR lv; //And a direction
VECTOR
} LINE;
/*
Note: In geometry a line has no inherit "direction" but in our "vector"
description
the VECTOR will allow us to consider the LINE to be "flowing" in one of two
directions
as well as to distinguish between the 'left' and 'right' sides of the LINE.
*/
typedef struct { //A line may also be defined
by...
POINT lp1; //Any two
points
POINT lp2; //On the
line
} line_2pt;
typedef struct { //A PLANE's preferred definition is in Hessian
form...
VECTOR pn; //A (possibly unit) VECTOR normal to the
plane
double d; //And d, the distance of the plane from the
origin
} PLANE;
typedef struct { //Or, in pure scalar (coordinate) terms, as...
double A,B,C,D; //The general equation
Ax+By+Cz+D=0
} plane_coord;
typedef struct { //Or, as...
POINT pp; //Any point on the
plane
VECTOR pv1,pv2; //And 2 non-parallel direction
vectors
} plane_pt2v;
typedef struct { //A TRIANGLE is defined preferably
by...
POINT tp; //A
POINT
VECTOR tv1,tv2; //And 2 non-parallel VECTORs which define the other two
corners
} TRIANGLE;
typedef struct { //A TRIANGLE may also be defined
by...
POINT tp1,tp2,tp3; //Three distinct non-colinear
POINTs
} triangle_3pt;
//- - - - - - - - - - - - FUNCTION PROTOTYPES OF BASIC GEOMETRIC OPERATIONS - -
- - - - - - - - -
int vector_1_point (POINT *p,VECTOR *v);
int vector_2_points (POINT *p1,POINT *p2,VECTOR *v);
double distance_2_points (POINT *p1,POINT *p2);
double vector_length (VECTOR *v);
int vector_add (VECTOR *v1,VECTOR *v2,VECTOR *v3);
int vector_sub (VECTOR *v1,VECTOR *v2,VECTOR *v3);
int vector_mpy (VECTOR *v,double scale);
int vector_resize (VECTOR *v,double vector_len);
double dot_product (VECTOR *v1,VECTOR *v2);
int cross_product (VECTOR *v1,VECTOR *v2,VECTOR *v3);
int plane_3_points (POINT *a,POINT *b,POINT *c,PLANE *pl);
int intersect_line_plane (LINE *l,PLANE *pl,POINT *p);
int plane_from_pt_2v (POINT *p,VECTOR *v1,VECTOR *v2,PLANE *pl);
int lines_2_closest_points (LINE *l1,LINE *l2,POINT *pa,POINT *pb);
double distance_2_lines (LINE *l1,LINE *l2,LINE *l3);
int intersect_2_planes (PLANE *pl1,PLANE *pl2,LINE *l);
int triangle_from_3pt (POINT *p1,POINT *p2,POINT *p3,TRIANGLE *t);
int line_2_points (POINT *p1,POINT *p2,LINE *l);
int vector_rotate (VECTOR *X,VECTOR *U,double theta);
//- - - - - - - - - - - - - FUNCTION CODE FOR BASIC GEOMETRIC OPERATIONS - - -
- - - - - - - -
/*
A VECTOR is constructed from the origin to a POINT. If the point is the
origin itself a null
vector will result, but this is not a problem. Zero is always returned.
*/
int vector_1_point(POINT *p,VECTOR *v) {
v->x=p->x;
v->y=p->y;
v->z=p->z;
return(0);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
A VECTOR is constructed from two POINTs. If the two points are the same then
a null vector
will result, but this is not a problem. Zero is always returned.
*/
int vector_2_points(POINT *p1,POINT *p2,VECTOR *v) {
v->x=p2->x-p1->x;
v->y=p2->y-p1->y;
v->z=p2->z-p1->z;
return(0);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
The distance between two points is calculated via the Pythagorean theorem.
If the points are
identical the value 0.0 will be returned so no check is needed.
*/
double distance_2_points(POINT *p1,POINT *p2) {
return sqrt((p2->x - p1->x) * (p2->x - p1->x) +
(p2->y - p1->y) * (p2->y - p1->y) +
(p2->z - p1->z) * (p2->z - p1->z));
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
The length of a VECTOR is returned. No check for a null vector is needed.
*/
double vector_length(VECTOR *v) {
return sqrt(v->x * v->x + v->y * v->y + v->z * v->z);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
Add two VECTORs v1 and v2, giving a third VECTOR v3. No check for null
vectors is needed.
Zero is always returned.
*/
int vector_add(VECTOR *v1,VECTOR *v2,VECTOR *v3) {
v3->x=v1->x+v2->x;
v3->y=v1->y+v2->y;
v3->z=v1->z+v2->z;
return(0);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
Subtract VECTOR v2 from v1, giving a third VECTOR v3. No check for null
vectors is needed.
Zero is always returned.
*/
int vector_sub(VECTOR *v1,VECTOR *v2,VECTOR *v3) {
v3->x=v1->x-v2->x;
v3->y=v1->y-v2->y;
v3->z=v1->z-v2->z;
return(0);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
Multiply the lengths of each of the three components of a VECTOR 'v' by
'scale'. The vector
'v' is modified as a result. 'scale' may be 0.0 resulting in 'v' becoming
the null vector.
This is not a problem. Zero is always returned.
*/
int vector_mpy(VECTOR *v,double scale) {
v->x*=scale;
v->y*=scale;
v->z*=scale;
return(0);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
Resize the VECTOR 'v' so its length is 'vector_len'. The vector 'v' is
modified as a result.
'vector_len' may be 0.0 resulting in 'v' becoming the null vector. This is
not a problem.
If 'vector_len' is less than zero then the direction of 'v' will be reversed
because the
sign of each of its component vectors will be reversed. In this case the
length of 'v' will
be set to the absolute value of 'vector_len'. Again, this is not a problem.
*/
int vector_resize(VECTOR *v,double vector_len) {
double factor;
/* if (vector_len<=0.0) return(0); Old logic forbidding 'vector_len' to be
negative */
factor=vector_length(v); //Get original
length of 'v'
if (factor==0.0) {
//If it's zero
if (vector_len==0.0) { //But caller actually wants
final length zero
v->x=0.0; v->y=0.0; v->z=0.0; //Then make sure they really are
absolute zero
}
else { //Otherwise, caller wants a non-zero final length, but
that's impossible
return(1); //So
return failure
}
}
else { //Otherwise, original
length was not zero
factor=vector_len/factor; //So calc scale factor (may be zero if
'vector_len' is zero)
v->x*=factor; v->y*=factor; v->z*=factor; //And apply it to all
the components
}
return(0);
//Return success
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
The dot product of the two vectors v1 and v2 is returned. If the two
vectors are *essentially* perpendicular return 0.0.
*/
double dot_product(VECTOR *v1,VECTOR *v2) {
double temp,dp;
if (v1->x==0.0 && v1->y==0.0 && v1->z==0.0 || //If
either
v2->x==0.0 && v2->y==0.0 && v2->z==0.0) { //vector is
NULL
/* printf("v=null "); */
return(0.0); //Pretend they are perpendicular and return
0.0
}
dp=v1->x * v2->x + v1->y * v2->y + v1->z * v2->z; //Compute dot
product
if (dp==0.0) { //If dot product has gone to
zero
/* printf("dp=0.0 "); */
return(0.0); //Vectors must be perpendicular so return
0.0
}
/*
Since the magnitude of the dot product of two vectors is the product of
their magnitudes times the cosine of the angle between them I can solve
for the cosine of the angle between them by dividing the already computed
dot product's magnitude by the product of their individual magnitudes. If
this cosine is close to zero then the angle must be close to 90 degrees.
If so, the vectors are 'essentially' perpendicular.
*/
temp=vector_length(v1)*vector_length(v2); //Get product of vector
lengths
if (temp==0.0) { //If product of vector lengths has gone to
zero
/* printf("v*v=null "); */
return(0.0); //Vectors must be
perpendicular
}
if (fabs(dp)/temp<MIN_COSINE) //If cosine of angle close to
zero
return(0.0); //Vectors must be
perpendicular
return(dp); //Otherwise, return normal dot
product
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
The cross product of the two VECTORS v1 and v2 is returned as a VECTOR in
v3. If either of the
two vectors are null *or* the two vectors are essentially parallel return 1
(since there is no
cross product), otherwise return 0.
*/
int cross_product(VECTOR *v1,VECTOR *v2,VECTOR *v3) {
double divisor;
if (v1->x==0.0 && v1->y==0.0 && v1->z==0.0 ||
//If either
v2->x==0.0 && v2->y==0.0 && v2->z==0.0) {
//vector is null
v3->x=0.0; //Set the output
vector to null
v3->y=0.0;
v3->z=0.0;
return(1); //Pretend parallel
and return 1
}
v3->x=v1->y * v2->z - v2->y * v1->z;
v3->y=v1->z * v2->x - v2->z * v1->x;
v3->z=v1->x * v2->y - v2->x * v1->y;
if (v3->x==0.0 && v3->y==0.0 && v3->z==0.0) { //If computed
vector is null
return(1); //The inputs must have been
parallel, return 1
}
/*
THIS OPTIONAL CODE IMPLEMENTS THE "ESSENTIALLY PARALLEL" LOGIC IF YOU WANT IT
Since the magnitude of the cross product of two vectors is the product of
their magnitudes
times the sine of the angle between them I can solve for the sine of the
angle between
them by dividing the already computed cross product's magnitude by the
product of their
individual magnitudes. If this sine is close to zero then the angle must be
close to zero
degrees. If so, the vectors are 'essentially' parallel.
divisor=vector_length(v1)*vector_length(v2); //Calc product of
input magnitudes
if (divisor==0.0) //If divisor has gone to zero input vectors too small
to work with
return(1); //Return the result as
if parallel
if ((vector_length(v3)/divisor)<MIN_SINE) //If sine of angle is
below minimum
return(1); //Return 'essentially'
parallel result
*/
return(0); //Return
normal result
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
Find the plane formed by three POINTs via function for finding a plane from
one POINT
and two VECTORS.
*/
int plane_3_points(POINT *a,POINT *b,POINT *c,PLANE *pl) {
VECTOR temp_v1,temp_v2;
temp_v1.x = b->x - a->x; //Create a VECTOR from POINT
'a' to POINT 'b'
temp_v1.y = b->y - a->y;
temp_v1.z = b->z - a->z;
temp_v2.x = c->x - a->x; //Create a VECTOR from POINT
'a' to POINT 'c'
temp_v2.y = c->y - a->y;
temp_v2.z = c->z - a->z;
if (plane_from_pt_2v(a,&temp_v1,&temp_v2,pl)==0)
//Transform to a PLANE
return(0); //If subcontract
went OK return 0
else //Otherwise, something's wrong with our input
(two points same?)
return(1);
//So return error
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
To calculate the intersection of a line and a plane use the vector
definitions of each...
_ _
n . p = d (normal vector and distance def. of PLANE)
_ _ _
p = p1 + t * l (point and vector def. of LINE)
_ _ _ _ _ _
n . (p1 + t * l) = d (subst p1 + t * l for p in PLANE equation)
_ _ _ _
n . p1 + t * n . l = d (using distributive law)
_ _ _ _
t = (d - n . p1) / (n . l) (solve for t)
_
Plug value for t into LINE eq. to get p
*/
int intersect_line_plane(LINE *l,PLANE *pl,POINT *p) {
double t,temp;
VECTOR p1v;
if ((temp=dot_product(&l->lv,&pl->pn))==0.0) //If line and plane parallel
return(1); //Return bad status 1
p1v.x=l->lp.x; p1v.y=l->lp.y; p1v.z=l->lp.z;
t=(pl->d-dot_product(&pl->pn,&p1v))/temp;
p->x = l->lp.x + t * l->lv.x;
p->y = l->lp.y + t * l->lv.y;
p->z = l->lp.z + t * l->lv.z;
return(0); //Return normal status 0
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/* Compute the Hessian form of the plane from the 2-vectors-and-point form.
_ _
n . p = d (normal vector and distance def. of PLANE)
_ _ _ _ _
n = (v1 X v2) / |v1 X v2| (def. of unit normal vector to PLANE)
_ _ _ _ _
d = ((v1 X v2) / |v1 X v2|) . p (solve for d)
*/
int plane_from_pt_2v(POINT *p,VECTOR *v1,VECTOR *v2,PLANE *pl) {
double temp;
VECTOR tempv;
if (cross_product(v1,v2,&pl->pn)==1) //If 2 vectors parallel
return(1); //Return error
temp=vector_length(&pl->pn);
pl->pn.x/=temp; pl->pn.y/=temp; pl->pn.z/=temp; //Make norm vector unit
tempv.x=p->x; tempv.y=p->y; tempv.z=p->z;
pl->d=dot_product(&pl->pn,&tempv);
return(0); //Normal return
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -
/*
Given two lines 'l1' and 'l2' determine the point 'pa' on line 'l1' and the
point 'pb' on line 'l2' that are the two closest points between the lines.
The normal return is the value 0. This occurs even when the two lines
intersect and, therefore, 'pa' and 'pb' are the same point. This is not an
error. If the given lines are parallel then there are no unique points you
can call 'closest'. In that case a 1 is returned as an error indication and
the resulting values for 'pa' and 'pb' are undefined.
*/
int lines_2_closest_points(LINE *l1,LINE *l2,POINT *pa,POINT *pb) {
VECTOR c;
PLANE pl1,pl2;
if (cross_product(&l1->lv,&l2->lv,&c)==1) { /* If lines are parallel */
return(1);
}
/*
In the next 4 function calls (plane_from_pt_2v and intersect_line_plane) the
check for an error '1' return is for safety only. Because the two input lines
have already been shown to be non-parallel these calculations *should* return
a normal result. If, though, for any reason these functions *do* return an
error result it might be because of precision problems and so we must then
return error '1' just as if our first test for parallel lines had been true.
It should never happen if we get this far.
*/
/* form plane containing line l1, parallel to vector c */
if (plane_from_pt_2v(&l1->lp,&l1->lv,&c,&pl1)==1)
return(1);
/* form plane containing line l2, parallel to vector c */
if (plane_from_pt_2v(&l2->lp,&l2->lv,&c,&pl2)==1)
return(1);
/* closest point on l1 is where l1 intersects pl2 */
if (intersect_line_plane(l1,&pl2,pa)==1)
return(1);
/* closest point on l2 is where l2 intersects pl1 */
if (intersect_line_plane(l2,&pl1,pb)==1)
return(1);
/*
Interesting, but not necessary, code for distinguishing between intersecting
or skew lines.
if (distance_2_points(pa,pb)==0.0) //Intersecting
else //Skew
*/
return(0); //Normal return
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - -
double distance_2_lines(LINE *l1,LINE *l2,LINE *l3) {
POINT p1,p2;
if (lines_2_closest_points(l1,l2,&p1,&p2)==1) //If the lines are parallel
return(-1.0); //Return an error (in future could return singular value)
l3->lp.x=p1.x; l3->lp.y=p1.y; l3->lp.z=p1.z;
l3->lv.x=p2.x-p1.x;
l3->lv.y=p2.y-p1.y;
l3->lv.z=p2.z-p1.z;
return(distance_2_points(&p1,&p2));
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - -
/*
The intersection of the two planes pl1 and pl2 is returned as a line in
point-vector format in 'l'. If the planes are parallel 1 is returned.
Otherwise, 0.
*/
int intersect_2_planes(PLANE *pl1,PLANE *pl2,LINE *l) {
VECTOR *v1,*v2;
double determinant,c1,c2,d1,d2;
if (cross_product(&pl1->pn,&pl2->pn,&l->lv)==1)
return(1);
/* OK, they're not parallel, use Paul Bourke algorithm. */
v1=&pl1->pn; d1=pl1->d;
v2=&pl2->pn; d2=pl2->d;
determinant=dot_product(v1,v1)*dot_product(v2,v2)
-dot_product(v1,v2)*dot_product(v1,v2);
/* Next check catches boundary condition where cross product not zero but
determinant is, so must still call the planes parallel. */
if (determinant==0.0)
return(1);
c1=(d1*dot_product(v2,v2)-d2*dot_product(v1,v2))/determinant;
c2=(d2*dot_product(v1,v1)-d1*dot_product(v1,v2))/determinant;
l->lp.x=c1*v1->x+c2*v2->x;
l->lp.y=c1*v1->x+c2*v2->y;
l->lp.z=c1*v1->x+c2*v2->z;
return(0);
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - -
int triangle_from_3pt(POINT *p1,POINT *p2,POINT *p3,TRIANGLE *t) {
VECTOR temp;
t->tp.x = p1->x;
t->tp.y = p1->y;
t->tp.z = p1->z;
t->tv1.x = p2->x - p1->x;
t->tv1.y = p2->y - p1->y;
t->tv1.z = p2->z - p1->z;
t->tv2.x = p3->x - p1->x;
t->tv2.y = p3->y - p1->y;
t->tv2.z = p3->z - p1->z;
if (cross_product(&t->tv1,&t->tv2,&temp)==1) //If vector from p1->p2
parallel to vector from p1->p3
return(1); //Return error, points are colinear or
two points are the same
else //Otherwise
return(0); //Normal return
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - -
/*
Return the preferred 'point-vector' definition of a line given two points
on the line. This routine always returns POINT p1 as the 'point' and the
directed vector going from POINT p1 to POINT p2 as the 'vector'. A check
is made, however, to ensure that both p1 and p2 are not the same POINT.
*/
int line_2_points(POINT *p1,POINT *p2,LINE *l) {
l->lp.x=p1->x; l->lp.y=p1->y; l->lp.z=p1->z; //LINE 'l'
defined by POINT 'p1'
l->lv.x=p2->x-p1->x; l->lv.y=p2->y-p1->y; l->lv.z=p2->z-p1->z; //And vector
p1->p2
if (l->lv.x==0.0 && l->lv.y==0.0 && l->lv.z==0.0) //If the vector is null,
points were the same
return(1); //So return bad status
else //Otherwise, points were
not the same
return(0); //Return normal status
}
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - -
/*
Rotate the VECTOR X about the VECTOR U through an angle of 'theta' radians.
Note that if X and U
are pointing in the same direction their cross product returns the null
vector. This is not a problem
because it simply causes X to be returned unchanged as it ought to.
*/
int vector_rotate(VECTOR *X,VECTOR *U,double theta) {
VECTOR temp_u,temp_x,temp_v1,temp_v2;
double temp;
temp_u=*U; //Make a copy of U so we can
resize it
vector_resize(&temp_u,1.0); //Make it a unit vector for the
following calculations
temp_x=*X; //Make a copy of X to work with
vector_mpy(&temp_x,cos(theta)); //Now have X*cos(theta)
[in temp_x]
temp=dot_product(&temp_u,&temp_x)*(1.0-cos(theta)); //Now have
(U.X)(1-cos(theta)) [in temp]
temp_v1=temp_u; //Prepare to multiply U by result
of previous line
vector_mpy(&temp_v1,temp); //Now have U*(U.X)(1-cos(theta))
[in temp_v1]
cross_product(&temp_u,X,&temp_v2); //Now have (UxX)
vector_mpy(&temp_v2,sin(theta)); //Now have (UxX)*sin(theta)
[in temp_v2]
vector_add(&temp_v1,&temp_v2,X); //Now have
U*(U.X)(1-cos(theta))+(UxX)*sin(theta) [in X]
vector_add(&temp_x,X,X); //Now have
X*cos(theta)+U*(U.X)(1-cos(theta))+(UxX)*sin(theta) [in X]
return(0);
}
//= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
= = = = = = = = = = = = = = = = = =
//= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
= = = = = = = = = = = = = = = = = =
|
|