#define Planets 15 #define NELEMS 10 // number of elements // element indices referring to PELS struct #define INC 1 // inclination #define AOP 2 // argument of perihelion or perigee #define LAN 3 // longitude of ascending node #define MAJ 4 // semi-major axis #define ECC 5 // eccentricity #define MAN 6 // mean anomaly #define TRA 7 // true anomaly #define RAD 10 // radius // Orbital velocity indices #define VANOMALY 0 #define ORBITVEL 1 #define ENOMALY 2 #define ACCEL 3 // General flags #define ON 1 #define OFF 0 // Frame names char * framename[] = {"Heliocentric","Geocentric","Equatorial","Topocentric","Cartesian", "Selenocentric","Appulse","Planetocentric"} ; // The Body name-string array char *plname[ ] = { "NO NAME ","Earth ","Moon ","Mercury ","Venus ","Mars ", "Jupiter ","Saturn ","Uranus ","Neptune ","GeNodal ","Sun ", "EqNodal ","Perigee ","Motion ","Fixed Pt" } ; // month name string array char *month[] = { "January ", "February ","March ", "April ", "May ", "June ", "July ", "August ","September ","October ", "November ","December " } ; // days in months int dimonth[ ] = { 0, // Is no month numbered zero 31, // Jan 28, // Feb 31, // Mar 30, // Apr 31, // May 30, // Jun 31, // Jul 31, // Aug 30, // Sep 31, // Oct 30, // Nov 31 // Dec } ; // element description-string array char *elname[ ] = { "0: No name ", "Inclination ", "Argument of Perihelion ", "Longitude of Asc Node ", "Semi-major Axis ", "Eccentricity ", "Mean Anomaly ", "True Anomaly ", "8: No name ", "9: No name ", "Radius " } ; // The list of bodies other than the Moon, Earth and Pluto, // currently included in the "system". // The default is None Included. int BodyList[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Geocentric lunar nodal (and Pluto) 0, // Sun 0, // Equatorial lunar nodal long or RA 0, // Perigee's longitude or RA 0, // Equatorial lunar nodal lat or dec. 0, // Geocentric Systemic Motion 0, // Fixed point 0 // Perigee } ; // list of lunar alignments included in Raster display int RastBodyList[] = { 0, // None ("None" is an unused dummy entry) 0, // Earth 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; // list of bodies for which directions are being charted or displayed int DirBody[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; // list of bodies for which accelerations are being charted or displayed int AccBody[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; // list of bodies currently showing in the viewer int ViewBody[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; int ShowShift[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; int ShowSysDir[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; int ShowMassDir[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; int ShowSysLineTP[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; int ShowInterBodyTP[] = { 0, // None ("None" is an unused dummy entry 0, // Earth to keep the indices in order.) 0, // Moon 0, // Mercury 0, // Venus 0, // Mars 0, // Jupiter 0, // Saturn 0, // Uranus 0, // Neptune 0, // Pluto (not used) 0 // Sun } ; // Bodies by program indices ordered outward from Sun. // No body at program index zero, flagged as -1 // Sun taken as zeroeth, Moon is not properly in list, but taken as highest. // Pluto not included, flagged as -1. // 0 E Mo Me Ve Ma Ju Sa Ur Ne Pl Su // 1 2 3 4 5 6 7 8 9 10 11 int bList[] = { -1, 3, 9, 1, 2, 4, 5, 6, 7, 8, -1, 0 } ; // cumulative days in months int monat[] = { 0, 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 500 } ; // Array of earth-relative masses, in same order as planet indices. // From "Norton's Star Atlas". double masses[] = { 0.0, 1.0, 0.0123001, 0.0560, 0.8150, 0.1074, 317.890, 95.140, 14.520, 17.250, 0.100, 334787.40 } ; // Array of body masses in kilograms, in same order as planet indices. // - EAR MOO MER VEN MAR JUP SAT URA NEP PLU SUN double mass_kg[] = {0.0, 5.98e24, 7.36e22, 3.18e23, 4.88e24, 6.42e23, 1.90e27, 5.68e26, 8.68e25, 1.03e26, 1.40e22, 1.99e30} ; // various initial values and constants double julianyear = 365.25, julian1900 = 2415020.0, leapyear = 366.0 ; double julian1970 = 2440587.5, ephemeris = 2415019.5 ; double twentyfour = 24.0, ezero = 23.45229444 ; double ediff = -0.35626283e-6, gzero = 99.19815967 ; double gdiff = 0.9856473354, parallax = -4.266480268e-5 ; double longitude = -7.27220522e-3, latitude = 0.902626112 ; // The above two refer to Kings Langley, Nick Thomas' home town. // It's pretty close to Greenwich, as may be seen! double mylong = -2.1490, mylat = 57.088000 ; // Aberdeen (actually Tigh a' Chomainn, Peterculter) double sylong = 151.1666666, sylat = -33.916666 ; // Sydney double stlong = -5.5715, stlat = 56.696666 ; // Strontian double tflong = -121.0333333, tflat = 39.283333 ; //'Terra Flora' CA USA double geoid = 0.9933071239 ; // epoch, per planet, // provided for cases in which the planets are not all // at same epoch for some programming reason double pe[ Planets+1 ] ; // various globals double halfpi, pi1_2, pi, quarpi, pi_4, twopi, thrhalfpi, pi3_2, r2d ; double e0, edot, en, g0, gdot, g ; double radperhr, oblong, oblat ; // Earth's rot. rate, Observer's place. double inplong, inplat ; // User-entered geographic location ; double tplong, tplat ; // General Topocentric place. double AltAz[3] ; // For Azimuth, Altitude, and Range. int sc ; // second of minute, int mn ; // minute of hour, int hr ; // hour of day, int day ; // day of month, int mo ; // month of year, and int yr ; // year, all of current epoch, as integers double x0, y_0, z0 ; // cartesian position double x1, y_1, z1 ; // co-ordinates in various frames double x2, y2, z2 ; double x3, y3, z3 ; double xva, yva, zva ; // topocentric velocity components double xv0, yv0, zv0, xv1, yv1, zv1 ; // velocity components in double xv2, yv2, zv2, xv3, yv3, zv3 ; // various frames double xa0, ya0, za0, xa1, ya1, za1 ; // acceleration components in double xa2, ya2, za2, xa3, ya3, za3 ; // various frames double ep ; // current Julian epoch. double xearth, yearth ; // Earth's heliocetric co-ordinates double vxearth, vyearth ; // Earth's heliocentric velocity components double axearth, ayearth ; // Earth's heliocentric acceleration components char dayt[] = "None " ; // default date string double xp[13], yp[13], zp[13] ; // position co-ordinates double vx[13], vy[13], vz[13] ; // velocity-vector components int First, xpl, ypl, counter = 0 ; int VELS = ON ; // 'using velocities' flag // ************************* Zero-crossings ********************************************** // Per Mean Body direction relative to Earth's Mean Body direction int wantMBDreEarZCticks = 0 ; //int MBzeroEps = 0 ; int iMBDticks[Planets] ; int firstMBDtick[Planets] ; int iMBDup[Planets] ; // upward zero-crossing count per body int iMBDdn[Planets] ; // downward zero-crossing count per body double EpMBDreEarUp[Planets][MAXALS] ; // ep per body per upward zero-crossing double EpMBDreEarDn[Planets][MAXALS] ; // ep per body per downward zero-crossing // Per Body Mass direction relative to Earth's Mass direction int wantMSDreEarZCticks = 0 ; int MDzeroEps = 0 ; int iMSDup[Planets] ; // upward zero-crossing count per body int iMSDdn[Planets] ; // downward zero-crossing count per body double EpMSDreEarUp[Planets][MAXALS] ; // ep per body per upward zero-crossing double EpMSDreEarDn[Planets][MAXALS] ; // ep per body per downward zero-crossing double epZCdiff[Planets][MAXALS] ; // zero-crossing time differences // ************************************************************************************** double Eqarray[10] ; // Struct for elements and variations, plus per-body mass // and velocity info. typedef struct tagPELS { double el[NELEMS+1] ; double var[NELEMS+1] ; double mdot ; double vel ; // velocity in a.u./year double accel ; // relative acceleration to primary in a.u./year^2 double theta ; // longitude of velocity vector double mu ; // angle to vernal equinox of velocity vector double pd ; // orbital period of body double mass ; // earth-relative mass of body } PELS, * LPPELS ; LPPELS elvar[Planets+1] ; double vop ; // a precessional velocity // Quantities pertaining to the equatorial lunar nodal typedef struct tagEqnod { double ra ; // radius of node double ta ; // true anomaly of node (wrt perigee) double ra2 ; // radius of node two double ta2 ; // true anomaly of node two (wrt perigee) double len ; // direction of equatorial ascending lunar node in range 0 to 2*pi double len1 ; // the above in range -pi/2.0 to pi/2.0 double let ; // the latitude of the node double ael ; // angle between equatorial and lunar orbital planes double vanom; // angle of velocity vector to major lunar axis double vang ; // angle of velocity vector to the equatorial nodal double alpha; // elevation of vector wrt equatorial plane double veq ; // velocity component in equatorial plane } EQNOD, * LPEQNOD ; EQNOD eqnod ; // structs for Alignment date, time, lats, longs & alignment type. typedef struct tagDATETIME { double epal ; // julian epoch double daycount ; // days from start epoch of data double alhour ; // alhour/24.0 gives fraction of day double amlong ; // moon longitude double amlat ; // moon latitude double aplong ; // body longitude double aplat ; // body latitude int alyear ; // civil year int almonth ; // civil month int alday ; // civil day int invert ; // range fix int conj ; // alignment type:1 => conj,0 => oppo char aldate[256] ; // Text-form civil date char altime[256] ; // Text-form hour } DATETIME, * LPDATETIME; LPDATETIME Date_Time[MAXALS][Planets] ; LPDATETIME modDate_Time[MAXALS][Planets] ; // Struct for lunar alignment search typedef struct tagALIGN { double MeanAlIval ; // mean IAP double epfirst ; // first epoch of search double eplast ; // final epoch of search double interval ; // of sampling double period ; // of the search int found ; // number of alignments found int centre ; // frame of reference int body ; // as it says! } ALIGN, * LPALIGN; LPALIGN align ; // General and Lunar-Axial Alignments, and // grouping by Lunar-Axial Alignment. typedef struct tagAls { ULONG nAls ; // number of alignments in grand period ULONG nAx ; // number of axial alignments in grand period ULONG nGrps ; // number of groups in grand period double EpAx[MAXALS] ; // epoch per axial alignment double EpPg[MAXALS] ; // epoch per perigeal alignment UINT Axtyp[MAXALS] ; // type per axial alignment (1 => perigeal) double AxLong[MAXALS] ; // dir. per axial alignment double grEpAx[MAXALS] ; // initial axial epoch per group double Alepof[MAXALS] ; // epoch per in-series alignment double Alongof[MAXALS] ; // long. per in-series alignment int Bodyof[MAXALS] ; // body per in-series alignment int Gridx[MAXALS] ; // index in group per in-series alignment (part-inverse of Sidx) int Grpof[MAXALS] ; // group per in-series alignment (part-inverse of Sidx) int AlsperBody[Planets] ; // count per body of in-series alignments ULONG nGr[MAXALS] ; // number of alignments per group int Sidx[MAXALS][Planets] ; // in-series index per alignment per group (inverse of Gridx) int Body[MAXALS][Planets] ; // body per alignment per group int inGrp[MAXALS][Planets] ; // alignment(s) per body per group (0 => None) ULONG Bplc[MAXALS][Planets] ; // place per al. per group double Aep[MAXALS][Planets] ; // epoch per al. per group double Along[MAXALS][Planets] ; // dir. per al. per group double alnint[MAXALS][Planets] ; double clusterality[MAXALS] ; double clusterepoch[MAXALS] ; } ALS, * LPALS ; ALS Als, modAls ; typedef struct tagPhases { int xphase[200] ; // original x data from file int yphase[200] ; // original y data from file double epPhase[200] ; // ep's of phases double Advance[200] ; // phase shifts in days (taken as zero on 1/1/1983, else negative) } PHASES, * LPPHASES ; PHASES ph ; int nphases = 0 ; #define MAXRSP 40000 // 200 by 200 lattice typedef struct tagRSP { double epo ; // epoch double xin ; // helioc. x ordinate of position in ecliptic double yin ; // helioc. y-ordinate of position in ecliptic double dir ; // direction of systemic bisectoid at position double aln ; // alignment of systemic bisectoid at position double dsp ; // dispersion of systemic bisectoid at position } RSP, * LPRSP ; LPRSP pRsp[MAXRSP] ; // pointers to FSPs // We have a line anywhere with a longitude (lo) and a latitude (la) // wrt the heliocentric/ecliptic frame. Using this struct for input // and output, we refer this line as a tangent to a given orbit // (that of 'body'), and find the two diametrically opposite points of contact. typedef struct tagTANG { int body ; // input: body of orbit double lo ; // input: heliocentric longitude of tangent double la ; // input: heliocentric latitude of tangent int quadrant ; // output double x1, y1, z1 ; // output: first contact point, heliocentric double x2, y2, z2 ; // output: second contact point, heliocentric } TANG, * LPTANG ; TANG tang ; LPTANG ptang = &tang ; /* Core functions */ void setepoch( void ) ; /* convert civil date to epoch */ void civildate( void ) ; /* epoch to civil date */ void calendar( void ) ; /* current date as string */ double quip ( double y1, /* quadratic interpolation*/ double y2, double y3 ) ; double anglemod( double a ) ; /* return 0 <= a <= 2*pi */ double range( double a ) ; /* return 0 <= a <= pi */ double range90( double a ) ; /* return 0 <= a <= pi/2.0 */ void angles( double x, double y, double z, double *L) ; /* angles from ordinates */ int LunarAlign( int ) ; void Geodirs( LPSYSTEMIC, int ) ; // get all geocentric directions void Heliodirs( LPSYSTEMIC ) ; // get all heliocentric directions void Pairdirs( LPSYSTEMIC, int ) ; // get all body-pair directions /* Point of intersection of two straight lines with A(xa,ya) and B(xb,yb) as "origins", and ua and ub as gradients */ void cross( double xa, double ya, double xb, double yb, double ua, double ub, double r[] ) ; double TanOrb(double a, // semi-major axis double e, // eccentricity double xo, // orbitocentric x of orbital point double yo // orbitocentric y of orbital point ) ; void corba( void ) ; // contact points of tangents to orbit by angle to major axis double radius( int p ) ; /* return radius of planet p */ double RADIUS( double, /* return radius from given elements */ double, double ) ; void kepler( int p ) ; /* calculate true anomaly of p, from actual elements */ double KEPLER( double, /* calculate true anomaly of p, given an epoch, the period and the eccentricity */ double, double ) ; void VELOCITY( double, /* generic orbital velocity function */ double, double, double, double * ) ; /* angular velocity of body in frame */ double PERIOD( double v, /* period of body, from eccentricity, semi-major axis, */ double a, /* eccentric anomaly and velocity. */ double e, double th ) ; void veloc( int p ) ; /* orbital velocity vector */ double AngVel( int p, /* angular velocity and fequency of body in frame */ int c, double * u) ; void heliotogeo( int p ) ; /* calculate planet ordinates */ void update( void ) ; /* update elems to epochal values */ void disp( int p ) ; /* display elems for current epoch */ double azim( double epoch, double oblo, double obla ) ; /* altazimuth: returns azimuth */ double rnext( HFILE ) ; /* read next */ void getels( void ) ; /* get the elements */ int IsLeapYr( int ) ; /* is this a leap year? (input is full 4 digit year) */ void SortBodies( LPSYSTEMIC ) ; /* sort directions into ascending order */ void Systemic( LPSYSTEMIC, /* quantities */ int /* frame */ ) ; void SystemicPrimary( LPSYSTEMIC, int ) ; // find Systemic Primary int Quadrant( double ) ; // return quadrant of angle void ConjugateBody( int bod, // identifies target body and orbit double x, // ecliptic position of reference point double y, // through which chord also passes double *j // ecliptic position of conjugate body, and ratio ) ; void IBturn(int M1, int M2, double beta, LPSYSTEMIC S) ; int GetAlignmentData( int, double *, double, double, LPFOURIER ) ; int GetSpectrum( LPFOURIER, int, double *, double, double ) ; int All_Aligns( int ) ; // alignment epochs of all included bodies, // grouped by lunar interaxial period double EqNodal( void ) ; // geocentric longitude of lunar nodal in equatorial plane void EqProps( void ) ; // various properties in relation to the equatorial plane void lperi( int ) ; // position of perihelion or of perigee double Tanomaly( double, // time from perihelion to a given anomaly, and radius. double *, int ) ; double Lanomaly( int ) ; // lunar anomaly of a body double Lanomaly2( int ) ; // lunar anomaly of a body (style 2) double Banomaly( double, // anomaly in a body's orbit of an ecliptic point. double, double *, int ) ; double Ancum( double *, // cumulative angle, function returns signed angular step double *, double * ) ; double OrbXY( double a, double e, double anomaly, double *x, double *y ) ; double RandSP( RSP *, int ) ; // systemics on a randomly chosen location in centre int Tangents( int, int, int, double, double, double, double * ) ; int GetPhases( void ) ; int NextInt( HFILE ) ; //void MassSysDironCoM( LPSYSTEMIC ) ; //void MassSysDironMB( LPSYSTEMIC ) ; //void MassSysDironBods( LPSYSTEMIC ) ; void GenSysDirandAl( LPSYSTEMIC ) ; //void SysPerIncBody( LPSYSTEMIC ) ; void MassSysDironBodsdCos( LPSYSTEMIC ) ; void SysAngVel( LPSYSTEMIC, int ) ; void GeoSys( LPSYSTEMIC ) ; void HelioSys( LPSYSTEMIC ) ; void RadiusofMeanCircle( LPSYSTEMIC ) ; void Heliosystem( LPSYSTEMIC ) ; //void SysDirectiononMB( LPSYSTEMIC, int ) ; double TPofBisectoids( LPSYSTEMIC ) ; void ModAlign( LPSYSTEMIC ) ; /********************************************************************** Altitude, Azimuth and Range of current Body. Supplied with: the current sample epoch copied from global variable ep, or epa. the location's geographic latitude and longitude. the body's co-ordinates in the equatorial frame at the sample epoch. Fills global array AltAz[]. Returns Azimuth. *********************************************************************/ double azim( double epoch, // the current sample epoch double loclo, // location's geographic longitude double locla // location's geographic latitude ) { double cs, sn, th, t, az ; double x ; // Ordinates of body, in various frames. double y ; // (These were pointers, now changed to variables. double z ; // There was probably no need to do this.) double gha ; // get the time of day, in hours, // corresponding to the current // sample ep. t = (epoch - floor(epoch)) * 24.0 ; // get the GHA directly, not via update() ; gha = anglemod( anglemod( gdot * floor(epoch) ) + g0 ) ; // get the location's hour angle. th = gha + radperhr * t + loclo ; // x0, y_0 and z0, having first been obtained // from heliotogeo(body) at the sample ep, should // correspond to the ep to which t corresponds. // rotate round earth axis cs = cos( th ) ; sn = sin( th ) ; x = x0 * cs + y_0 * sn ; // OX in ob's meridian and equ. plane, +ive outwards y = y_0 * cs - x0 * sn ; // OY in ob's position and lat. plane, +ive eastwards if (VELS) { xva = xv0 * cs + yv0 * sn ; // +ive outwards yva = yv0 * cs - xv0 * sn ; // +ive eastwards } // tilt to observer's geographic latitude cs = cos( locla ) ; sn = sin( locla ) ; z = x * cs + z0 * sn + parallax ; // OZ in zenith/nadir line. +ive zenithward x = z0 * cs - x * sn ; // OX in plane of ob's meridian, +ive northward if (VELS) { zva = xva * cs + zv0 * sn ; // +ive zenithwards xva = zv0 * cs - xva * sn ; // +ive northwards } // azimuth +ive clockwise from north point AltAz[0] = az = itan( y, x ) ; // altitude +ive zenithwards from equatorial plane AltAz[1] = atan(z / sqrt(pow(x, 2.0) + pow(y, 2.0))) ; // range AltAz[2] = sqrt(x0*x0 + y_0*y_0 + z0*z0) ; return az ; } /********************************************************************* The Civil Date as a String, dayt[]. ********************************************************************/ void calendar( void ) { char s1[17] ; static char BigDate[] = "30th September 2000 " ; int rr ; UINT bigdate ; bigdate = strlen(BigDate) ; dayt[0] = '\0' ; strcat( dayt, " " ) ; // concat leading space _itoa( day, s1, 10 ) ; // day of month strcat( dayt, s1 ) ; // concat it rr = day % 10 ; // remainder on div by 10 // sort out and concat suffix if ( rr > 3 || rr == 0 || day == 11 || day == 12 || day == 13 ) strcat( dayt, "th " ) ; if ( rr == 1 && day != 11 ) strcat( dayt, "st " ) ; if ( rr == 2 && day != 12 ) strcat( dayt, "nd " ) ; if ( rr == 3 && day != 13 ) strcat( dayt, "rd " ) ; strcat( dayt, month[mo - 1] ) ; // concat month-name _itoa( yr, s1, 10 ) ; // the year if ( strlen( s1 ) < 2 ) strcat( dayt,"0" ) ; // for two digits only strcat( dayt, s1 ) ; // concat the year while ( strlen(dayt) < bigdate ) strcat(dayt," ") ; } int IsLeapYr( int Yr ) { int leap ; leap = 0 ; // 1st assume ordinary year if (!(Yr % 4)) // If year is exactly divisible by 4. then { leap = 1 ; // it is normally a leap year, // *** BUT *** if (!(Yr % 100)) // if year is also div by 100, { if ( Yr % 400) // but not by 400, leap = 0 ; // it is not a leap year } } return leap ; } /********************************************************************* Convert the Civil Date to the Julian Date ********************************************************************/ void setepoch(void) { long epo ; int YearsFrom1900 ; double JulianDaysFrom1900 ; // At the start of a leap year, ep and civil days are 1 day different, // ep being one day more than civildays. // Correction is applied at the end of February, and thereafter ep and // civildays match. // Thus if we want an ep from a civildate, and the date falls in a leap // year, then if the date falls in January of February, the ep should be // one less than civildays, else it should be the same. YearsFrom1900 = yr - 1900 ; // these are the julian days in the above years JulianDaysFrom1900 = YearsFrom1900 * julianyear ; // cast this double to an integer: any fraction of a day // is simply truncated. epo = (long)JulianDaysFrom1900 ; // whole days if ( IsLeapYr(yr) && // and (mo==1 || mo==2) ) // if also january or february epo-- ; // reduce ep by one day. // Now add in the days corresponding to the part-year. ep = (double)( epo + day + monat[mo] ) ; calendar() ; // and obtain the string version } /****************************************************************** Convert a Julian Epoch (ep) to the Civil Date quantities (day, mo, yr) ******************************************************************/ void civildate( void ) { int x, year ; double days, hours, minutes ; yr = (int)(ep / julianyear) ; // whole yrs of julian days since 01/01/1900 year = yr + 1900 ; // extend to four digits x = (int)(ep - floor((double)yr * julianyear)) ; // whole days left over mo = 1 ; while (monat[mo + 1] < x) mo++ ; // find the nominal calendar month day = x - monat[mo] ; // find the nominal calendar day days = floor(ep) ; // whole days since 1900, irrespective of leap/non-leap years hr = (int)( hours = 24.0 * (ep - days)) ; // hour of day mn = (int)(minutes = 60.0 * (hours - (double)hr)) ; // minute of hour sc = (int)( 60.0 * (minutes - (double)mn)) ; // second of minute // Check for leap year and adjust accordingly. // At the start of a leap year, ep and civil days are 1 day different, // ep being one day more than civildays. // Correction is applied at the end of February, and thereafter ep and // civildays match. // Thus, if we want a civil date and ep falls in a leap year, then if // ep falls in january or february, civildays should be one more than // ep, else it should be the same. // So, if leap year and jan or feb, apply correction. if (IsLeapYr(year)) { if ( mo < 3 ) day++ ; if ( day == 32 ) { day = 1 ; mo = 2 ; } } yr += 1900 ; // extend to four digits calendar() ; } // ******************************************** // Fetch Elements and Variations from files. // Angles are converted from degrees to radians. // Variations are in units per day. // ******************************************** void getels( void ) { LPPELS s ; int first, last, i, j, k ; HFILE inners, outers, which ; char fn1[] = "innerels.c", fn2[] = "outerels.c" ; double ep0 ; inners = _lopen( fn1, OF_READ ) ; outers = _lopen( fn2, OF_READ ) ; for ( k = 0 ; k < 2 ; k++ ) { if ( k == 0 ) // inner planets (except Mars) { which = inners ; first = 1 ; last = 5 ; ep0 = 0.5 ; } if ( k == 1 ) // outer planets { which = outers ; first = 6 ; last = 9 ; ep0 = julian1970 - ephemeris ; } i = first ; while ( i <= last ) { s = elvar[i] ; pe[i] = ep0 ; for ( j = 1 ; j <= 3 ; j++ ) s->el[j] = rnext( which ) / r2d ; for ( j = 4 ; j <= 5 ; j++ ) s->el[j] = rnext( which ) ; s->el[MAN] = rnext( which ) / r2d ; s->mass = masses[i] ; // orbital periods from Kepler III if ( i == 2 ) // Moon s->pd = 27.55455 / julianyear ; // anomalistic month from Norton // s->pd = sqrt( masses[SUN] * exp( 3.0 * log( s->el[MAJ] ) ) ) ; else { if ( i == 1 ) // Earth s->pd = 1.0 ; else // the others s->pd = exp( 3.0*log( s->el[MAJ] )/2.0 ) ; } i++ ; } i = first ; while ( i <= last ) { s = elvar[i] ; for ( j = 1 ; j <= 3 ; j++ ) s->var[j] = rnext( which ) / r2d ; s->var[MAJ] = rnext( which ) ; s->var[ECC] = rnext( which ) ; s->var[MAN] = rnext( which ) / r2d ; s->mdot = rnext( which ) / r2d ; i++ ; } } s = elvar[SUN] ; s->pd = 0.0 ; s->theta = 0.0 ; s->vel = 0.0 ; s->mu = 0.0 ; s->mass = masses[SUN] ; _lclose( inners ) ; _lclose( outers ) ; } int GetPhases( void ) { HFILE f ; OFSTRUCT of ; int n, i, j ; double ypix1, xpix , ypix ; double eptmp, e1, e2, period, days ; char File[] = "phases1.c" ; eptmp = ep ; day = 1 ; mo = 1 ; yr = 1983 ; setepoch() ; ep1983 = e1 = ep ; day = 31 ; mo = 12 ; yr = 1999 ; setepoch() ; e2 = ep ; ep = eptmp ; period = e1 - e2 ; xpix = 377.0 ; ypix1 = 2.0 ; ypix = 636.0 ; days = -35.0 ; if ( (f = OpenFile(File, &of, OF_READ)) == HFILE_ERROR ) { msgbox = 1 ; MessageBox(NULL, "!!! Sorry, Phase Shifts Not Available !!!", "Phase Shift Fetch", MB_ICONEXCLAMATION|MB_OK) ; msgbox = 0 ; _is.PhasesExist = 0 ; return 0 ; } else { i = j = 0 ; while ((n = NextInt( f )) >= 0) { if (!i || !(i % 2)) // i is even { ph.xphase[j] = n ; i++ ; ph.Advance[j] = (double)n * days / xpix ; } else { ph.yphase[j] = n ; i++ ; ph.epPhase[j++] = period * ((double)n - ypix1) / ypix + e2 ; } } _is.PhasesExist = 1 ; _lclose( f ) ; return j ; } } int NextInt(HFILE f) { int n, i ; char c, L[10] ; i = 0 ; while (n = _lread(f,&c,1)) { if (c != 13) // if not a CR, append character L[i++] = c ; else { L[i] ='\0' ; // make ASCIIZ return (atoi(L)) ; } } return -1 ; // signal end of file } /********************************************************************** Find Cartesian Ordinates of position and velocity of body p, for the current epoch, in Heliocentric (x2, y2, z2), Geocentric (x1, y1, z1) and Equatorial (x0, y0, z0) Frames. *********************************************************************/ void heliotogeo( int p ) { LPPELS s ; double b, h, ap, ta, ra, ic, ln, ec, major, co, sn ; double m, v, ac ; s = elvar[p] ; ap = s->el[AOP] ; ta = s->el[TRA] ; ra = s->el[RAD] ; ic = s->el[INC] ; ln = s->el[LAN] ; if ( VELS ) { /* Velocities and accelerations */ m = s->mu ; /* vel ang to maj axis ("vel anomaly") */ v = s->vel ; /* the magnitude */ ac = s->accel ; /* the acceleration to primary */ if ( p > EAR && p < SUN && p != PLU ) /* skip the earth and sun - */ { /* already done in update() */ /* All else heliocentric except the moon */ h = ap + m ; /* use arg of peri */ xv3 = v * cos( h ) ; /* and vel anomaly */ h = v * sin( h ) ; yv3 = h * cos( ic ) ; /* use inclination */ zv3 = h * sin( ic ) ; co = cos(ln) ; sn = sin(ln) ; /* use long of asc node */ xv2 = xv3 * co - yv3 * sn ; yv2 = xv3 * sn + yv3 * co ; /* vel now ecliptic */ zv1 = zv2 = zv3 ; /* ecliptic zv's indifferent to centre */ // acceleration h = ap + ta + pi ; xa3 = ac * cos( h ) ; h = ac * sin( h ) ; ya3 = h * cos( ic ) ; za3 = h * sin( ic ) ; xa2 = xa3 * co - ya3 * sn ; ya2 = xa3 * sn + ya3 * co ; za1 = za2 = za3 ; } /* Earth and Sun zv's and za's not assigned, so assign them on calls! */ switch ( p ) { case EAR: /* The earth.. */ xv1 = yv1 = zv1 = 0.0 ; /* geoc: earth at rest */ xv0 = yv0 = zv0 = 0.0 ; xv2 = vxearth ; /* helioc: earth moving.. */ yv2 = vyearth ; zv2 = 0.0 ; /* ..but only in ecliptic */ xa1 = ya1 = za1 = 0.0 ; xa0 = ya0 = za0 = 0.0 ; xa2 = axearth ; ya2 = ayearth ; za2 = 0.0 ; break ; case SUN: /* The sun... */ xv2 = yv2 = zv2 = 0.0 ; /* helioc: sun at rest */ xv1 = - vxearth ; /* geoc: sun moving .. */ yv1 = - vyearth ; zv1 = 0.0 ; /* .. but only in ecliptic */ xa2 = ya2 = za2 = 0.0 ; /* helioc: sun at rest */ xa1 = - axearth ; /* geoc: sun moving .. */ ya1 = - ayearth ; za1 = 0.0 ; /* .. but only in ecliptic */ break ; case MOO: /* because moon worked geocentrically */ xv1 = xv2 ; yv1 = yv2 ; /* its geoc comps */ xv2 += vxearth ; yv2 += vyearth ; /* its helioc comps */ xa1 = xa2 ; ya1 = ya2 ; /* its geoc comps */ xa2 += axearth ; ya2 += ayearth ; /* its helioc comps */ break ; default: /* for geocentric */ xv1 = xv2 - vxearth ; yv1 = yv2 - vyearth ; xa1 = xa2 - axearth ; ya1 = ya2 - ayearth ; break ; } if ( p > EAR && p != NOD ){ /* skipped on earth & lunar node */ co = cos(en) ; sn = sin(en) ; yv0 = yv1 * co - zv1 * sn ; zv0 = yv1 * sn + zv1 * co ; xv0 = xv1 ; /* velocity now equatorial */ ya0 = ya1 * co - za1 * sn ; za0 = ya1 * sn + za1 * co ; xa0 = xa1 ; /* acceleration now equatorial */ } } /* end if ( VELS ) */ /***** Positions *****/ if ( p > EAR && p < SUN && p != NOD ) { b = ap + ta ; /* skipped for E & Sun & lunar node */ x3 = ra * cos( b ) ; b = ra * sin( b ) ; y3 = b * cos( ic ) ; z3 = b * sin( ic ) ; co = cos( ln ) ; sn = sin( ln ) ; x2 = x3 * co - y3 * sn ; y2 = x3 * sn + y3 * co ; z1 = z2 = z3 ; } switch ( p ) { case EAR: x1 = y_1 = z1 = 0.0 ; x0 = y_0 = z0 = 0.0 ; x2 = xearth ; y2 = yearth ; z2 = 0.0 ; break ; case MOO: /* because moon worked geocentrically */ x1 = x2 ; y_1 = y2 ; /* its geoc comps */ x2 += xearth ; y2 += yearth ; /* its helioc comps */ break ; case NOD: /* position of ascending lunar node */ s = elvar[MOO] ; major = s->el[MAJ] ; ec = s->el[ECC] ; ap = s->el[AOP] ; ic = s->el[INC] ; ln = s->el[LAN] ; ta = twopi - ap ; /* true anomaly of ascending node */ ra = RADIUS( major, ec, ta ) ; co = cos( ln ) ; sn = sin( ln ) ; x1 = ra * co ; y_1 = ra * sn ; /* node's geocentric components */ z1 = 0.0 ; x2 = x1 + xearth ; y2 = y_1 + yearth ; /* node's helioc comps */ z2 = z1 ; if (VELS) // nodal velocity components { double w, a ; w = s->var[LAN] ; // nodal AV in radians per day a = ln - pi/2.0 ; // directed perpendicular to radius xv1 = w * ra * cos( a ) ; // comp along 1st of Aries yv1 = w * ra * sin( a ) ; // comp in ecliptic perp. to 1st of Aries zv1 = 0.0 ; // comp perp. to ecliptic is zero xv2 = xv1 + vxearth ; // heliocentric: sun at rest,earth moving yv2 = yv1 + vyearth ; zv2 = zv1 ; co = cos(en) ; sn = sin(en) ; yv0 = yv1 * co - zv1 * sn ; zv0 = yv1 * sn + zv1 * co ; xv0 = xv1 ; // velocity now equatorial } break ; case SUN: x2 = y2 = z2 = 0.0 ; x1 = x2 - xearth ; y_1 = y2 - yearth ; z1 = 0.0 ; break ; default: x1 = x2 - xearth ; y_1 = y2 - yearth ; break ; } if ( p > EAR ) { co = cos( en ) ; sn = sin( en ) ; y_0 = y_1 * co - z1 * sn ; z0 = y_1 * sn + z1 * co ; x0 = x1 ; } } /* * Position of the ascending lunar node in * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * the equatorial plane. * ~~~~~~~~~~~~~~~~~~~~~ * * This Node lies on the cut of the equatorial plane by the Moon's * orbital plane, exactly where the moon passes northwards * through the equatorial plane. It is the equivalent in the * equatorial plane of the Ascending Lunar Node in the ecliptic. * * This function supplies the cartesian co-ordinates of this Node * in Geocentric, Heliocentric and Equatorial frames, and its * distance from the Earth, and returns the angle (varying annually) * between the Moon's orbital plane and the equatorial plane. */ double EqNodal( void ) { LPPELS h ; double xs, ys, zs ; double xs1, ys1, zs1 ; double ys2, zs2 ; double ys3, zs3 ; double xn, yn, zn ; double xn1, yn1, zn1 ; double xq, yq, zq ; double xq1, yq1, zq1 ; double yq2, zq2 ; double c, ln, rtln, ap, ra, ra2 ; double mu, ic, a, ec ; double u, wp, cs, sn ; double len, len1, let, an, ael, th, ta, ta2 ; double tmp, beta1, beta2 ; double result[5] ; int shortleg ; h = elvar[MOO] ; ln = h->el[LAN] ; // long of asc. lunar node ap = Eqarray[5] = h->el[AOP] ; // argument of asc. lunar node ic = h->el[INC] ; // inclination of lunar orbit a = h->el[MAJ] ; // semi-major axis of lunar orbit ec = h->el[ECC] ; // eccentricity of lunar orbit // We work with a plane parallel to the ecliptic and elevated // by 1.0 a.u. "above" it. The equatorial and lunar orbital planes // intersect this plane in lines parallel to the vernal equinox and // the ecliptic lunar nodal respectively. Joining the intersection // of these elevated lines to the Earth gives the line of // intersection of the equatorial and lunar orbital planes - which // must be the lunar nodal in the equatorial plane, as the lunar // orbit passes through it twice. u = 1.0 / tan(ic) ; // distance on ecliptic to elevated S(xs,ys,0) // line u is at a right angle to the ecliptic lunar nodal // We project a line through the Earth, on the lunar orbital plane and at // right angles to the ecliptic lunar nodal, to meet the elevated plane in point S. // As the inclination of the lunar orbital plane is constant, and as the // ecliptic lunar nodal rotates, taking the plane with it, the line from // Earth to S must sweep out a cone, and the locus of elevated S must // be a circle of radius u, centred at (0,0,1). It is the cut of the // cone by the elevated plane. //rtln = anglemod(ln + pi/2.0) ; // line to S is 90 degrees AHEAD of the long of ascending node rtln = ln + pi/2.0 ; // line to S is 90 degrees AHEAD of the long of ascending node xs = u * cos(rtln) ; // x in ecliptic of S ys = u * sin(rtln) ; // y in ecliptic of S zs = 1.0 ; // z of S, elevated above ecliptic mu = tan(ln) ; // slope of tangent at S to elevated circle: this // line is parallel to the ecliptic lunar nodal. c = ys - mu * xs ; // y intercept of tangent at S to elevated circle // We now project the 1st of Aries on to the elevated plane along the equatorial plane. // That is, we find the cut of the equatorial plane with the elevated plane. // Then we find the elevated point of intersection, N(xn, yn, zn), of the above tangent // (the elevated ecliptic lunar nodal) with the elevated vernal equinox. // en is the tilt of the ecliptic to the earth's equatorial plane // about the vernal equinox. It is a positive value. // // To be the tilt of the earth's equatorial plane to the ecliptic, // en must be taken as negative. yn = 1.0 / tan(-en) ; // intersection point N(xn,yn,zn) in the elevated plane xn = (yn - c) / mu ; // of the "elevated ecliptic nodal" with the "elevated zn = 1.0 ; // vernal equinox". tmp = sqrt(u*u - yn*yn) ; beta1 = itan(yn, -tmp) ; // these are the longitude boundaries, in beta2 = itan(yn, tmp) ; // the ecliptic, of the equatorial nodal. // if the line on the ecliptic through the Earth at right angles to the // ecliptic lunar nodal lies between these boundaries, then S in on the // shorter leg of its journey between A and B, else S is on the longer leg. shortleg = (rtln >= beta1 && rtln < beta2) ? 1 : 0 ; // geocentric longitude of the equatorial lunar nodal // (this is its longitude in the _ecliptic_) //len = len1 = itan(yn, xn) ; //len = len1 = atan2(yn, xn) ; len = len1 = atan(yn/xn) ; // There is a problem with convention here: when xn goes from positive to negative, // N switches instantaneously from "right" to "left", and vice versa, so we get // 180 degree angular steps as N passes through infinity. Shall we take the // equatorial nodal always to "point right", i.e., have a heading, like a ship? // The ecliptic lunar nodal conventionally has a heading: it points from Earth // through the ascending node. Should the equatorial nodal have also? // if (xn > 0) len1 -= twopi ; // if (xn < 0) len1 -= pi ; eqnod.len1 = Eqarray[0] = len1 ; // geocentric latitude of the lunar nodal in the equatorial plane //(with respect to the ecliptic) // Ok, pan to len1 /* cs = cos(len1) ; sn = sin(len1) ; xn1 = xn*cs + yn*sn ; let = itan (zn, xn1) ; */ let = atan(zn / sqrt(xn*xn + yn*yn)) ; if (xn < 0.0) let *= -1.0 ; eqnod.let = Eqarray[1] = let ; // We next want the lunar anomaly of the equatorial nodal -- // pan to the longitude of the lunar ascending node in ecliptic, and // read N(xn, yn, zn) as N(xn1, yn1, zn1) in the new frame cs = cos(ln) ; sn = sin(ln) ; xn1 = xn*cs + yn*sn ; yn1 = yn*cs - xn*sn ; zn1 = zn ; // get the argument of the equatorial nodal th = Eqarray[3] = itan(sqrt(yn1*yn1 + zn1*zn1) , xn1) ; // subtract the argument of the perigee from this to obtain // the lunar anomaly of the equatorial nodal eqnod.ta = Eqarray[4] = ta = adiff2(th, ap) ; // lunar anomaly of equatorial nodal 1 eqnod.ta2 = ta2 = anglemod(ta + pi) ; // lunar anomaly of equatorial nodal 2 // now for the ordinates of the equatorial node. eqnod.ra = ra = RADIUS(a, ec, ta ) ; // radius of equatorial node 1 eqnod.ra2 = ra2 = RADIUS(a, ec, ta2) ; // radius of equatorial node 2 wp = julianyear * h->var[AOP] ; // angular velocity of perigee in radians/j.y. VELOCITY(a, ec, ta, h->pd, result) ; // just to get mu mu = result[VANOMALY] ; // angle (to maj axis) of tangent to orbit at equatorial node vop = wp * ra * cos(ta + pi/2.0 - mu) ; // tangential velocity at node due to precession // Ouch! Another silly mistake here! //x1 = ra * cos(len1) ; // geocentric x ordinate of node //y_1 = ra * sin(len1) ; // geocentric y ordinate of node x1 = ra * cos(let) * cos(len1) ; // geocentric x ordinate of node y_1 = ra * cos(let) * sin(len1) ; // geocentric y ordinate of node z1 = ra * sin(let) ; // geocentric z ordinate of node x2 = x1 + xearth ; // heliocentric x ordinate of node y2 = y_1 + yearth ; // heliocentric y ordinate of node z2 = z1 ; // heliocentric z ordinate of node cs = cos(en) ; sn = sin(en) ; x0 = x1 ; // equatorial x ordinate of node y_0 = y_1 * cs - z1 * sn ; // equatorial y ordinate of node z0 = y_1 * sn + z1 * cs ; // equatorial z ordinate of node (should be zero!) // Find angle between lunar and equatorial planes: // Q is the point at which a line through the Earth, on the equatorial plane // at right angles to the vernal equinox, cuts the elevated plane. xq = 0.0 ; // x of Q yq = yn ; // y of Q zq = 1.0 ; // z of Q // pan to long of equatorial nodal, reading S(xs,ys,zs) and Q(xq,yq,zq) // in the elevated plane as S(xs1,ys1,zs1) and Q(xq1,yq1,zq1), respectively. cs = cos(len1) ; sn = sin(len1) ; xs1 = xs*cs + ys*sn ; ys1 = ys*cs - xs*sn ; zs1 = zs ; xq1 = xq*cs + yq*sn ; yq1 = yq*cs - xq*sn ; zq1 = zq ; // tilt to lat of nodal, re-reading S and Q. cs = cos(let) ; sn = sin(let) ; ys2 = ys1 ; zs2 = zs1*cs - xs1*sn ; yq2 = yq1 ; zq2 = zq1*cs - xq1*sn ; // find angle of Q(yq2,zq2) and rotate around // the x axis to it, reading S(ys2,zs2) an = itan(zq2, yq2) ; cs = cos(an) ; sn = sin(an) ; ys3 = ys2*cs + zs2*sn ; zs3 = zs2*cs - ys2*sn ; // get angle to 0S: this is the angle between the equatorial // and lunar planes. ael = itan(zs3, ys3) ; eqnod.ael = Eqarray[2] = ael ; return ael ; // return angle between equatorial and lunar planes } void EquProps( void ) // properties relative to the equatorial nodal { // [update() to ep is assumed done] LPPELS h ; double i, elo, ga, vm, mu, al, c, tn ; h = elvar[MOO] ; i = EqNodal() ; // angle between lunar and equatorial planes elo = eqnod.len1 ; // angle of equatorial node to 1st of Aries in current frame veloc( MOO ) ; // velocity properties of moon vm = h->vel ; // magnitude of velocity eqnod.vanom = mu = h->mu ; // angle to major axis of velocity ga = h->el[AOP] - elo ; // angle of the major axis to the equatorial nodal c = pow(cos(i), 2.0) ; tn = pow(tan(pi - (ga + mu)), 2.0) ; eqnod.alpha = al = itan(tan(i), sqrt(1.0/(c * tn) + 1.0)) ; eqnod.veq = vm * cos( al ) ; } /******************************************** Position of perihelion or of perigee ********************************************/ void lperi( int p ) // { LPPELS s ; double b, sn, cs, ln, ap, a, ra, ic, ec ; s = elvar[p] ; ap = s->el[AOP] ; ln = s->el[LAN] ; ic = s->el[INC] ; ec = s->el[ECC] ; a = s->el[MAJ] ; // ra = a * ec ; // !!!!!!! BAD mistake !!!!!!! This is (c, 0), ra = a * (1.0 - ec ) ; // and should be ((a - c), 0) x3 = ra * cos( ap ) ; b = ra * sin( ap ) ; y3 = b * cos( ic ) ; z3 = b * sin( ic ) ; // co-ordinates take the given body as primary, // so values are correct for heliocentric, but // must be adjusted for geocentric, except in // the Moon's case, which is the other way round! cs = cos( ln ) ; sn = sin( ln ) ; x2 = x3 * cs - y3 * sn ; y2 = x3 * sn + y3 * cs ; z1 = z2 = z3 ; switch ( p ) { case MOO: /* because moon is worked geocentrically */ x1 = x2 ; y_1 = y2 ; /* its geoc comps */ x2 += xearth ; y2 += yearth ; /* its helioc comps */ break ; default: /* all other bodies are worked heliocentrically */ x1 = x2 - xearth ; y_1 = y2 - yearth ; break ; } // Equatorial, all bodies. // Only the geocentric data is used here. cs = cos( en ) ; sn = sin( en ) ; y_0 = y_1 * cs - z1 * sn ; z0 = y_1 * sn + z1 * cs ; x0 = x1 ; } /********************************************************************** Obtain True Anomaly of Body p by iterated approximation *********************************************************************/ void kepler( int p ) { LPPELS s ; double e1, ea, ec, ma ; static double Limit = 1.0e-6 ; // precision s = elvar[ p ] ; // elements of body p e1 = ma = s->el[MAN] ; // mean anomaly at ep of p (also 1st guess of eccentric anomaly) ec = s->el[ECC] ; // eccentricity of p // this iteration finds the eccentric anomaly // corresponding to the mean anomaly with // the precision specified by 'Limit'. do { ea = e1 ; e1 = ma + ec * sin( ea ) ; } while ( ( e1 - ea ) > Limit ) ; // convert to true anomaly s->el[TRA] = anglemod(2.0*atan(sqrt((1+ec)/(1-ec))*tan(e1/2.0))) ; } /***************************************** Given interval from periapsis, and period and eccentricity of an orbit, find the corresponding true anomaly. *****************************************/ double KEPLER( double t, // interval double T, // period double e // eccentricity ) { double e1, ea, ma ; static double Limit = 1.0e-6 ; // precision e1 = ma = twopi * t / T ; // mean anomaly corresponding to the interval // this iteration finds the eccentric anomaly // corresponding to the mean anomaly with // the precision specified by 'Limit'. do { ea = e1 ; e1 = ma + e * sin( ea ) ; } while ( ( e1 - ea ) > Limit ) ; return anglemod(2.0*atan(sqrt((1+e)/(1-e))*tan(e1/2.0))) ; // true anomaly for interval } /********************************************************************** Quadratic Interpolation: returns Turning Point. (Assumes unit intervals on x axis, and fits parabola to three (x,y) points) *********************************************************************/ double quip ( double y_1, double y2, double y3 ) { double a, b ; a = ( y3 + y_1 )/2.0 - y2 ; b = y2 - y_1 - a ; return ( -b/(2.0*a) ) ; } /********************************************************************** Point of Intersection of straight line a with gradient ua through point A(xa,ya), with straight line b with gradient ub through point B(xb,yb). x of intersection returned in r[0] y of intersection returned in r[1] *********************************************************************/ void cross( double xa, double ya, double xb, double yb, double ua, double ub, double r[] ) { double u, v, sign ; u = ya - ua*xa - yb + ub*xb ; v = ub - ua ; if ( u == 0.0 ) { r[0] = 0.0 ; goto sweetFA ; } if ( fabs(v) < 1.0E-300 ) { sign = ( v < 0.0 ) ? -1.0 : 1.0 ; r[0] = u / (sign * 1.0E-250 ) ; } else r[0] = u / v ; /* x ordinate of intersection */ sweetFA: r[1] = ua * ( r[0] - xa ) + ya ; /* y ordinate of intersection */ } /********************************************************************* Given a set (brush) of lines all passing through the same origin, each line having its own "long" and "lat", find the long and lat of the axoid of the brush. The method is:- - find the planes in the Ox axis in which the lines lie. - find the bisectoidal plane of these planes. - find the planes in the Oy axis in which the lines lie. - find the bisectoidal plane of these planes. - find the line of intersection of these bisectoidal planes. This is the axoid. *********************************************************************/ void axoid( void ) { ; } /********************************************************************** Return Radius of Body p, using that body's elements at epoch ep. *********************************************************************/ double radius( int p ) { LPPELS s ; double major, ec, ta ; s = elvar[ p ] ; major = s->el[MAJ] ; ec = s->el[ECC] ; ta = s->el[TRA] ; return RADIUS(major, ec, ta) ; } /********************************************************************** Return Radius. given semi-major axis, eccentricity and true anomaly. This is independent of epoch, and need not pertain to an actual body, or orbit. The radius of any point in any orbit may be found. **********************************************************************/ double RADIUS( double major, double ec, double ta ) { return ( major * ( 1.0 - ec * ec )/( 1.0 + ec * cos( ta ) ) ) ; } /********************************************************************** Restrict angle to range 0 to pi *********************************************************************/ double range( double a ) { while ( a < 0.0 ) a += pi ; while ( a > pi ) a -= pi ; return ( a ) ; } /********************************************************************** Restrict angle to range 0 to pi/2.0 *********************************************************************/ double range90( double a) { double p = pi/2.0 ; while (a < 0.0 ) a += p ; while (a > p ) a -= p ; return a ; } /********************************************************************** Restrict angle to range 0 to 2*pi *********************************************************************/ double anglemod( double a ) { while ( a < 0.0 ) a += twopi ; while ( a > twopi ) a -= twopi ; return ( a ) ; } /********************************************************************** Return pointer to Long and Lat from Cartesian Position *********************************************************************/ void angles ( double x, double y, double z, double L[] ) { L[0] = itan( y, x ) ; L[1] = atan( z/(sqrt( x*x + y*y )) ) ; } /********************************************************************** Line input from ASCII, CRLF delineated file: read in chars to L[] until line-feed or end of file. If line feed, back up over the CR and substitute '\0'. Convert L[] string to Double and return it. If end of file, return -1.0. **********************************************************************/ double rnext(HFILE hf) { int cbRead, i ; char c, L[60] ; i = 0 ; while (cbRead = _lread(hf, &c, 1)) { if (c != 10) L[i++] = c ; else { L[--i] = '\0' ; // zero terminate, ASCIIZ style return (atof(L)) ; // convert to and return double } } return (-1.0) ; } /********************************************************************** Update all Elements and Variations to current epoch *********************************************************************/ void update( void ) // void update(int n, int list[]) { LPPELS s, t ; double x, ra, ve, dt ; int i, j, deltai ; en = e0 + edot * ep ; g = anglemod ( anglemod( gdot * ep ) + g0 ) ; deltai = 1 ; for ( i = EAR ; i <= MAR ; i += deltai ) { s = elvar[ i ] ; dt = ep - pe[i] ; pe[i] = ep ; for ( j = 1 ; j <= 5 ; j++ ) { s->el[j] += s->var[j] * dt ; } for ( j = 1 ; j <= 3 ; j++ ) s->el[j] = anglemod( s->el[j] ) ; x = anglemod( s->var[MAN] * dt ) ; s->el[MAN] = anglemod( s->el[MAN] + x + s->mdot * dt ) ; kepler( i ) ; s->el[RAD] = radius( i ) ; if ( VELS ) veloc( i ) ; } s = elvar[EAR] ; t = elvar[SUN] ; x = s->el[AOP] + s->el[TRA] ; ra = s->el[RAD] ; xearth = ra * cos( x ) ; yearth = ra * sin( x ) ; if ( VELS ) { x = s->el[AOP] + s->mu ; ve = s->vel ; vxearth = ve * cos( x ) ; vyearth = ve * sin( x ) ; t->vel = ve ; t->theta = s->theta ; x = s->el[AOP] + s->el[TRA] + pi ; axearth = s->accel * cos( x ) ; // heliocentric components of acceleration ayearth = s->accel * sin( x ) ; // of earth in ecliptic wrt 1st Aries } for ( j = 1 ; j <= 10 ; j++ ) { t->el[j] = s->el[j] ; t->var[j] = s->var[j] ; } t->el[AOP] += pi ; if ( t->el[AOP] > twopi ) t->el[AOP] -= twopi ; t->mdot = s->mdot ; civildate() ; for ( i = JUP ; i <= NOD ; i += deltai ) { s = elvar[i] ; dt = ep - pe[i] ; pe[i] = ep ; for ( j = 1 ; j <= 5 ; j++ ) s->el[j] += s->var[j] * dt ; x = anglemod(( s->var[MAN] + s->mdot * dt/2.0 ) * dt ) ; s->el[MAN] = anglemod( s->el[MAN] + x ) ; kepler( i ) ; s->el[RAD] = radius( i ) ; if ( VELS ) veloc( i ) ; } } // **************************************************************** // Angle wrt Major Axis of Tangent to orbit at a point in the orbit // **************************************************************** double TanOrb( double a, // semi-major axis double e, // eccentricity double xo, // orbitocentric x of orbital point double yo // orbitocentric y of orbital point ) { double c, b, y, k, th, van ; c = a * e ; // focus b = sqrt(a*a - c*c) ; // semi-minor axis k = -b/a ; //"eliptifier" y = a * yo / b ; // y ordinate on circumambient circle th = itan(y, xo) ; //"eccentric anomaly" van = itan(k, y/xo) ; // argument of tangent wrt semi-major axis if ( th < halfpi || th > thrhalfpi ) // adjust if needed van = anglemod(van + pi) ; return van ; } /************************************************************************* Find points of contact of tangents at known angles to the major axis of a body. update() for an epoch is assumed done. *************************************************************************/ void corba( void ) // contact point of tangent to orbit by angle to major axis { TANG s ; LPPELS h ; double a, b, c, e, lan, inc, aop ; double xa, ya, za ; double xb, yb, zb ; double xc, yc, zc ; double xd, yd, zd ; double beta, rho, mu, theta, cs, sn ; double xo, yo ; double x1, y1 ; double x2, y2 ; s = tang ; h = elvar[s.body] ; a = h->el[MAJ] ; e = h->el[ECC] ; c = e * a ; b = sqrt(a*a - c*c) ; rho = a / b ; lan = h->el[LAN] ; inc = h->el[INC] ; aop = h->el[AOP] ; xa = cos(s.lo) ; ya = sin(s.lo) ; za = sin(s.la) ; // pan to LAN round OZ (affects x, y) cs = cos(lan) ; sn = sin(lan) ; xb = xa * cs + ya * sn ; yb = ya * cs - xa * sn ; zb = za ; // tilt to INC round new OX (affects y, z) cs = cos(inc) ; sn = sin(inc) ; xc = xb ; yc = yb * cs + zb * sn ; zc = zb * cs - yb * sn ; // yaw to AOP round new OZ (affects x, y) cs = cos(aop) ; sn = sin(aop) ; xd = xc * cs + yc * sn ; yd = yc * cs - xc * sn ; zd = zc ; mu = itan(yd, xd) ; beta = anglemod(atan(rho * yd / xd)) ; // 0 <= beta <= twopi theta = beta - pi/2.0 ; xo = fabs(a * cos(theta)) ; yo = fabs(b * sin(theta)) ; switch (s.quadrant = Quadrant(mu)) { case 0: case 2: x1 = xo ; y1 = -yo ; x2 = -xo ; y2 = yo ; break ; case 1: case 3: x1 = xo ; y1 = yo ; x2 = -xo ; y2 = -yo ; break ; } // first contact point // unyaw by AOP round new OZ (affects x, y) cs = cos(-aop) ; sn = sin(-aop) ; xd = x1 * cs + y1 * sn ; yd = y1 * cs - x1 * sn ; zd = 0.0 ; // untilt by INC round new OX (affects y, z) cs = cos(-inc) ; sn = sin(-inc) ; xc = xd ; yc = yd * cs + zd * sn ; zc = zd * cs - yd * sn ; // unpan by LAN round OZ (affects x, y) cs = cos(-lan) ; sn = sin(-lan) ; s.x1 = xc * cs + yc * sn ; s.y1 = yc * cs - xc * sn ; s.z1 = zc ; // second contact point // unyaw by AOP round new OZ (affects x, y) cs = cos(-aop) ; sn = sin(-aop) ; xd = x2 * cs + y2 * sn ; yd = y2 * cs - x2 * sn ; zd = 0.0 ; // untilt by INC round new OX (affects y, z) cs = cos(-inc) ; sn = sin(-inc) ; xc = xd ; yc = yd * cs + zd * sn ; zc = zd * cs - yd * sn ; // unpan by LAN round OZ (affects x, y) cs = cos(-lan) ; sn = sin(-lan) ; s.x2 = xc * cs + yc * sn ; s.y2 = yc * cs - xc * sn ; s.z2 = zc ; } /************************************************************************ Obtain Velocity Vector of Body p, magnitude in a.u. per julian year, argument in range 0 to twopi with respect to the semi-major axis (zero at perihelion). Also supply a value proportional to acceleration. ************************************************************************/ void veloc( int p ) { LPPELS s ; double res[5] ; s = elvar[ p ] ; VELOCITY( s->el[MAJ], s->el[ECC], s->el[TRA], s->pd, res ) ; s->theta = res[ENOMALY] ; s->mu = res[VANOMALY] ; s->vel = res[ORBITVEL] ; s->accel = res[ACCEL] ; } /************************** Generic Velocity function. May be used for fictitious as well as kosher bodies! **************************/ void VELOCITY ( double a, // semi-major axis double e, // eccentricity double ta, // true anomaly double T, // orbital period double result[] ) { double c, b, k, r, x, y, th, u, vo ; double H, h, beta ; c = a * e ; // focus b = sqrt(a*a - c*c) ; // semi-minor axis k = -b/a ; //"eliptifier" r = RADIUS(a, e, ta) ; // radius x = c + r * cos( ta ) ; // x ordinate y = a * r * sin( ta ) / b ; // y ordinate on circumambient circle result[ENOMALY] = th = itan(y, x) ; //"eccentric anomaly" result[VANOMALY] = H = itan(k, y/x) ; // argument of velocity wrt semi-major axis if ( th < halfpi || th > thrhalfpi ) // adjust if needed result[VANOMALY] = H = anglemod(result[VANOMALY] + pi) ; u = e * cos(th) ; result[ORBITVEL] = vo = (twopi * a / T) * sqrt((1 + u)/(1 - u)) ; // the velocity in a.u. per year beta = ta + halfpi - H ; h = vo * cos(beta) ; u = 1.0 / r ; k = -e * cos(ta)/(a * (1.0 - e*e)) ; // 2nd derivative of 1/r wrt true anomaly result[ACCEL] = h*h * (u + k) ; // acceleration in a.u. per yr per yr } double PERIOD ( double v, // velocity, in a.u. per julian year double a, // semi-major axis in a.u. double e, // eccentricity double th // eccentric anomaly ) { double u, T ; u = e * cos(th) ; T = twopi * a * sqrt((1 + u) / (1 - u)) / v ; return T ; // in julian years } /************************************************** Lunar Anomaly of a Body. This is the angle, wrt the major lunar axis, of a line taken through the Earth (the focus of the lunar orbit) and the orthogonal projection of a body into the lunar orbital plane. **************************************************/ double Lanomaly(int body) { double an, lan, ic, aop, sn, cs ; double xa, ya ; double yb ; double xc, yc ; LPPELS s ; s = elvar[MOO] ; lan = s->el[LAN] ; ic = s->el[INC] ; aop = s->el[AOP] ; heliotogeo(body) ; // use the geocentric values // pan to long of ascending lunar node sn = sin(lan) ; cs = cos(lan) ; xa = x1 * cs + y_1 * sn ; ya = y_1 * cs - x1 * cs ; // tilt to inclination of lunar orbit yb = ya * cos(ic) + z1 * sin(ic) ; // pan to argument of the perigee sn = sin(aop) ; cs = cos(aop) ; xc = xa * cs + yb * sn ; yc = yb * cs - xa * sn ; an = itan(yc, xc) ; an = an ; // dummy breakpoint return an ; // the required anomaly } /************************************************** Anomaly in a Body's orbit of a point in the Ecliptic This finds the angle wrt the body's major axis of the cut of the body's orbital plane of a plane normal to the ecliptic through a body and the Body's primary focus. The cutting plane contains the "directing point" (x, y) in the ecliptic. The radius in the ecliptic of the point of intersection of the orbit is returned in * rad. **************************************************/ double Banomaly(double x, double y, double * rad, int body) { double ic, lan, aop ; double sn, cs ; double xa, ya, za ; double u, h, b, an, pg, an1 ; double xb, yb, r ; LPPELS s ; s = elvar[body] ; ic = s->el[INC] ; lan = s->el[LAN] ; aop = s->el[AOP] ; sn = sin(lan) ; cs = cos(lan) ; xa = x * cs + y * sn ; ya = y * cs - x * sn ; za = ya * tan(ic) ; u = sqrt(ya*ya + za*za) ; b = itan(u, xa) ; an = adiff2(b, aop) ; //check sense lperi(body) ; pg = itan(y2, x2) ; cs = cos(pg) ; sn = sin(pg) ; xb = x * cs + y * sn ; yb = y * cs - x * sn ; an1 = itan(yb, xb) ; h = adiff2(an1, an) ; if (h > pi/2.0 && h < 3.0*pi/2.0) an += pi ; // find radius r = RADIUS(s->el[MAJ], s->el[ECC], an) ;// orbital radius at anomaly xa = r * cos(an) ; ya = r * sin(an) ; // ordinates of r in orbit, xa along maj. axis. cs = cos(-aop) ; sn = sin(-aop) ; // pan back to lan xb = xa * cs + ya * sn ; yb = ya * cs - xa * sn ; xa = xb ; // drop on to ecliptic ya = yb * cos(ic) ; *rad = sqrt(xa*xa + ya*ya) ; // radius in ecliptic return an ; } /************************************************** Lunar Anomaly of a body, style 2. This finds the angle wrt the lunar major axis of the cut of the Lunar orbital plane of a plane normal to the ecliptic through a body and the Earth. **************************************************/ double Lanomaly2(int body) { double ant, an, sn, cs ; double mbod, mlan, ph, u, c ; double ic, lan, aop ; double pg ; // direction of perigee wrt 1st Aries in ecliptic double xp, yp ; double xp1, yp1 ; double xs, ys, zs ; double xt, yt, zt ; double xa, ya, za ; double xb, yb ; double xc, yc ; LPPELS s ; s = elvar[MOO] ; ic = s->el[INC] ; lan = s->el[LAN] ; aop = s->el[AOP] ; heliotogeo(body) ; xp = x1 ; yp = y_1 ; // geocentric ordinates of body P in ecliptic ph = itan(y_1, x1) ; mbod = tan(ph) ; // takes care of any zero divisor zs = 1.0 ; u = zs / tan(ic) ; mlan = tan(lan) ; cs = cos(lan + pi/2.0) ; // lan plus ninety degrees sn = sin(lan + pi/2.0) ; xs = u * cs ; // geocentric S(xs, ys, zs) ys = u * sn ; c = ys - mlan * xs ; // intercept xt = c / (mbod - mlan) ; // geocentric T(xt, yt, zt) yt = mbod * xt ; zt = zs ; // pan to lan, reading T cs = cos(lan) ; sn = sin(lan) ; xa = xt * cs + yt * sn ; ya = yt * cs - xt * sn ; za = zt ; // tilt to ic cs = cos(ic) ; sn = sin(ic) ; xb = xa ; yb = ya * cs + za * sn ; // was error here. had xa rather than ya. // pan to aop cs = cos(aop) ; sn = sin(aop) ; xc = xb * cs + yb * sn ; yc = yb * cs - xb * sn ; an = itan(yc, xc) ; // The above does not take proper account of the longitude of an, as // it only finds the direction wrt 1st Aries of the plane normal to the // ecliptic in which the body lies. // So get the longitude of the perigee in the ecliptic for comparison // with an. This angle is quite near to the sought longitude. lperi(MOO) ; pg = itan(y_1, x1) ; // longitude of perigee in ecliptic // pan to PG in major axis of Moon, // as projected orthogonally into the // ecliptic, reading body's position, // also as projected into the ecliptic. cs = cos(pg) ; sn = sin(pg) ; xp1 = xp * cs + yp * sn ; yp1 = yp * cs - xp * sn ; ant = itan(yp1, xp1) ; u = adiff2(ant, an) ; // an should not be very different from ant. If it is, then // an is 180 degrees out, so adjust an. if (u > pi/2.0 && u < 3.0*pi/2.0) an += pi ; return an ; } // ************************************************** // orbitocentric ordinates of a body at given anomaly // ************************************************** double OrbXY(double a, double e, double anomaly, double *x, double *y) { double c, r ; r = RADIUS(a, e, anomaly) ; // radius at anomaly c = a * e ; // focus *x = r * cos(anomaly) + c ; // x ordinate *y = r * sin(anomaly) ; // y ordinate return r ; } /************************************************** Time-interval from Perihelion of an Anomaly being the time for a body to pass from perihelion to some point in the orbit, given by the angle (the "anomaly"), wrt the major axis, of a line through that point and the primary. This is the time the body _would_ take, starting from the perihelion, with the elements held in the states they have at the instant of starting. That is, it uses "dead reckoning". It is not necessarily a time actually taken. **************************************************/ double Tanomaly(double anomaly, double * radius, int body) { double a, e, c, b, T, th, x, y, t ; LPPELS s ; s = elvar[body] ; a = s->el[MAJ] ; // semi major axis in a.u. e = s->el[ECC] ; // eccentricity T = s->pd ; // period in years // cartesian co-ordinates (in a.u.)with origin 0 // at centre of ellipse and OX positive directed // through the perihelion. * radius = OrbXY(a, e, anomaly, &x, &y) ; c = a * e ; b = sqrt(a*a - c*c) ; y = a * y / b ; th = itan(y, x) ; // the "eccentric anomaly". // the time, in days, from perihelion t = julianyear * T * (th - e * sin(th)) / twopi ; return t ; } // ************************************************************** // Apparent Angular Velocity of Body P with respect to Centre. // VELS must be on prior to update(). // ************************************************************** double AngVel(int P, int Centre, double *s) { double x, y, z, r, w ; double vx, vy, vz, v ; double phi, sn, cs ; heliotogeo( P ) ; switch (Centre) { case HE: x = x2 ; y = y2 ; z = z2 ; vx = xv2 ; vy = yv2 ; vz = zv2 ; break ; case EQ: x = x0 ; y = y_0 ; z = z0 ; vx = xv0 ; vy = yv0 ; vz = zv0 ; break ; case GE: default: x = x1 ; y = y_1 ; z = z1 ; vx = xv1 ; vy = yv1 ; vz = zv1 ; break ; } phi = itan( y, x ) ; // longitude or RA of body sn = sin( phi ) ; cs = cos( phi ) ; v = vy*cs - vx*sn ; // velocity comp. at right angles to radius. r = sqrt( x*x + y*y ) ; // effective radius of body w = v / r ; // apparent angular velocity in radians per julian year w /= julianyear ; // apparent angular velocity in radians per day. *s = w / twopi ; // frequency in cycles per day return w ; } /****************************************************** The Systemic Differential. TO DO!! "Differential" is used here in its mechanical sense. A mechanical differential has "planetary" gears interfacing two "primary" gears. The axle of the planetary gears rotates at the mean angular velocity of the interfaced gears, and its angular position is the integration wrt time of this average. The initial position of the planetary axle is arbitrary. We here extend this concept to several "primary gears". We average and integrate the angular velocities wrt Centre of the included bodies, ******************************************************/ double SysDifferential(double *av, int Centre) { return 0.0 ; // for now!! } //***************************************************** // Obtain Systemic Motion wrt Centre. // // VELS must be on prior to update(). //***************************************************** double SMotion(double *a, int Centre) { int p, start, end ; double n, s[20], ssum, csum, ph ; double vx, vy, vz ; if (Centre == HE) { start = EAR ; end = NEP ; } else { start = MOO ; end = SUN ; } for ( p = start, ssum = csum = n = 0.0 ; p <= SUN ; p++ ) { if (p == PLU || !BodyList[p]) continue ; heliotogeo( p ) ; switch (Centre) { case HE: vx = xv2 ; vy = yv2 ; vz = zv2 ; break ; case EQ: vx = xv0 ; vy = yv0 ; vz = zv0 ; break ; case GE: default: vx = xv1 ; vy = yv1 ; vz = zv1 ; break ; } s[p] = itan(vy, vx) ; ssum += sin(2.0 * s[p]) ; csum += cos(2.0 * s[p]) ; n++ ; } ph = 0.5 * itan(ssum, csum) ; *a = 0.5 * sqrt(ssum*ssum + csum*csum) / n ; return ph ; } // ********************************************* // Find Alignments of Body bod // with Moon in given interval // (the "Grand Period"). // // Starting ep is the current ep on the system // // Alignments are found in respect of a frame's "Base plane": // // Geocentric, by longitude in ecliptic or RA in the equatorial plane // Topocentric, by azimuth in the plane tangential to Earth's surface // Selenocentric, by "longitude" in the lunar orbital plane // Cartesian, by "longitude" in the cartesian regression plane // // except in the case of Appulse, // // for which the angle of closest approach is used, // ********************************************* int LunarAlign( int Centre ) { double c1, c2, c3 ; // peak-finder sequence of results double x, y, z ; // body ordinates double xm, ym, zm ; // moon ordinates double i , d, nd ; // loop control indices double lp ; // body longitude or RA double lm ; // moon longitude or RA double epa ; // epoch of an alignment double epDat1 ; // first epoch of data. double epDat2 ; // final epoch of data. double eptemp ; // used to hold current entry epoch double ept, eploop ; double s, u, hr ; // intermediates. double L[2], th ; // lat and long of alignment double value1, value2 ; int bod ; // the body under examination for alignments int hour ; // hour of day of alignment int minute ; // minute of hour of alignment int n ; // alignment counter char buff[100] = {'\0'} ; n = 0 ; align->found = 0 ; eptemp = ep ; // hold entry epoch if ( Centre != GE && Centre != EQ ) { sprintf(buff,"Bad Centre.\n\nMust be Geocentric or Equatorial") ; MessageBox(NULL, buff, "Lunar Alignment", MB_ICONEXCLAMATION|MB_OK) ; return ((UINT)0) ; } c1 = c2 = c3 = 0.0 ; // zero the peak detection quantities // Modification 14/09/2000: // // It turns out we need the IAP between the alignment _before_ the // data begins and the one after it begins, if the first data epoch // falls between two alignments (which it generally will). This is // for Agraphing purposes, so that it is known into which DAYs the // data falling before the first in-series alignment should be summed. // // We therefore back off by a number of days at least 3 larger than // the largest IAP possible (15) to be sure of catching it // (in fact we make it 5 larger), and we tack that number of days on // to the data period scanned for alignments. // // While scanning, we adjust the indexing so that the immediately // pre-series alignment has index zero, and the first in-series // alignment has index one. We still return the number of in-series // alignments found. d = align->interval ; // stepping at nd = align->period ; // over period bod = align->body ; // for body epDat1 = align->epfirst ; // first data epoch epDat2 = align->eplast ; // final data epoch ep -= 15.0 ; update() ; // back off then update to start, this ep is the loop ep for ( i = 0.0, n = 0 ; i <= nd + 35.0 ; i += d ) // d is usually 0.01 days. { if ( i > 0.0 ) { ep += d ; update() ; } // use the applied ep to find alignments heliotogeo( MOO ) ; // Get all ordinates of moon at sample ep switch ( Centre ) // Select out according to Centre { case GE: xm = x1 ; ym = y_1 ; zm = z1 ; break ; case EQ: xm = x0 ; ym = y_0 ; zm = z0 ; break ; default: break ; } heliotogeo( bod ) ; // Get ordinates of body at sample ep switch ( Centre ) // Select out according to Centre. { case GE: x = x1 ; y = y_1 ; z = z1 ; break ; case EQ: x = x0 ; y = y_0 ; z = z0 ; break ; default: break ; } lp = itan( y, x ) ; // planet longitude or RA lm = itan( ym, xm ) ; // lunar longitude or RA c1 = c2 ; c2 = c3 ; // duckshuffle old results left u = cos( lp - lm ) ; // difference of longs or RA's or azims c3 = pow( u, 2.0 ) ; // cos squared long or RA or azim diff if ( c1 < c2 && c3 < c2 && i > 2.0 * d) // is middle greater than flankers? { // if yes, we have a maximum. // refine by quadr. interp. epa = ep - d * (2.0 - quip(c1, c2, c3)) ; // this is a current date // record only the pre-series alignment // closest to the data start, in case // there is more than one. if (epa < epDat1) n = 0 ; else n++ ; // increment found-in-series count if ( n >= MAXALS ) { MessageBox(NULL, "Too Many Alignments", "LunarAlign()",MB_OK) ; return 0 ; } align->found = n ; // record index of aligns found eploop = ep ; // save current looping ep ep = epa ; // substitute ordinary alignment epoch epa update() ; // gives Elements at epa, and calls civildate(), // which gives day, mo, yr, and calls calendar(), // which gives dayt[] // get the alignment time. hour = (int)(hr = 24.0 * (epa - floor(epa))) ; if ( (minute = (int)(60.0 * (hr - hour))) == 60 ) { hour++ ; minute = 0 ; } // record the ordinary alignment epoch Date_Time[bod][n]->alyear = yr ; Date_Time[bod][n]->almonth = mo ; Date_Time[bod][n]->alday = day ; Date_Time[bod][n]->alhour = hr ; Date_Time[bod][n]->daycount = epa - epDat1 ; Date_Time[bod][n]->epal = epa ; // record julian epoch of alignment strcpy(Date_Time[bod][n]->aldate,dayt) ; // record the civil date sprintf(Date_Time[bod][n]->altime,"%2.2d:%2.2d", hour, minute) ; Date_Time[bod][n]->conj = (u > 0.0) ? 1 : 0 ; heliotogeo( MOO ) ; // get lunar ordinates at epa switch ( Centre ) { case GE: xm = x1 ; ym = y_1 ; zm = z1 ; break ; case EQ: xm = x0 ; ym = y_0 ; zm = z0 ; break ; default: break ; } angles(xm, ym, zm, L) ; // convert to long and lat or RA & dec. Date_Time[bod][n]->amlong = L[0] ; // record lunar longitude, or RA Date_Time[bod][n]->amlat = L[1] ; // record lunar latitude, or dec heliotogeo( bod ) ; // get body ordinates at epa switch ( Centre ) { case GE: x = x1 ; y = y_1 ; z = z1 ; break ; case EQ: x = x0 ; y = y_0 ; z = z0 ; break ; default: break ; } angles(x, y, z, L) ; // convert to long and lat or RA & dec. Date_Time[bod][n]->aplong = L[0] ; // record body longitude, or RA. Date_Time[bod][n]->aplat = L[1] ; // record body latitude, or dec. /********************************************************************************************* At this point we have the ordinary epoch, epa, of an alignment of body, bod, with the Moon. If we call modAlign() now, we get the systemics for this epoch, including the various orientations. We calculate the directivities according to the selected prediction model, and use it to adjust the ep, for which we obtain the elements, and civil date. This is the predicted alignment ep. **********************************************************************************************/ //if (wantModanom) { ModAlign(pSy) ; // get systemics for current ordinary alignment epoch, epa. switch (PredictionModel) // use values corresponding to model selection { case 0: value1 = pSy->mMVectDir + halfpi ; value2 = pSy->hmean_lo ; break ; case 1: value1 = pSy->mMaVectDir ; value2 = pSy->hmean_lo ; break ; case 2: value1 = pSy->hLongs[JUP] ; value2 = pSy->hmean_lo ; break ; default : value1 = pSy->hLongs[JUP] ; value2 = pSy->hmean_lo ; break ; } moda = adiff2(value1, value2) ; // difference of angles th = range90(moda/2.0 + ashift/r2d) ; // corrected and ranged s = 1.0 - pow(cos(th), 2.0) ; // get directivity ept = ep ; // hold ordinary ep ep -= Semi_lun * s ; // calculate corresponding shifted ep update() ; // substitute elements of the shifted ep for those of the current ep // get the predicted alignment time. hour = (int)(hr = 24.0 * (ep - floor(ep))) ; if ( (minute = (int)(60.0 * (hr - hour))) == 60 ) { hour++ ; minute = 0 ; } modDate_Time[bod][n]->alyear = yr ; modDate_Time[bod][n]->almonth = mo ; modDate_Time[bod][n]->alday = day ; modDate_Time[bod][n]->alhour = hr ; modDate_Time[bod][n]->daycount = ep - epDat1 ; modDate_Time[bod][n]->epal = ep ; // record julian epoch of alignment strcpy(modDate_Time[bod][n]->aldate,dayt) ; // record the date sprintf(modDate_Time[bod][n]->altime,"%2.2d:%2.2d", hour, minute) ; modDate_Time[bod][n]->invert = 0 ; modDate_Time[bod][n]->conj = (u > 0.0) ? 1 : 0 ; heliotogeo( MOO ) ; // get lunar ordinates at predicted alignment ep switch ( Centre ) { case GE: xm = x1 ; ym = y_1 ; zm = z1 ; break ; case EQ: xm = x0 ; ym = y_0 ; zm = z0 ; break ; default: break ; } angles(xm, ym, zm, L) ; // convert to long and lat or RA & dec. modDate_Time[bod][n]->amlong = L[0] ; // record lunar longitude, or RA modDate_Time[bod][n]->amlat = L[1] ; // record lunar latitude, or dec heliotogeo( bod ) ; // get body ordinates at epa switch ( Centre ) { case GE: x = x1 ; y = y_1 ; z = z1 ; break ; case EQ: x = x0 ; y = y_0 ; z = z0 ; break ; default: break ; } angles(x, y, z, L) ; // convert to long and lat or RA & dec. modDate_Time[bod][n]->aplong = L[0] ; // record body longitude, or RA. modDate_Time[bod][n]->aplat = L[1] ; // record body latitude, or dec. ep = ept ; } // if (wantModanom)... ep = eploop ; // restore current loop ep update() ; } //if ( c1 < c2 && .... // advance progress-bar if (hWndPB1 != NULL) SendMessage(hWndPB1, PBM_SETPOS, (WPARAM)( (UINT)( 10000.0 * (ep - epDat1)/(epDat2 - epDat1) ) ), 0) ; } if (hWndPB1 != NULL) // all done { SendMessage(hWndPB1, PBM_SETPOS, (WPARAM) 0, 0) ; SetWindowText(hWndFetch,"") ; } ep = eptemp ; // restore entry epoch return n ; // if zero, none found ; } /* * * Find quadrant of an angle in range 0 to twopi * * Returns zero-based index of the quadrant - that is, it * returns 0 for the 1st, 1 for the second, 2 for the third * and 3 for the fourth quadrant. * */ int Quadrant( double angle ) { angle = anglemod(angle) ; // should cope with any angle if (angle >= 0 && angle < halfpi) return 0 ; if (angle >= halfpi && angle < pi) return 1 ; if (angle >= pi && angle < thrhalfpi) return 2 ; if (angle >= thrhalfpi && angle < twopi) return 3 ; return -1 ; // on failure } /******************************** Obtain list of Geocentric Alignments of all bodies for the current grand period, grouped by inter-axial interval. ********************************/ int All_Aligns( int Centre ) { double i, nd, d ; double epa, epg, eptmp ; // ept, sum ; double lp, lm, lpg, vm, vp ; double c1[Planets], c2[Planets], c3[Planets] ; double v1[Planets], v2[Planets], v3[Planets] ; double g1, g2, g3 ; double csqsum, csum, ssum, cs, theta, twotheta ; int p ; static is_pg ; ULONG q ; // n ; ULONG gr ; // group index LPPELS s ; //memset(&Als, 0, sizeof(struct tagAls)) ; for ( p = 0 ; p < Planets ; p++) { c1[p] = v1[p] = c2[p] = v2[p] = c3[p] = v3[p] = 0.0 ; Als.AlsperBody[p] = 0 ; } g1 = g2 = g3 = 0.0 ; gr = 0L ; Als.nAls = Als.nAx = Als.nGrps = 0L ; for (q = 0L ; q < MAXALS ; q++) { Als.nGr[q] = 0L ; Als.grEpAx[q] = 0.0 ; for (p = MER ; p <= SUN ; p++) { if (p == PLU) continue ; Als.inGrp[q][p] = 0 ; Als.Bplc[q][p] = 0 ; // index of occurrence of Body in group ("Body place") } } AllAlsDone = 0 ; s = elvar[MOO] ; // elements and variations of Moon for this ep d = 0.05 ; // resolution (days) nd = AstGP ; // grand period in days for (i = 0.0 ; i <= nd ; i += d) // search loop { if ( i ) {ep += d ; update() ;} // Detect and record epochs of Moon in axis, and // the directions of the Perigee at these moments. // We use these moments to arrange the general // alignments in time-groups - that is, the // alignments of all the bodies, however ordered, // must lie on or between two lunar axial // alignments. (CAVEAT: sometimes more than one // alignment per body lies, and sometimes none lies, // between consecutive axial alignments, as the IAPs can // be greater than or less than the inter-axial period.) g1 = g2 ; g2 = g3 ; g3 = pow(cos(s->el[TRA]), 2.0) ; // we use the actual lunar anomaly // Detect peaks in cos_squared(anomaly): these occur when the anomaly // passes through 0 and 180 degrees - that is, when the Moon passes // through the major axis. Refine the estimate by quadratic // interpolation. if ( i > 2.0 * d && ( (g1 <= g2 && g2 > g3) || (g1 < g2 && g2 >= g3) ) ) { // get epoch of Moon in axis Als.grEpAx[gr] = Als.EpAx[Als.nAx] = epg = ep - d * (2.0 - quip(g1, g2, g3)) ; // get exact direction of Perigee at this epoch eptmp = ep ; ep = epg ; update() ; // detect quadrant switch (Quadrant(s->el[TRA])) { case 0: case 3: Als.Axtyp[Als.nAx] = is_pg = 1 ; // Moon is in quadrant 0 or 3 (close to perigee) break ; case 1: case 2: Als.Axtyp[Als.nAx] = is_pg = 0 ; // Moon is in quadrant 1 or 2 (close to apogee) break ; } // now calculate the longitude, or RA, of the perigee lperi(MOO) ; switch (Centre) { case EQ: lpg = itan(y_0, x0) ; // the geocentric equatorial longitude of the perigee break ; case GE: default: lpg = itan(y_1, x1) ; // the geocentric ecliptic longitude of the perigee break ; } // if this axial alignment is apogeal, add pi to the direction of the perigee if (!is_pg) {lpg += pi ; anglemod(lpg) ;} // store direction and increment count of Moon-in-Major-Axis Als.AxLong[Als.nAx++] = lpg ; // an axial alignment has been found, and so we will move on to new group, // increment group count //(alignments found before the first axial // alignment should go in group zero). //compute the temporal clustering //Als.clusterality[gr] = csqsum/Als.nGr[gr] ; //Als.clusterepoch[gr] = Semi_lun * 0.5 * itan(ssum, csum) / twopi ; gr++ ; csqsum = csum = ssum = 0.0 ; // restore search ep and corresponding elements ep = eptmp ; update() ; } heliotogeo(MOO) ; // get the Moon's direction switch (Centre) { case EQ: lm = itan(y_0, x0) ; if (wantVs && VELS) // get velocity direction vm = itan(yv0, xv0) ; break ; case GE: default: lm = itan(y_1, x1) ; if (wantVs && VELS) // get velocity direction vm = itan(yv1, xv1) ; break ; } // Check out each Body in turn, starting with Mercury. // This will be for the current group, and for the // current search ep for ( p = MER ; p <= SUN ; p++) { if (p == PLU) continue ; heliotogeo(p) ; // Get the Body's direction switch (Centre) { case EQ: lp = itan(y_0, x0) ; if (wantVs && VELS) // get velocity direction vp = itan(yv0, xv0) ; break ; case GE: default: lp = itan(y_1, x1) ; if (wantVs && VELS) // get velocity direction vp = itan(yv0, xv0) ; break ; } // peak detector c1[p] = c2[p] ; c2[p] = c3[p] ; c3[p] = pow(cos(lm - lp), 2.0) ; if (i > 2.0 * d && ( (c1[p] <= c2[p] && c2[p] > c3[p]) || (c1[p] < c2[p] && c2[p] >= c3[p]) ) ) { // an Alignment has been found: // refine estimate of its epoch by quadratic interpolation. epa = ep - d * (2.0 - quip(c1[p], c2[p], c3[p])) ; { long n ; // the zero-based ordinal index of an alignment in the series, as found long u ; // the zero-based ordinal index of an alignment in a group, as found u = Als.nGr[gr] ;// current index (== count - 1) for this group n = Als.nAls ; // current index (== count - 1) for the series Als.nGr[gr]++ ; // increment count of alignments of any kind in this group Als.Sidx[gr][u] = n ; // in-series index matching the group index of this alignment Als.Body[gr][u] = p ; // body of this alignment in this group Als.Aep[gr][u] = epa ; // epoch of this alignment in this group // temporal cluster theta = twopi * (epa-epg)/Semi_lun ; // equivalent angle cs = cos(theta) ; csqsum += pow(cs, 2.0) ; twotheta = 2.0 * theta ; csum += cos(twotheta) ; ssum += sin(twotheta) ; Als.Bplc[gr][p] = u ; // index of occurrence of Body in group ("Body place") Als.inGrp[gr][p]++ ; // there is another alignment with this body in the GROUP Als.Grpof[n] = gr ; // the group of this current in-series index Als.Gridx[n] = u ; // the matching group-index of the current in-series index Als.Bodyof[n] = p ; // body of the current in-series index Als.Alepof[n] = epa ; // the epoch of this current in-series index Als.nAls++ ; // increment total count of alignments Als.AlsperBody[p]++ ; // there is another alignment with this body in the SERIES eptmp = ep ; // hold the current search epoch ep = epa ; // update elements to the found epoch update() ; heliotogeo(p) ; // get the Body's direction at the found epoch switch (Centre) { case EQ: Als.Alongof[n] = Als.Along[gr][u] = itan(y_0, x0) ; break ; case GE: default: Als.Alongof[n] = Als.Along[gr][u] = itan(y_1, x1) ; break ; } ep = eptmp ; // restore the search epoch update() ; // and its corresponding elements } } // end of alignment found code } // for (p =... } // end of search-loop //Store count of groups Als.nGrps = gr ; AllAlsDone = 1 ; return 1 ; } double PredShift( int bod, int i, double * result ) { double mbT, mdT, dT ; mdT = modDate_Time[bod][i]->epal ; // shifted epoch of alignment mbT = Date_Time[bod][i]->epal ; // actual epoch of alignment dT = mdT - mbT ; // if shifted comes sooner than actual, dT is negative result[0] = dT / Semi_lun + 1.0 ; result[1] = dT ; return (mbT + dT) ; // return predicted alignment epoch } //********************************************************** //* Collect systemic data for Fourier and other analysis * //********************************************************** //********************************************************** //********************************************************** //* Get Spectrum of Alignment "energy" for comparison * //* with corresponding Lambda "energy" -- and a few * //* other things! * //********************************************************** int GetAlignmentData( int Centre, double * res, double hi_order, double lo_order, LPFOURIER fou ) { double i, d, nd ; // loop control double lp[20] ; // body longitude or RA or azimuth double lm ; // moon longitude or RA or azimuth double AVbod[20], rate ; // angular velocity double mxclr, mnclr, clr ; //double cs ; //double sn ; double x, y, z, xm, ym, zm ; // for appulse calculation //double slop, sloe ; // for selenocentric calculation double lan, inc ; LPPELS Q ; int bod ; static room, was ; static c, j, k, h ; static hi_coflo = 0, hi_cofhi = 0 ; static long sum = 0L, prod = 0L ; static double oT = 0.0 ; // Most recent implied period static double oep ; // Most recent ep static double N, csum, ssum, rsum ; static double xmean, ymean ; // heliocentric cartesian Mean Body static double bx[20], by[20] ; // heliocentric cartesian positions of bodies static double lb[20] ; // heliocentric dirs of bodies static double sb[20] ; // dirs on mean body Q = elvar[MOO] ; // for selenocentric. bod = align->body ; // for body. nd = align->period ; // over grand period. room = MAXDATA - 1 ; // room available. d = nd / (double)room ; // sample as room allows--- // we should get a nearly constant MAXDATA samples if (hi_order > (double)room / 2.0) // keep top order within bounds { msgbox = 1 ; MessageBox(NULL,"High order bound adjusted.", "High Order Bound", MB_OK) ; msgbox = 0 ; fou->hiAlOrd = hi_order = (double)room / 2.0 ; * res = (hi_order - lo_order) / align->period ; } else fou->hiAlOrd = hi_order ; fou->loAlOrd = lo_order ; fou->OrderStep = * res ; if (pIMFou != NULL) { pIMFou->hiAlOrd = hi_order ; pIMFou->loAlOrd = lo_order ; pIMFou->OrderStep = * res ; } for ( // the business! i = 0.0, mnclr = 1.0e10, // min closure rate mxclr = -1.0e10, // max closure rate j = 0 ; i <= nd ; i += d ) { if ( i == 0.0 ) oep = ep ; if ( i > 0.0 ) { ep += d ; update() ; } // this is always the "current, physical" epoch if (wantModanom) { double u, th, ept, value1, value2 ; ModAlign(pSy) ; // calculate systemics for current, "physical" epoch switch (PredictionModel) // choose from among them { case 0: value1 = pSy->mMVectDir + halfpi ; value2 = pSy->hmean_lo ; break ; // dir of velocity case 1: value1 = pSy->mMaVectDir ; value2 = pSy->hmean_lo ; break ; // dir of acceleration case 2: value1 = pSy->hLongs[JUP] ; value2 = pSy->hmean_lo ; break ; // dir of jupiter default : value1 = pSy->hLongs[JUP] ; value2 = pSy->hmean_lo ; break ; } moda = adiff2(value1, value2) ; // this is for the current ep th = range90(moda / 2.0 + ashift/r2d) ; // corrected and ranged u = 1.0 - pow(cos(th), 2.0) ; // get directivity ept = ep ; // hold current ep ep -= Semi_lun * u ; // calculate corresponding shifted ep update() ; // substitute elements for the shifted ep for those of the current ep ep = ept ; // restore current ep } lan = Q->el[LAN] ; inc = Q->el[INC] ; heliotogeo( MOO ) ; // Get all ordinates of moon at epoch ep switch ( Centre ) // Select out according to Centre { double u ; case GE: lm = itan( ym = y_1, xm = x1) ; // Geocentric lunar longitude zm = z1 ; if (VELS) {AVbod[MOO] = AngVel(MOO, GE, &u) ;} break ; case EQ: lm = itan(ym = y_0, xm = x0) ; // lunar RA zm = z0 ; if (VELS) {AVbod[MOO] = AngVel(MOO, EQ, &u) ;} break ; default: break ; } // Riemann Sum of Geocentric Systemic Directivity // I'm not yet using all this!! //for (h = MER ; h <= SUN ; h++) // not Earth, not Moon //{ // if(h == PLU) continue ; // heliotogeo( h ) ; // Get ordinates of body at epoch ep heliotogeo( bod ) ; // Get ordinates of body at epoch ep switch ( Centre ) // Select out according to Centre. { double u ; case GE: lp[bod] = itan( y = y_1, x = x1 ) ; // planet longitude z = z1 ; if (VELS) AVbod[bod] = AngVel(bod, GE, &u) ; break ; case EQ: lp[bod] = itan( y = y_0, x = x0 ) ; // planet RA z = z0 ; if (VELS) AVbod[bod] = AngVel(bod, EQ, &u) ; break ; default: break ; } // the angular difference, and the cos-squared alignment itself, defn I //if (h == bod) { double k ; k = pow(cos(lm - lp[bod]), 2.0) ; if (wantModanom) fou->moddata[j] = k ; else fou->data[j] = k ; } //} // for (h = MER... if (VELS) { /* * "Closure rate" is just the difference between the angular velocities * (viewed geocentrically and in the ecliptic, equatorial plane or, in * the topocentric case, the plane tangential to the earth) of the Moon * and the selected body. If the rate were exactly constant, the sine or * cosine of the position-difference would be a perfect sinusoid with a * constant period, T = twopi/AV. * * The rate is not constant, so we may speak of an "instantaneous" period, * i.e., the period corresponding to a particular value of velocity-difference * at a particular instant. For want of a briefer term, we call this period * the "implied" period. * * Given this, we may reckon out what, were it sustained, its starting, * middle and ending epochs would be. Were a bud "sensitive" to rate * in some way, it conceivably might compile a sort of dossier * of such significant moments, and use it to predict such moments. * * The following code represents an attempt to ascertain what information * could be available from such a compilation, and (eventually, perhaps) * whether it could account for the "phase shift". * * The implied periods are calculated,"cluster histograms" of the * sigificant moments, and a general histogram of these moments, covering * the grand period, are compiled. */ fou->avbod[j] = AVbod[bod] ; // AV of selected body fou->avlun[j] = AVbod[MOO] ; // AV of Moon clr = AVbod[MOO] - AVbod[bod] ; // difference of AVs fou->angclr[j] = clr ; // closure rate mxclr = (clr >= mxclr)? clr : mxclr ; // max rate mnclr = (clr <= mnclr)? clr : mnclr ; // min rate fou->maxclr = mxclr ; fou->minclr = mnclr ; // time rate of change of alignment fou->dAbydt[j] = rate = - clr * sin(2.0 * (lm - lp[bod])) ; /*** if (wanttImps && !_is.ImpsExist) { double T, Tm, E, e, ip, angdiff, rimpt ; int Ns, Ne ; imp->lunipd[j] = Tm = twopi / AVbod[MOO] ; // rate-implied period of moon fou->clrpd[j] = T = twopi / clr ; // rate-implied period of moon & body if (wanttImps && pIMFou != NULL) pIMFou->data[j] = T ; // for spectrum // greatest and least implied periods of moon and body imp->maxImpPd = (T > imp->maxImpPd)? T : imp->maxImpPd ; imp->minImpPd = (T < imp->minImpPd)? T : imp->minImpPd ; // compile general implied period histogram e = 6.0 * T ; // work to 4-hour class intervals c = (modf(e, &ip) < 0.5) ? (int)ip : (int)ip + 1 ; // sort out indices imp->lo_c = c < imp->lo_c ? c : imp->lo_c ; // lowest class interval index imp->hi_c = c >= imp->hi_c ? c : imp->hi_c ; // highest class interval index imp->Timp[c]++ ; // accumulate implied period instances imp->MxTimps = // find most populated class interval imp->Timp[c] > imp->MxTimps ? imp->Timp[c] : imp->MxTimps ; imp->Ttot += imp->Timp[c] ; // total of implied period instances // Cluster histograms // For this, the j'th, sample.... angdiff = adiff2(lm, lp[bod]) ; // absolute angular position difference rimpt = angdiff / clr ; // rate-implied time into cycle E = ep - rimpt ; // start ep of rate-implied period // Number of cycles of period T from the start of the grand period Ns = (int)((E - align->epfirst) / T) ; // Number of cycles of period T to the end of the grand period Ne = (int)(( align->eplast - E) / T) ; // at the start of a cycle, cosine is +1 ("hi"), representing implied conjunction e = imp->epnexthi[j] = E - T ; // start ep of rate-implied period (the one before the current) e -= align->epfirst ; // sort out indices if ( e >= 0.0 && e <= align->eplast) { e *= 6.0 ; // work to 4-hour class intervals c = (modf(e, &ip) < 0.5) ? (int)ip : (int)ip + 1 ; // sort out indices hi_cofhi = c > hi_cofhi ? c : hi_cofhi ; imp->epClusterHi[c]++ ; // build cluster histo } // half way through a cycle, cosine is -1 ("lo"), representing implied opposition e = imp->epnextlo[j] = E - T / 2.0 ; // middle ep of rate-implied period //(the one before the current) e -= align->epfirst ; // sort out indices if ( e >= 0.0 && e <= align->eplast) // keep in range { e *= 6.0 ; // work to 4-hour class intervals c = (modf(e, &ip) < 0.5) ? (int)ip : (int)ip + 1 ; // sort out indices hi_coflo = c > hi_coflo ? c : hi_coflo ; imp->epClusterLo[c]++ ; // build cluster histo } // at the start of a cycle, cosine is +1 ("hi"), representing implied conjunction e = imp->epnexthi[j] = E ; // start ep of rate-implied period (the current one) e -= align->epfirst ; // sort out indices if ( e >= 0.0 && e <= align->eplast) // keep in range { e *= 6.0 ; // work to 4-hour class intervals c = (modf(e, &ip) < 0.5) ? (int)ip : (int)ip + 1 ; // sort out indices hi_cofhi = c > hi_cofhi ? c : hi_cofhi ; imp->epClusterHi[c]++ ; // build cluster histo } // half way through a cycle, cosine is -1 ("lo"), representing implied opposition e = imp->epnextlo[j] = E + T / 2.0 ; // middle ep of rate-implied period (the current one) e -= align->epfirst ; // sort out indices if ( e >= 0.0 && e <= align->eplast) { e *= 6.0 ; // work to 4-hour class intervals c = (modf(e, &ip) < 0.5) ? (int)ip : (int)ip + 1 ; // sort out indices hi_coflo = c > hi_coflo ? // highest index of low so far c : hi_coflo ; imp->epClusterLo[c]++ ; // build cluster histo } // at the start/end of a cycle, cosine is +1 ("hi"), representing implied conjunction e = imp->epnexthi[j] = E + T ; // start ep of rate-implied period // (end of the current one, and start of the next one) e -= align->epfirst ; // sort out indices if ( e >= 0.0 && e <= align->eplast) { e *= 6.0 ; // work to 4-hour class intervals c = (modf(e, &ip) < 0.5) ? (int)ip : (int)ip + 1 ; // sort out indices hi_cofhi = c > hi_cofhi ? // highest index of high so far c : hi_cofhi ; imp->epClusterHi[c]++ ; // build cluster histo } } // end of if (wanttImps && !_is.ImpsExist) ***/ } // if (VELS) j++ ; // next sample // advance the progress-bar if (hWndPB1 != NULL) SendMessage(hWndPB1, PBM_SETPOS, (WPARAM)( (UINT)(10000.0 * i / nd) ), 0) ; } // for (i --- End of Sampling fou->items = j ; // There were j samples // This is always done, if VELS are on, // so mark as existing, if they are. // It should be cancelled on a change of // grand period dates, or of the body for // which alignments are found. if (VELS) _is.AvelsExist = 1 ; // This is always done, so mark as existing. // It should be cancelled on a change of // grand period dates, or of the body for // which alignments are found. _is.CsqRateExists = 1 ; if (pIMFou != NULL) pIMFou->items = j ; if (hWndPB1 != NULL) // all done { SendMessage(hWndPB1, PBM_SETPOS, (WPARAM) 0, 0) ; SetWindowText(hWndFetch,"") ; } /**** if (VELS && wanttImps && !_is.ImpsExist) // Do implied-period processing { // find weighted mean of all implied periods for ( c = imp->lo_c, sum = 0L, prod = 0L ; c <= imp->hi_c ; c++ ) { UINT n ; sum += (n = imp->Timp[c]) ; prod += (n * c) ; } // this is an INDEX value!! imp->wmTimp = (double)prod/(double)sum ; // isolate and count clusters, then find their // weighted mean ep's. // This is NOT foolproof, as it depends on clusters // actually being isolated (i.e. by class intervals // of zero content), which they sometimes are // not. Better method needed. Could we associate // clusters with the actual alignments....? was = imp->bigloidx = 0 ; for ( k = 0 ; k <= hi_coflo ; k++ ) { if (imp->epClusterLo[k] > 0) // is non zero now { if (was == 0) // if was zero { was = 1 ; // we start a cluster imp->locluststart[imp->bigloidx] = k ; // on this one } } else // is zero now { if (was) // if was non zero { was = 0 ; // we end a cluster imp->loclustend[imp->bigloidx++] = k - 1 ; // on the one before } } } // we should have the lo clusters isolated now. // find their weighted-mean ep's for ( k = 0 ; k <= imp->bigloidx ; k++ ) { for ( c = imp->locluststart[k], sum = 0L, prod = 0L ; c <= imp->loclustend[k] ; c++ ) { UINT n ; sum += (n = imp->epClusterLo[c]) ; prod += n * c ; } imp->Lomeanep[k] = (double)prod/(double)sum ; } was = imp->bighiidx = 0 ; for ( k = 0 ; k <= hi_cofhi ; k++ ) { if (imp->epClusterHi[k] > 0) // is non zero now { if (was == 0) // if was non zero { was = 1 ; // we start a cluster imp->hicluststart[imp->bighiidx] = k ; // on this one } } else // is zero now { if (was) // if was non zero { was = 0 ; // we end a cluster imp->hiclustend[imp->bighiidx++] = k - 1 ; // on the one before } } } // we should have the hi clusters isolated now. // find their weighted-mean ep's. for ( k = 0 ; k <= imp->bighiidx ; k++ ) { for ( c = imp->hicluststart[k], sum = 0L, prod = 0L ; c <= imp->hiclustend[k] ; c++ ) { UINT n ; sum += (n = imp->epClusterHi[c]) ; prod += n * c ; } imp->Himeanep[k] = (double)prod/(double)sum ; } _is.ImpsExist = 1 ; } // if (VELS && wanttImps && !_is.ImpsExist)...End of implied period processing. ***/ if (wantAlnSp && !_is.ALNSpExists) { SetWindowText(hWndFetch, "Alignment Spectrum") ; // Supply the spectrum of the alignment _is.ALNSpExists = GetSpectrum( pAlFou, j, res, hi_order, lo_order ) ; } if (VELS) { if (wanttImps && _is.ImpsExist && !_is.IMPSpExists) { SetWindowText(hWndFetch,"IMP Spectrum") ; // Supply the spectrum of the implied period _is.IMPSpExists = GetSpectrum( pIMFou, j, res, hi_order, lo_order ) ; } } if (hWndPB1 != NULL) // all done { SendMessage(hWndPB1, PBM_SETPOS, (WPARAM) 0, 0) ; SetWindowText(hWndFetch,"") ; } return 1 ; } int GetSpectrum(LPFOURIER f, int n, double * r, double hi, double lo) { Spectrum( f, n, r, hi, lo, RAW ) ; return 1 ; } double Ancum(double * newa, double * olda, double * acum) { double diff ; * acum += (diff = adiff(* newa, * olda)) ; * olda = * newa ; return diff ; } /***************************************************************** Given a chord of an orbit through the body of that orbit, and through a point (x, y) in the orbital plane, find the location of the chord's other cut of the orbit. This is the location of the "Conjugate Body". *****************************************************************/ /*** void ConjugateBody( int bod, // identifies target body and orbit double x, // ecliptic position of reference point double y, // through which chord also passes double *j // ecliptic position of conjugate body, and ratio ) { double a, b, c, aa, bb, th, u, T ; double aop, lan, inc, e, ta, ra ; double cs, sn ; double xa, ya, za ; double xb, yb, zb ; double xc, yc, zc ; double xt, yt, xu ; double mu, k, ks, g, h, mk, sq ; double fx, tst1, tst2 ; double xv, yv, d1, d2 ; double mT, iT, mC, iC ; // re tangents double xp, yp ; // for pole (join of tangents) LPPELS s ; s = elvar[bod] ; // we are working directly in the orbital plane of the target body // get the elements of the orbit and body for the given epoch a = s->el[MAJ] ; aa = a * a ; e = s->el[ECC] ; c = a * e ; b = sqrt(aa - c*c) ; bb = b * b ; aop = s->el[AOP] ; lan = s->el[LAN] ; inc = s->el[INC] ; ta = s->el[TRA] ; ra = s->el[RAD] ; T = s->pd ; // body's period // refer reference point (x,y) to orbital frame // if target is Moon & Moon orbit, transform heliocentric // ecliptic reference point to geocentric ecliptic if (bod == MOO) { x -= xearth ; y -= yearth ; } cs = cos(lan) ; // pan round OZ by lan sn = sin(lan) ; xa = x * cs + y * sn ; ya = y * cs - x * sn ; za = 0.0 ; cs = cos(inc) ; // tilt around OXa by inc sn = sin(inc) ; xb = xa ; yb = ya * cs + za * sn ; zb = za * cs - ya * sn ; cs = cos(aop) ; // yaw around OZb by aop sn = sin(aop) ; zc = zb ; xc = xb * cs + yb * sn ; xc += c ; // shift origin from orbit focus to orbit centre yc = yb * cs - xb * sn ; // (xc,yc) "orbitocentric" co-ordinates of reference point xt = ra * cos(ta) + c ; // "orbitocentric" co-ordinates of true body yt = ra * sin(ta) ; fx = fabs(xt) ; mu = (yt - yc)/(xt - xc) ; // slope of chord k = yt - mu * xt ; // intercept of chord xu = -k / mu ; // chord crosses major axis here (yu = 0) ks = k * k ; // k^2 mk = mu * k ; // mu * k h = mu*mu + bb/aa ; // mu^2 + (b/a)^2 g = ks - bb ; // k^2 - b^2 sq = sqrt(mk*mk - h * g) ; tst1 = (-mk + sq) / h ; // two xt values, of which one should be that of the true body tst2 = (-mk - sq) / h ; // use the value most different from xt (should cope with any rounding error) xv = (fabs(fx - fabs(tst1)) >= fabs(fx - fabs(tst2))) ? tst1 : tst2 ; yv = b * (u = sqrt(aa - xv*xv)) / a ; // root always evaluates to zero or a positive value // nevertheless, we also have two possible yv values if (-a < xu && xu < a) { // chord crosses major axis inside the orbit if (yt >= 0.0) {yv *= -1.0 ; u *= -1.0 ;} // so yv has opposite sign to yt } else { // chord crosses major axis outside the orbit if (yt < 0.0) {yv *= -1.0 ; u *= -1.0 ;} // so yv has same sign as yt } th = itan(u, xv) ; // eccentric anomaly of conjugate body j[6] = T * (th - e * sin(th)) / twopi ; // time in julian years from perihelion of conjugate body // ratio of distances from reference point to actual and conjugate body d1 = sqrt(pow(xt - xc, 2.0) + pow(yt - yc,2.0)) ; d2 = sqrt(pow(xv - xc, 2.0) + pow(yv - yc,2.0)) ; j[2] = (d2 >= d1) ? d2 / d1 : d1 / d2 ; // keep on the chart // of the tangent to orbit at the true body (is also the line of the true body's velocity) mT = tan(TanOrb(a, e, xt, yt)) ; // slope iT = yt - mT * xt ; // intercept // of the tangent to orbit at the conjugate body mC = tan(TanOrb(a, e, xv, yv)) ; // slope iC = yv - mC * xv ; // intercept // orbitocentric position of intersection of tangents (pole of chord wrt orbit) xp = (iC - iT)/(mT - mC) ; // this works! yp = mT * xp + iT ; // distance from pole to true body j[5] = sqrt(pow(xt - xp, 2.0) + pow(yt - yp,2.0) ) ; // transform (xv, yv) back to heliocentric ecliptic xv -= c ; // return origin to focus cs = cos(-aop) ; // unyaw sn = sin(-aop) ; xb = xv * cs + yv * sn ; yb = yv * cs - xv * sn ; zb = 0.0 ; cs = cos(-inc) ; // untilt sn = sin(-inc) ; xa = xb ; ya = yb * cs + zb * sn ; za = zb * sn - yb * sn ; cs = cos(-lan) ; // unpan sn = sin(-lan) ; j[0] = xa * cs + ya * sn ; j[1] = ya * cs - xa * sn ; if (bod == MOO) { j[0] += xearth ; j[1] += yearth ; } // transform (xp, yp) back to heliocentric ecliptic xp -= c ; // return origin to focus cs = cos(-aop) ; // unyaw sn = sin(-aop) ; xb = xp * cs + yp * sn ; yb = yp * cs - xp * sn ; zb = 0.0 ; cs = cos(-inc) ; // untilt sn = sin(-inc) ; xa = xb ; ya = yb * cs + zb * sn ; za = zb * sn - yb * sn ; cs = cos(-lan) ; // unpan sn = sin(-lan) ; j[3] = xa * cs + ya * sn ; j[4] = ya * cs - xa * sn ; if (bod == MOO) { j[3] += xearth ; j[4] += yearth ; } } void ellip( double a, double b, double xp, double yp, double *j ) // **************************************************************** // Supply co-ordinates of contact points of tangents to an ellipse // with semi-major axis a, and semi-minor axis b, from coplanar, // _external_ point P(xp,yp). If P not external, exit at once, and // return -1.0 in j[4], else return 0.0 in j[4]. // **************************************************************** { double k, s, m, n, p, q ; double xt, aa, bb ; //double xo ; // assign some often-used values aa = a*a ; bb = b*b ; j[4] = 0.0 ; // find value under the square root k = ( n = (s = bb*xp) * xp + aa*yp*yp ) - aa*bb ; // if negative, P is internal and the solutions are complex, so exit if ( k < 0.0 ) { j[4] = -1.0 ; return ; } // else, assign some more values m = aa / n ; p = b / a ; q = yp * sqrt( k ) ; // and compute coordinates of tangent contact points j[0] = xt = m * ( s + q ) ; // 1st xt j[1] = p * sqrt( aa - xt*xt) ; // 1st yt // find intersection xo of would-be tangent with 0X. // If it cuts inside the ellipse, take y as negative. // xo = xp - j[1]*(xp - xt)/(yp - j[1]) ; // if (xo > -a && xo < a ) j[1] *= -1.0 ; j[2] = xt = m * ( s - q ) ; // 2nd xt j[3] = p * sqrt( aa - xt*xt) ; // 2nd yt // find intersection xo of would-be tangent with 0X. // If it cuts inside the ellipse, take y as negative. // xo = xp - j[3]*(xp - xt)/(yp - j[3]) ; // if (xo > -a && xo < a ) j[3] *= -1.0 ; } ***/ //********************* // Tangents to Orbit //********************* int Tangents( int bdy, // if non-zero, use P_ref following, else use input co-ordinates xr, yr int P_ref, // reference body if bdy is non-zero, else zero int P_targ, // body of target orbit double xr, double yr, double zr, double * j // output ) { LPPELS s ; // double h[5] ; double a, b, c, e, lan, aop, inc ; double aa, bb, R, Rsq ; double x, y, z, yh, g, u, v ; double xb, yb, zb ; double xa, ya, za ; double xt1, yt1 ; double xt2, yt2 ; double cs, sn ; // elements of target body // at current epoch s = elvar[P_targ] ; a = s->el[MAJ] ; aa = a * a ; e = s->el[ECC] ; c = a * e ; b = sqrt(aa - c*c) ; bb = b * b ; aop = s->el[AOP] ; lan = s->el[LAN] ; inc = s->el[INC] ; if (bdy) // use heliocentric position of planet P_Ref { // at current epoch as reference point heliotogeo(P_ref) ; xr = x2 ; yr = y2 ; zr = z2 ; } // else use incoming (xr, yr) directly // transform (xr, yr) to target's orbital plane if (P_targ == MOO) // convert to geocentric { xr -= xearth ; yr -= yearth ; } // pan round OZr to long. of target's asc. node cs = cos(lan) ; sn = sin(lan) ; xa = xr * cs + yr * sn ; ya = yr * cs - xr * sn ; za = zr ; // tilt round OXa into the target's orbital plane // with P_ref projected into target's orbital plane // at right angles to ecliptic. // yb = ya / cos(inc) ; // xb = xa ; cs = cos(inc) ; sn = sin(inc) ; yb = ya * cs + za * sn ; zb = za * cs - ya * sn ; xb = xa ; // Aha! This orbit has the FOCUS as its origin, // but the tangent-finder works to its true centre! // So the x-ordinate of ref point is displaced by c. // yaw round OZb to the argument of the target's // perihelion, adding c to the x ordinate to put // the origin at the orbit's true centre. cs = cos(aop) ; sn = sin(aop) ; x = xb * cs + yb * sn ; x += c ; y = yb * cs - xb * sn ; // <----- error was here. I had sn for cs!!! z = zb ; /************ alternative method, using 'ellip()' above, but slower ************** ellip(a, b, x, y, h) ; if (h[4] < 0.0) return 0 ; if ( y < 0.0 && x > -a && x < a ) { h[1] = -h[1] ; h[3] = -h[3] ; } if ( x < -a ) { if ( y >= 0.0 ) { if ( h[0] >= h[2] ) h[3] = -h[3] ; else h[1] = -h[1] ; } else { if( h[0] >= h[2] ) h[1] = -h[1] ; else h[3] = -h[3] ; } } if ( x >= a ) { if ( y >= 0.0 ) { if ( h[0] >= h[2] ) h[1] = -h[1] ; else h[3] = -h[3] ; } else { if ( h[0] >= h[2] ) h[3] = -h[3] ; else h[1] = -h[1] ; } } xt1 = h[0] ; yt1 = h[1] ; xt2 = h[2] ; yt2 = h[3] ; ********/ // skip and return zero if ref. pt. lies inside the target orbit. if (fabs(x) < a && fabs(y) < b * sqrt(aa - x*x) / a ) return 0 ; // compute coordinates of tangent contact points // first in own frame yh = a * y / b ; // y relative to circumambient circle R = sqrt(x*x + yh*yh) ; // R (use the positive root) Rsq = R * R ; // R squared u = aa / R ; v = a * sqrt(1.0 - aa/Rsq) ; // use the positive root g = itan(yh, x) ; // gamma // rotate that frame cs = cos(-g) ; sn = sin(-g) ; xt1 = u * cs + v * sn ; yt1 = b * (v * cs - u * sn) / a ; v *= -1.0 ; xt2 = u * cs + v * sn ; yt2 = b * (v * cs - u * sn) / a ; // transform contact points back to heliocentric ecliptic, // having subtracted c from the x-ordinate to return the // origin to the focus. // contact point one // unyaw xt1 -= c ; cs = cos(-aop) ; sn = sin(-aop) ; xb = xt1 * cs + yt1 * sn ; yb = yt1 * cs - xt1 * sn ; zb = 0.0 ; // untilt cs = cos(-inc) ; sn = sin(-inc) ; ya = yb * cs + zb * sn ; za = zb * cs - yb * sn ; xa = xb ; // ya = yb * cos(inc) ; // xa = xb ; // unpan cs = cos(-lan) ; sn = sin(-lan) ; j[0] = xa * cs + ya * sn ; j[1] = ya * cs - xa * sn ; // contact point two // unyaw xt2 -= c ; cs = cos(-aop) ; sn = sin(-aop) ; xb = xt2 * cs + yt2 * sn ; yb = yt2 * cs - xt2 * sn ; zb = 0.0 ; // untilt cs = cos(-inc) ; sn = sin(-inc) ; ya = yb * cs + zb * sn ; za = zb * cs - yb * sn ; xa = xb ; // ya = yb * cos(inc) ; // xa = xb ; // unpan cs = cos(-lan) ; sn = sin(-lan) ; j[2] = xa * cs + ya * sn ; j[3] = ya * cs - xa * sn ; if (P_targ == MOO) // convert results to ecliptic { j[0] += xearth ; j[1] += yearth ; j[2] += xearth ; j[3] += yearth ; } return 1 ; } //***************************************************************** // Find the Systemics at Epoch on a Random Position in the Ecliptic // --- I/O via struct --- //***************************************************************** double RandSP(RSP * rsp, int C) { int p ; RSP s ; double u, ao, csm, ssm, N ; static double oldep = 0.0 ; s = *rsp ; ep = s.epo ; // get epoch if (ep != oldep) // update elements only if changed epoch update() ; oldep = ep ; for (p = EAR, csm = ssm = N = 0.0 ; p <= SUN ; p++) { if (p == PLU || !BodyList[p]) continue ; heliotogeo(p) ; switch (C) { case GE: u = 2.0 * itan(y_1 - s.yin, x1 - s.xin) ; break ; // geocentric case EQ: u = 2.0 * itan(y_0 - s.yin, x0 - s.xin) ; break ; // equatorial case HE: default: u = 2.0 * itan(y2 - s.yin, x2 - s.xin) ; break ; // heliocentric } csm += cos(u) ; ssm += sin(u) ; N++ ; } s.dir = 0.5 * itan(ssm, csm) ; // 0.0 <= bisectoid direction < pi s.aln = ao = sqrt(csm*csm + ssm*ssm) / (2.0 * N) ; // -0.5 <= bisectoid alignment <= 0.5 s.dsp = arccos(2.0 * ao) ; // 0.0 <= bisectoid dispersion <= pi/2 return s.dir ; } // **************************************************************************** // Get all the heliocentric co-ordinates, directions and radii in the ecliptic. // **************************************************************************** void Heliodirs( LPSYSTEMIC S ) { int p ; for ( p = EAR ; p <= NEP ; // Not the Sun p++ ) { if (p == PLU) continue ; heliotogeo( p ) ; S->x_on_Sun[p] = x2 ; S->y_on_Sun[p] = y2 ; S->z_on_Sun[p] = z2 ; S->hLongs[p] = itan( y2, x2 ) ; // range 0 to twopi S->hLats[p] = atan(z2 / sqrt(x2*x2 + y2*y2)) ; // range -pi/2 to pi/2 S->hRads[p] = sqrt(x2*x2 + y2*y2) ; } }// Heliodirs( LPSYSTEMIC S ) // ****************************************************************************************** // Get all the geocentric co-ordinates, directions and radii on ecliptic or equatorial plane // ****************************************************************************************** void Geodirs(LPSYSTEMIC S, int Centre) { int p ; double x, y, z, L[3], xsum, ysum, zsum, N ; for ( p = MOO ; // not the Earth p <= SUN ; p++ ) { if (p == PLU) continue ; heliotogeo( p ) ; // get body ordinates at ep switch ( Centre ) { case GE: S->gx_on_Earth[p] = x = x1 ; S->gy_on_Earth[p] = y = y_1 ; S->gz_on_Earth[p] = z = z1 ; break ; case EQ: S->ex_on_Earth[p] = x = x0 ; S->ey_on_Earth[p] = y = y_0 ; S->ez_on_Earth[p] = z = z0 ; break ; default: break ; } angles(x, y, z, L) ; // convert to long & lat, or RA & dec. S->gLongs[p] = L[0] ; S->gLats[p] = L[1] ; S->gRads[p] = sqrt(x*x + y*y) ; } for ( p = MOO, N = xsum = ysum = zsum = 0.0 ; p <= SUN ; p++) { double a, b, c ; switch ( Centre ) { case GE: x = S->gx_on_Earth[p] ; y = S->gy_on_Earth[p] ; z = S->gz_on_Earth[p] ; break ; case EQ: x = S->ex_on_Earth[p] ; y = S->ey_on_Earth[p] ; z = S->ez_on_Earth[p] ; break ; default: break ; } if (p != PLU && BodyList[p]) // if this one is on the list { xsum += x ; // sum of ordinates ysum += y ; zsum += z ; N++ ; // number of directions totalled } S->gmean_x = a = xsum / N ; // mean geocentric/equatorial cartesian position of selected bodies S->gmean_y = b = ysum / N ; S->gmean_z = c = zsum / N ; angles(a, b, c, L) ; S->gmean_lo = L[0] ; // long or RA of cartesian mean body S->gmean_la = L[1] ; // lat or dec of cartesian mean body S->gmean_ra = sqrt(a*a + b*b) ; // radius of cartesian mean body in base plane } }// Geodirs(LPSYSTEMIC S, int Centre) // ***************************************************************** // With each of the bodies in turn taken as centre, find the angles // (wrt 1st Aries on Ecliptic) of all non-centre bodies. // These angles are the same for Geocentric and Heliocentric frames, // but different for the Equatorial frame. // ***************************************************************** void Pairdirs(LPSYSTEMIC S, int Centre) { int p, q ; for (p = EAR ; p <= SUN ; p++) // for each body as centre { double xa, ya, xb, yb ; if (p == PLU) continue ; // Get position of centre body { heliotogeo(p) ; if (Centre == EQ) {xa = x0 ; ya = y_0 ;} else {xa = x2 ; ya = y2 ;} } // Get angles to remaining bodies for (q = EAR ; q <= SUN ; q++) { double u, v ; // Exclude body currently at centre if (q == PLU || q == p) continue ; // Get position of non-centric body { heliotogeo(q) ; if (Centre == EQ) {xb = x0 ; yb = y_0 ;} else {xb = x2 ; yb = y2 ;} } u = xb - xa ; v = yb - ya ; S->Pair[p][q] = itan(v, u) ; // the line PQ has the sense P->Q (i.e., P is the vertex) S->PairRad[p][q] = sqrt(u*u + v*v) ; // the radius to Q with P as centre } } }// Pairdirs(LPSYSTEMIC S, int Centre) void RadiusofMeanCircle( LPSYSTEMIC S ) { // needs Heliodirs() int p ; double N, xsum ; for ( // find the radius of the mean circle p = EAR, N = xsum = 0.0 ; p <= SUN ; // allow inclusion of the Sun p++ ) { if (p == PLU || !BodyList[p]) continue ; xsum += sqrt(pow(S->x_on_Sun[p] - S->x_on_Sun[MBD], 2.0) + pow(S->y_on_Sun[p] - S->y_on_Sun[MBD], 2.0)) ; N++ ; } S->rmc = xsum / N ; // radius of mean circle centred on the mean body } /*** void MassSysDironCoM(LPSYSTEMIC S) { // needs Heliodirs() int p ; double M, x, y, u, cs, sn, csum, ssum, cmsum, smsum ; for (p = EAR, M = csum = ssum = cmsum = smsum = 0.0 ; p <= SUN ; p++) { if (p == PLU || !BodyList[p]) continue ; x = S->x_on_Sun[p] - S->x_on_Sun[MAS] ; y = S->y_on_Sun[p] - S->y_on_Sun[MAS] ; u = itan(y, x) ; csum += (cs = cos(2.0 * u)) ; ssum += (sn = sin(2.0 * u)) ; if (p != SUN) { cmsum += masses[p] * cs ; smsum += masses[p] * sn ; } M++ ; } S->bSysDir[MAS] = 0.5 * itan( ssum, csum ) ; // bisectoid on mass centre S->bMSysDir[MAS] = 0.5 * itan( smsum, cmsum ) ; // mass bisectoid on mass centre } ***/ /**** void MassSysDironMB(LPSYSTEMIC S) { // needs Heliodirs() int p ; double M, x, y, u, csum, ssum ; for (p = EAR, M = csum = ssum = 0.0 ; p < SUN ; p++) { if (p == PLU || !BodyList[p]) continue ; x = S->x_on_Sun[p] - S->x_on_Sun[MBD] ; y = S->y_on_Sun[p] - S->y_on_Sun[MBD] ; u = itan(y, x) ; csum += masses[p] * cos(2.0 * u) ; ssum += masses[p] * sin(2.0 * u) ; M++ ; } S->bMSysDir[MBD] = 0.5 * itan( ssum, csum) ; // mass bisectoid on mean body } void GenSysDirandAl(LPSYSTEMIC S) { // needs Pairdirs() int p, q ; static double csum, ssum, u, N ; for ( p = EAR, csum = ssum = N = 0.0 ; p < SUN ; p++ ) { if (p == PLU || !BodyList[p]) continue ; for ( q = p + 1 ; q <= SUN ; q++ ) { double u ; if (q == PLU || !BodyList[q]) continue ; csum += cos(u = 2.0 * S->Pair[p][q]) ; ssum += sin(u) ; N++ ; } } S->bSysDir[GNR] = 0.5 * itan(ssum, csum) ; // Generalised Systemic Direction S->bSysAln[GNR] = sqrt(csum*csum + ssum*ssum) / (2.0 * N) ; // Generalised Systemic Alignment (def II) S->bSysDsp[GNR] = arccos(2.0 * S->bSysAln[GNR]) ; // Generalised Systemic Dispersion } void SysPerIncBody(LPSYSTEMIC S) // position systemics per body { // needs Pairdirs() int p ; for (p = EAR ; p <= SUN ; p++) // all of them { int b ; double a, ao, csm, ssm, N ; if (p == PLU || !BodyList[p]) continue ; { for (b = EAR, csm = ssm = N = 0.0 ; b <= SUN ; b++) { if (b == PLU || b == p || !BodyList[b]) continue ; // exclude body to self csm += cos(a = 2.0 * S->Pair[p][b]) ; ssm += sin(a) ; N++ ; } S->bSysDir[p] = 0.5 * itan(ssm, csm) ; // direction of bisectoid on p S->bSysAln[p] = ao = sqrt(csm*csm + ssm*ssm) / (2.0 * N) ; // alignment of bisectoid on p S->bSysDsp[p] = arccos(2.0 * ao) ; // dispersion of bisectoid on p } } } void MassSysDironBods(LPSYSTEMIC S) // mass bisectoid per body { // needs Pairdirs() int p, b ; static double a, r, mrsq, csm, css, ssm, m, use ; for (p = EAR ; p <= SUN ; p++) // body of origin (not Pluto, the Sun can be an origin ) { if (p == PLU || !BodyList[p]) continue ; // included target masses (not Pluto, not self and the Sun is not a target) // reset cumulators for (b = EAR, csm = css = ssm = 0.0 ; b < SUN ; b++) // target body { if (b == PLU || b == p || !BodyList[b]) continue ; m = masses[b] ; if (wantInertial) { r = S->PairRad[p][b] ; use = m * r * r ; } else use = m ; csm += use * cos(a = 2.0 * S->Pair[p][b]) ; ssm += use * sin(a) ; } S->bMSysDir[p] = 0.5 * itan(ssm, csm) ; // mass bisectoid on p } } void MassSysDironBodsdCos(LPSYSTEMIC S) // mass bisectoids by direction-cosine, per body { // needs Heliodirs() double xsm, ysm, zsm, msm, thx, thy, thz ; double A, B, m, xu, yu, zu, cu, sx, sy, sz ; double dx, dy, dz, mu ; double cx, cy, cz ; int p, b ; for (p = EAR ; p <= SUN ; p++) // Sun can be an origin { if (p == PLU || !BodyList[p]) continue ; // avoid Pluto // all other included bodies are targets for (b = EAR, xsm = ysm = zsm = msm = 0.0 ; b < SUN ; b++) // Sun not a target { if (b == PLU || b == p || !BodyList[b]) continue ; // Self not a target dx = S->x_on_Sun[b] - S->x_on_Sun[p] ; dy = S->y_on_Sun[b] - S->y_on_Sun[p] ; dz = S->z_on_Sun[b] - S->z_on_Sun[p] ; A = itan(dy, dx) ; B = itan(dz, dx / cos(A)) ; xu = cos(A) ; yu = sin(A) ; zu = sin(B) ; cu = cos(B) ; mu = itan( zu, yu ) ; thx = itan( (yu/cos(mu)), xu ) ; // angle to x-axis cx = cos(thx) ; mu = itan( zu, xu ) ; thy = itan( (xu/cos(mu)), yu ) ; // angle to y-axis cy = cos(thy) ; thz = itan( cu, zu ) ; // angle to z-axis cz = cos(thz) ; m = masses[b] ; // target b's relative masS-> msm += m ; // total up included target masseS-> // involving squares, these are always zero or positive, // so sign information is already lost at this stage. xsm += m * cx * cx ; ysm += m * cy * cy ; zsm += m * cz * cz ; } // because xsm/msm + ysm/msm + zsm/msm = 1.0, we get the _systemic_ direction cosines-- sx = sqrt(xsm / msm) ; // sqrt(SUM(mi*cosquared(thxi))/SUM(mi)) = cos(thx0) sy = sqrt(ysm / msm) ; // sqrt(SUM(mi*cosquared(thyi))/SUM(mi)) = cos(thy0) sz = sqrt(zsm / msm) ; // sqrt(SUM(mi*cosquared(thzi))/SUM(mi)) = cos(thz0) //if (S->bMSysDir[p] > pi/2.0) sx *= -1.0 ; // hack!! something better needed // these are 'vector' components in their respective axes, so-- S->malpha[p] = mu = itan(sy, sx) ; // longitude of systemic mass direction by dircos on body p S->mbeta[p] = itan(sz, sx / cos(mu)) ; // latitude of systemic mass direction by dircos on body p } } ***/ void SysAngVel(LPSYSTEMIC S, int Centre) { int p ; double csum, ssum, N ; for ( p = MOO, csum = ssum = N = 0.0 ; p < SUN ; p++ ) { double a, u ; if ( p == PLU || !BodyList[p] ) continue ; csum += cos(2.0 * (S->av[p] = a = AngVel(p, Centre, &u))) ; ssum += sin(2.0 * a) ; S->frq[p] = u ; // frequency in cycles per day N++ ; } S->sysav = 0.5 * itan(ssum, csum) ; S->malsysav = sqrt(csum*csum + ssum*ssum) / (2.0 * N) ; // defn II } void GeoSys(LPSYSTEMIC S) { // needs Geodirs() int p ; double ph, ss[17], cc[17] ; static double u, csum, ssum, N, ao, oep ; for ( p = MOO, csum = ssum = N = 0.0 ; // no Earth p <= SUN ; p++ ) { if (p == PLU || !BodyList[p]) continue ; csum += (cc[p] = cos(u = 2.0 * S->gLongs[p])) ; ssum += (ss[p] = sin(u)) ; N++ ; } ph = S->bSysDir[EAR] = 0.5 * itan( ssum, csum ) ; // Geocentric Systemic Direction (Systemic() returns this value) ao = S->bSysAln[EAR] = sqrt(csum*csum + ssum*ssum) / (2.0 * N) ; // Geocentric Systemic Alignment (defn.II). if (wantShifts) // skip otherwise { for ( p = MOO ; p <= SUN ; p++ ) { double a ; if (p == PLU || !BodyList[p]) continue ; a = 0.5 * itan(ssum - ss[p], csum - cc[p]) ; S->gshift[p] = adiff(a, ph) ; } for ( p = MOO, csum = ssum = 0.0 ; p <= SUN ; p++ ) { double a ; if (p == PLU || !BodyList[p]) continue ; csum += cos(a = 2.0 * S->gshift[p]) ; ssum += sin(a) ; } S->meangshift = 0.5 * itan( ssum, csum ) ; } // Individual geocentric system to body alignments (defn. II) for ( p = MOO ; p <= SUN ; p++ ) { if (p == PLU || !BodyList[p]) continue ; S->bSysAln[p] = ao * cos(2.0 * (ph - S->gLongs[p])) ; } // Get the moon's alignment to system in any event if (!BodyList[MOO]) S->bSysAln[MOO] = ao * cos(2.0 * (ph - S->gLongs[MOO])) ; oep = ep ; } // end of geocentric systemics void HelioSys( LPSYSTEMIC S ) { // needs Heliodirs() int p ; double M, ss[13], cc[13], csum, ssum, sd, sa ; for ( p = EAR, // do include the Earth csum = ssum = M = 0 ; p <= NEP ; // don't include the Sun p++ ) { double u ; if ( p == PLU || !BodyList[p] ) continue ; csum += (cc[p] = cos(u = 2.0 * S->hLongs[p])) ; ssum += (ss[p] = sin(u)) ; M++ ; } sd = S->bSysDir[SUN] = 0.5 * itan(ssum, csum) ; // Systemic Direction on the Sun sa = S->bSysAln[SUN] = 0.5 * sqrt(ssum*ssum + csum*csum) / M ;// Systemic Alignment on the Sun //S->bSysDsp[SUN] = arccos(2.0 * sa) ; // Systemic Dispersion on the Sun if (wantShifts) // "shifts" due to bodies, skip otherwise { for ( p = EAR ; p <= NEP ; // don't include the Sun p++ ) { double a ; if (p == PLU || !BodyList[p]) continue ; a = 0.5 * itan(ssum - ss[p], csum - cc[p]) ; S->hshift[p] = adiff(a, sd) ; } // mean heliocentric shift for ( p = EAR, csum = ssum = 0.0 ; p <= NEP ; // don't include the Sun p++ ) { double a ; if (p == PLU || !BodyList[p]) continue ; csum += cos( a = 2.0 * S->hshift[p] ) ; ssum += sin( a ) ; } S->meanhshift = 0.5 * itan(ssum, csum) ; } // if (wantShifts) } // end Heliocentric Systemics void Heliosystem(LPSYSTEMIC S) // Helioc MB and CoM, with per-body Vels and Accs, and means { // and corresponding coordinates, directions and radii int p ; double a, b, c, d, N, M, Mt ; static double u, v ; double xsum, ysum, zsum ; // pure position double xvsum, yvsum, zvsum ; // velocity double xasum, yasum, zasum ; // acceleration double xmsum, ymsum, zmsum ; // mass position double xmvsum, ymvsum, zmvsum ; // velocity double xmfsum, ymfsum, zmfsum ; // force for ( // find the heliocentric mean position of the selected bodies (The Mean Body) p = EAR, N = M = Mt = xsum = ysum = zsum = xvsum = yvsum = zvsum = xasum = yasum = zasum = xmvsum = ymvsum = zmvsum = xmfsum = ymfsum = zmfsum = xmsum = ymsum = zmsum = 0.0 ; p <= SUN ; // allow inclusion of the Sun p++ ) { if (p == PLU || !BodyList[p]) continue ; heliotogeo(p) ; xsum += x2 ; // position ordinates and sums ysum += y2 ; zsum += z2 ; //S->radsq[p] = x2*x2 + y2*y2 ; if (p != SUN) { // find sums of mass-weighted component positions xmsum += masses[p] * x2 ; ymsum += masses[p] * y2 ; zmsum += masses[p] * z2 ; Mt += masses[p] ; // total mass M++ ; } if (VELS && p != SUN) { // find sums of component velocities xvsum += (S->xvel[p] = xv2) ; yvsum += (S->yvel[p] = yv2) ; zvsum += (S->zvel[p] = zv2) ; // find sums of mass-weighted component velocities //(these are momenta, in fact, not velocities, since momentum = mv, // and become velocities only when divided by Mt.) xmvsum += (S->xmvel[p] = masses[p] * xv2) ; ymvsum += (S->ymvel[p] = masses[p] * yv2) ; zmvsum += (S->zmvel[p] = masses[p] * zv2) ; // find sums of component accelerations xasum += (S->xacc[p] = xa2) ; yasum += (S->yacc[p] = ya2) ; zasum += (S->zacc[p] = za2) ; S->BodAcc[p] = sqrt(xa2*xa2 + ya2*ya2 + za2*za2) ; // find sums of mass-weighted component accelerations //(these are forces, in fact, since f = ma, // and become accelerations only when divided by Mt.) xmfsum += (S->xmf[p] = masses[p] * xa2) ; ymfsum += (S->ymf[p] = masses[p] * ya2) ; zmfsum += (S->zmf[p] = masses[p] * za2) ; } N++ ; } // heliocentric ordinates of the Mean Body (means of ordinates) S->x_on_Sun[MBD] = a = xsum / N ; S->y_on_Sun[MBD] = b = ysum / N ; S->z_on_Sun[MBD] = c = zsum / N ; S->hmean_ra = d = sqrt(a*a + b*b) ; // heliocentric radius of MB S->hmean_lo = itan(b, a) ; // heliocentric direction of MB S->hmean_la = atan(c / d) ; // heliocentric latitude of MB // mass-weighted mean body MWMB (a.k.a. the Mass Centre, or Mass Centroid!) // coordinates of the mass centre S->x_on_Sun[MAS] = a = xmsum / Mt ; S->y_on_Sun[MAS] = b = ymsum / Mt ; S->z_on_Sun[MAS] = c = zmsum / Mt ; S->hMmean_ra = d = sqrt(a*a + b*b) ; // heliocentric radius of Mass Centre S->hMmean_lo = itan(b, a) ; // heliocentric direction of Mass Centre S->hMmean_la = atan(c / d) ; // heliocentric latitude of Mass Centre if (VELS) // find mean component velocities and accelerations { double alpha, beta ; // Re Mean Body (Centroid of Position, CoP) **** // find vector sum and direction of mean // component velocities in ecliptic. // they ought to be those of the Mean Body. S->xvel[MBD] = a = xvsum / N ; S->yvel[MBD] = b = yvsum / N ; S->zvel[MBD] = c = zvsum / N ; // angles of velocity vector of the MB alpha = S->mVectDir = itan(b, a) ; // Direction of Mean radius in the ecliptic on Mean Body S->mVectSum = sqrt(a*a + b*b) ; // Mean radius in the ecliptic on Mean Body beta = S->mVectLat = itan(c, (a / cos(alpha))) ; // Latitude of Mean radius in the ecliptic on Mean Body // Acceleration of Center of Position S->xacc[MBD] = a = xasum / N ; S->yacc[MBD] = b = yasum / N ; S->zacc[MBD] = c = zasum / N ; // angles of acceleration vector of CoP S->maVectDir = itan( b, a ) ; // Direction of acceleration of the Mean Body S->maVectSum = sqrt(a*a + b*b + c*c) ; // Acceleration of the Mean Body // Re Centre of Mass **** // find vector sum and direction of mean mass-weighted // component velocities (sum_of_momenta/Mt) in ecliptic. // they ought to be those of the centre of mass. S->xmvel[MAS] = a = xmvsum / Mt ; S->ymvel[MAS] = b = ymvsum / Mt ; S->zmvel[MAS] = c = zmvsum / Mt ; // angles of velocity vector of the CoM alpha = S->mMVectDir = itan(b, a) ; S->mMVectLat = itan(c, (a / cos(alpha))) ; S->mMVectSum = sqrt(a*a + b*b) ; // in the ecliptic // find vector sum and direction of mean mass-weighted // component accelerations in ecliptic (sum_of_forces/Mt). // they ought to be those of the centre of mass. S->xmacc[MAS] = a = xmfsum / Mt ; S->ymacc[MAS] = b = ymfsum / Mt ; S->zmacc[MAS] = c = zmfsum / Mt ; // angles of acceleration vector of the Mass Centre alpha = S->mMaVectDir = itan(b, a) ; // longitude of the acceleration S->mMaVectLat = itan(c, (a / cos(alpha))) ; // latitude of the acceleration S->mMaVectSum = sqrt(a*a + b*b) ; // magnitude in the ecliptic } } // Heliosystem() /******** double TPofBisectoids(LPSYSTEMIC S) { int p, k ; double cs, sn ; // Frame 1 (heliocentric, ecliptic, x1 along 1st Aries) static double oep ; // previous epoch static double opm[20] ; // previous mass systemic line angle to Aries static double opu[20] ; // previous position systemic line angle to Aries static double omm[20] ; // previous slope of normal to mass systemic line static double omu[20] ; // previous slope of normal to position systemic line static double ocm[20] ; // previous intercept of normal to mass systemic line static double ocu[20] ; // previous intercept of normal to position systemic line double dt ; // time interval between previous and current double vs ; // velocity of line's "origin", S double ps ; // angle to Aries of velocity, vs double pm ; // angle to Aries of systemic line, m double mo ; // slope of current systemic line double m ; // slope of normal to current systemic line double c ; // intercept of normal to current systemic line double w ; // estimated angular velocity of systemic line around S double W ; // angular velocity of systemic line around TP double dx, dy ; double th ; double xs1, ys1 ; // position of S double xk1, yk1 ; // position of point K unit distance from S on systemic line m double xV1, yV1 ; // position of tip of S's velocity vector double x, y ; // position of turning point T in m double xr, yr ; // position of centre of curvature double rho ; // radius of curvature // Frame 2 (origin shifted to K, ecliptic, rotated by pm, x2 along m) double xs2 ; // x ordinate of S double yV2 ; // y ordinate of vs vector tip double yr2 ; // y ordinate of K end of "turning line" double xt2 ; // x ordinate of T, the systemic line's turning point dt = ep - oep ; // time interval between current and previous epochs for (p = EAR ; p < SUN ; p++) // not Sun. (Really? Yes. This is heliocentric. It's all wrt the Sun "at rest") { if (p == PLU || !BodyList[p]) continue ; // just do all included bodies heliotogeo(p) ; // use heliotogeo() to get velocities VELS is on xs1 = x2 ; // current heliocentric position of S ys1 = y2 ; ps = itan(yv2, xv2) ; // longitude of vs vs = xv2 / cos(ps) ; // velocity vs of S (keeping the sign) !!! in a.u. per yr !!! vs /= julianyear ; // now in a.u. per day // two passes, the first using the mass-systemic line, the second the position-systemic line for (k = 0 ; k < 2 ; k++) { if (!k) { pm = S->bMSysDir[p] ; // use mass bisectoid mo = tan(pm) ; // slope of mass bisectoid m = -1.0 / mo ; // slope of normal to mass bisectoid w = adiff(pm, opm[p]) / dt ; // approximate angular velocity of mass systemic line around S } else { pm = S->bSysDir[p] ; // use position bisectoid mo = tan(pm) ; // slope of position bisectoid m = -1.0 / mo ; // slope of normal to position bisectoid w = adiff(pm, opu[p]) / dt ; // approximate angular velocity of systemic line around S } cs = cos(pm) ; sn = sin(pm) ; // position of K in Frame 1 xk1 = xs1 + cs ; yk1 = ys1 + sn ; // position in Frame 1 of vs vector tip xV1 = xs1 + vs * cos(ps) ; yV1 = ys1 + vs * sin(ps) ; // y ordinate in Frame 2 of vs vector tip yV2 = (yV1 - yk1) * cs - (xV1 - xk1) * sn ; // x ordinate of S in Frame 2, // and of the S end of "turning line" in Frame 2 (we have ys2 = yV2) xs2 = (xs1 - xk1) * cs + (ys1 - yk1) * sn ; // y ordinate of K end of "turning line" in Frame 2 (we have xr2 = 0) yr2 = w + yV2 ; // x ordinate of T, the turning point in Frame 2 (we have yt2 = 0) xt2 = xs2 - yV2 * xs2 / (yV2 - yr2) ; // angular velocity of bisectoid around TP W = yV2 / xt2 ; // position of T in Frame 1 (anti-rotate and anti-shift, with yt2 = 0) x = xk1 + xt2 * cos(-pm) ; y = yk1 - xt2 * sin(-pm) ; if (!k) { // re Mass S->TPx[p] = x ; // x of turning point in ecliptic of mass bisectoid S->TPy[p] = y ; // y of turning point in ecliptic of mass bisectoid S->bMSysW[p] = W ; // angular velocity around the turning point c = y - m * x ; // intercept of normal to mass bisectoid xr = -(c - ocm[p]) / (m - omm[p]) ; // x of centre of curvature yr = m * xr + c ; // y of centre of curvature dx = xr - x ; dy = yr - y ; switch(Quadrant(itan(dy, dx))) { case 0: case 3: if (W >= 0.0) th = pm ; else th = pm + pi ; break ; case 1: case 2: if (W >= 0.0) th = pm + pi ; else th = pm ; break ; } S->bMSysA[p] = th ; rho = sqrt(dx*dx + dy*dy) ; // radius of curvature S->bMSysVTP[p] = rho * W ; // "rolling" velocity ocm[p] = c ; omm[p] = m ; opm[p] = pm ; } else { // re Position S->TPBx[p] = x ; // x of turning point in ecliptic of position bisectoid S->TPBy[p] = y ; // y of turning point in ecliptic of position bisectoid S->bSysW[p] = W ; // angular velocity around the turning point c = y - m * x ; // intercept of normal to position bisectoid xr = -(c - ocu[p]) / (m - omu[p]) ; // x of centre of curvature yr = m * xr + c ; // y of centre of curvature dx = xr - x ; dy = yr - y ; switch(Quadrant(itan(dy, dx))) { case 0: case 3: if (W >= 0.0) th = pm ; else th = pm + pi ; break ; case 1: case 2: if (W >= 0.0) th = pm + pi ; else th = pm ; break ; } S->bSysA[p] = th ; rho = sqrt(dx*dx + dy*dy) ; // radius of curvature S->bSysVTP[p] = rho * W ; // "rolling" velocity ocu[p] = c ; omu[p] = m ; opu[p] = pm ; } } // for (k } // for (p oep = ep ; // current ep becomes old ep return th ; } // End turning points of bisectoids ****/ void ModAlign(LPSYSTEMIC S) { Heliodirs(S) ; Geodirs(S, UseCentre) ; Heliosystem(S) ; // SysDirectiononMB(S, UseCentre) ; // MassSysDironMB(S) ; // MassSysDironCoM(S) ; } // *************************************************************** // Systemic() // Obtain // the geocentric Longitudes, or RA's or Azimuths // of all bodies, both cumulative and non-cumulative // then, for the system of "included" bodies, // the Geocentric Systemic Direction // the Generalised Systemic Direction // the Geocentric Systemic Alignment (definition II) // the Generalised Systemic Alignment (definition II) // the Geocentric Alignment (defn.II) of the System to each Body // all at epoch ep. // Added at various times: // Direction of Perigee // Direction of Lunar Ascending Node in Ecliptic // Direction of Lunar Ascending Node in Equatorial Plane // Angle between Ecliptic and Equatorial Planes // Geocentric Systemic Motion Alignment, and its Direction // Generalised Systemic Alignment by all body-pairs, and its Direction // Ordinal sequence of Geocentric Body Positions wrt 1st of Aries // // 19 May 2002 // Storage, input and output by struct // // 2nd Oct 2002 // Angular velocities, systemic angular velocity, // Alignments of angular velocities to systemic angular velocity, // Mean alignment of angular velocities to systemic angular velocity. // // *************************************************************** void Systemic( LPSYSTEMIC S, // input & output struct int Centre // frame ) { int NonBod ; // using "Psuedobodies". int p, HIB ; double N, dsum, xsum, ysum, zsum ; static double sum, csum, ssum, csm, ssm ; Heliodirs(S) ; // Get all the heliocentric co-ordinates, directions and radii in the ecliptic. Geodirs(S, Centre) ; // Get all the geocentric co-ordinates, longitudes and latitudes // or RAs and Declinations appropriate to centre. NonBod = (PlotNEcLo || PlotNEqLo || PlotPeriLo || PlotMot)? 1 : 0 ; HIB = NonBod ? MOT : SUN ; // stop at SUN unless doing a "psuedo" body for ( p = MOO, // exclude Earth dsum = xsum = ysum = zsum = N = 0.0 ; p <= HIB ; // "high body" p++ ) { if ( p == EQN ) { EqNodal() ; // long (or RA) and lat (or dec) of the lunar nodal in the equatorial plane switch( Centre ) { case GE: // as registered wrt the ecliptic and 1st Aries S->gLongs[p] = itan(y_1, x1) ; //Long S->gLats[p] = atan(z1 / sqrt(x1*x1 + y_1*y_1)) ; //Lat if (x1 < 0.0) S->gLats[p] *= -1.0 ; break ; case EQ: // as registered wrt the equatorial plane and 1st Aries S->gLongs[p] = itan(y_0, x0) ; //RA S->gLats[p] = 0.0 ; //Dec break ; } } else if ( p == LEQ ) // angle between equatorial and lunar planes S->gLongs[p] = EqNodal() ; // is indifferent to Centre else if ( p == PER ) // longitude or RA of the lunar perigee { lperi(MOO) ; switch( Centre ) { case GE: // as registered wrt ecliptic S->gLongs[p] = itan(y_1, x1) ; //Long S->gLats[p] = atan(z1 / sqrt(x1*x1 + y_1*y_1)) ; //Lat break ; case EQ: // as registered wrt equatorial plane S->gLongs[p] = itan(y_0, x0) ; // RA S->gLats[p] = atan(z0 / sqrt(x0*x0 + y_0*y_0)) ; // Dec break ; } } else if ( VELS && p == MOT ) // systemic motions { double a ; int cent ; cent = (wantHeliocMots) ? HE : Centre ; S->gLongs[MOT] = SMotion(&a, cent) ; // Direction of Systemic Motion S->bSysAln[MOT] = a ; // Alignment of Motion Direction } else if (p == NOD) // Ascending lunar Node in Ecliptic { heliotogeo(p) ; switch( Centre ) { case GE: // as registered wrt ecliptic S->gLongs[p] = itan(y_1, x1) ; //Long S->gLats[p] = atan(z1 / sqrt(x1*x1 + y_1*y_1)) ; //Lat break ; case EQ: // as registered wrt equatorial plane S->gLongs[p] = itan(y_0, x0) ; // RA S->gLats[p] = atan(z0 / sqrt(x0*x0 + y_0*y_0)) ; // Dec break ; } } } GeoSys(S) ; // Geocentric Systemics (Systemic Direction, Alignment(defn. II) and Dispersion) HelioSys(S) ; // Heliocentric Systemics (Systemic Direction, Alignment(defn. II) and Dispersion) Heliosystem(S) ; // Systemic direction on Cartesian mean body. if (PlotMBhRadius) RadiusofMeanCircle(S) ; // radius of mean circle on MB // SysDirectiononMB(S, Centre) ; // systemic direction on the mean body // MassSysDironCoM(S) ; // systemic MASS direction on the Centre of Mass // MassSysDironMB(S) ; // systemic MASS direction on the MEAN BODY // GenSysDirandAl(S) ; // generalised systemic direction and its alignment // SysPerIncBody(S) ; // systemics per Body //if (PlotMassBscs) // MassSysDironBods(S) ; // systemic MASS direction per Body if (VELS && PlotSysAV) SysAngVel(S, Centre) ; // Systemic angular velocity // --------------------------- // Inter-Body Turning Points // --------------------------- if (VELS && wantInterbodyTPs) { int b1, b2 ; Pairdirs(S, Centre) ; // Get directions of all body-body pairs for (b1 = EAR ; b1 <= SUN ; b1++) { if (b1 == PLU || !ShowInterBodyTP[b1]) continue ; for (b2 = EAR ; b2 <= SUN ; b2++) { if (b2 == PLU || b2 == b1 || !ShowInterBodyTP[b2]) continue ; IBturn(b1, b2, S->Pair[b1][b2], S) ; } } } if (VELS && wantSystemicPrimary) SystemicPrimary(S, MAS) ; } // End of Systemic() //*********************************************************************** // Inter-Body Turning Points (Velocities of bodies should be known) //*********************************************************************** void IBturn(int M1, int M2, double beta, LPSYSTEMIC S) { double cs, sn, cs1, sn1 ; // // in frame 1 // ~~~~~~~~~~ double xm1_1, ym1_1 ; // position of M1 double Vx1_1, Vy1_1 ; // comps of M1 velocity, V1 double xm2_1, ym2_1 ; // position of M2 double Vx2_1, Vy2_1 ; // comps of M2 velocity, V2 // in frame 2 // ~~~~~~~~~~ // with origin of frame 2 shifted to M1(xm1_1, ym1_1) in frame 1 // and frame 2 X-axis in inter-body line m double xm1_2, ym1_2 ; // position of M1 double Vx1_2, Vy1_2 ; // comps of M1 velocity, V1 double xm2_2, ym2_2 ; // position of M2 double Vx2_2, Vy2_2 ; // comps of M2 velocity, V2 double xT_2 ; // position in frame 2 of turning point T (we have yT_2 = 0) beta = range(beta) ; cs = cos(beta) ; cs1 = cos(-beta) ; sn = sin(beta) ; sn1 = sin(-beta) ; xm1_1 = S->x_on_Sun[M1] ; // frame 1 position of M1 ym1_1 = S->y_on_Sun[M1] ; xm2_1 = S->x_on_Sun[M2] ; // frame 1 position of M2 ym2_1 = S->y_on_Sun[M2] ; Vx1_1 = S->xvel[M1] + xm1_1 ; // position of M1 vel vector-tip in frame 1 Vy1_1 = S->yvel[M1] + ym1_1 ; Vx2_1 = S->xvel[M2] + xm2_1 ; // position of M2 vel vector-tip in frame 1 Vy2_1 = S->yvel[M2] + ym2_1 ; // end points of velocity vectors as positions in frame 2 Vx1_2 = (Vx1_1 - xm1_1) * cs + (Vy1_1 - ym1_1) * sn ; // of V1_2 Vy1_2 = (Vy1_1 - ym1_1) * cs - (Vx1_1 - xm1_1) * sn ; Vx2_2 = (Vx2_1 - xm1_1) * cs + (Vy2_1 - ym1_1) * sn ; // of V2_2 Vy2_2 = (Vy2_1 - ym1_1) * cs - (Vx2_1 - xm1_1) * sn ; // positions in frame 2 of M1 and M2 xm1_2 = 0.0 ; ym1_2 = 0.0 ; xm2_2 = (xm2_1 - xm1_1) * cs + (ym2_1 - ym1_1) * sn ; ym2_2 = 0.0 ; // position in frame 2 of the turning point, T (we have yT_2 = 0.0) xT_2 = xm2_2 - Vy2_2 * (xm2_2 - xm1_2)/(Vy2_2 - Vy1_2) ; // position in frame 1 of the turning point, T S->xT[M1][M2] = xm1_1 + xT_2 * cs1 ; S->yT[M1][M2] = ym1_1 - xT_2 * sn1 ; } //***************************************************************** // Position of "Systemic Primary" // Rationale: // The acceleration vector of a body in orbit always points // through the primary body. We have a systemic acceleration // vector that may be considered always to pass through an // "effective" or "systemic" primary. Since primaries are "at // rest", and the radius vectors are turning in them, and since // turning points are always at rest, the effective primary // must be located at the turning point of the line of the // systemic acceleration vector. //***************************************************************** void SystemicPrimary(LPSYSTEMIC S, int type) { double cs, sn ; // Frame 1 (heliocentric, ecliptic, x1 along 1st Aries) static double oep ; // previous epoch static double opm ; // previous acceleration vector angle to Aries double dt ; // time interval between previous and current double vs ; // velocity of line's "origin", S double ps ; // angle to Aries of velocity vector, vs double pm ; // angle to Aries of acceleration vector, (line m) double w ; // estimated angular velocity of acceleration vector around S double W ; // angular velocity of acceleration vector around TP double xs1, ys1 ; // position of S double xk1, yk1 ; // position of point K unit distance from S on acceleration vector m double xV1, yV1 ; // position of tip of S's velocity vector double x, y ; // position of turning point T in m // Frame 2 (origin shifted to K, ecliptic, rotated by pm, x2 along m) double xs2 ; // x ordinate of S double yV2 ; // y ordinate of vs vector tip double yr2 ; // y ordinate of K end of "turning line" double xt2 ; // x ordinate of T, the acceleration vector's turning point dt = ep - oep ; // time interval between current and previous epochs if (type == MBD) { xs1 = S->x_on_Sun[MBD] ; // current heliocentric position of S ys1 = S->y_on_Sun[MBD] ; ps = itan(S->yvel[MBD], S->xvel[MBD]) ; // longitude of vs vs = S->xvel[MBD] / cos(ps) ; // velocity vs of S (preserving the sign) !!! in a.u. per yr !!! vs /= julianyear ; // now in a.u. per day pm = itan(S->yacc[MBD], S->xacc[MBD]) ; // longitude of acceleration vector } if (type == MAS) { xs1 = S->x_on_Sun[MAS] ; // current heliocentric position of S ys1 = S->y_on_Sun[MAS] ; ps = itan(S->ymvel[MAS], S->xmvel[MAS]) ; // longitude of vs vs = S->xmvel[MAS] / cos(ps) ; // velocity vs of S (keeping the sign) !!! in a.u. per yr !!! vs /= julianyear ; // now in a.u. per day pm = itan(S->ymacc[MAS], S->xmacc[MAS]) ; // longitude of acceleration vector } w = adiff(pm, opm) / dt ; // approximate angular velocity of acceleration vector around S opm = pm ; // save current for old cs = cos(pm) ; sn = sin(pm) ; // position of K in Frame 1 xk1 = xs1 + cs ; yk1 = ys1 + sn ; // position in Frame 1 of vs vector tip xV1 = xs1 + vs * cos(ps) ; yV1 = ys1 + vs * sin(ps) ; // y ordinate in Frame 2 of vs vector tip yV2 = (yV1 - yk1) * cs - (xV1 - xk1) * sn ; // x ordinate of S in Frame 2, // and of the S end of "turning line" in Frame 2 (we have ys2 = yV2) xs2 = (xs1 - xk1) * cs + (ys1 - yk1) * sn ; // y ordinate of K end of "turning line" in Frame 2 (we have xr2 = 0) yr2 = w + yV2 ; // x ordinate of T, the turning point in Frame 2 (we have yt2 = 0) xt2 = xs2 - yV2 * xs2 / (yV2 - yr2) ; // angular velocity of acceleration vector around TP W = yV2 / xt2 ; // position of T in Frame 1 (anti-rotate and anti-shift, with yt2 = 0) x = xk1 + xt2 * cos(-pm) ; y = yk1 - xt2 * sin(-pm) ; switch (type) { case MAS: S->SPMx = x ; S->SPMy = y ; break ; case MBD: default: S->SPx = x ; S->SPy = y ; break ; } oep = ep ; }