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 #include #include #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;yyr;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.yrmm+=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;ix=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_lenx*=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)/tempx==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)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); } //= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = //= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =