/**************************************************************************/ /* PROGRAM GEODESIC; */ /* The program "GEODESIC" and procedure "GEO_REDUCTION" are Copyright */ /* 1992 by Gary A. Crowell Sr. Permission is granted for unrestricted */ /* non-commerical use and distribution. */ /* This mainline program demonstrates a very simple application of */ /* the procedure GEO_REDUCTION, for Geodesic tracking data reduction. */ /* It is perfectly useful as it stands, but the user may wish to add */ /* some ruffles and lace, such as printer & file output, the ability to */ /* tag an ID onto the data, and calculation of tracking station height */ /* differences. GEODESIC is written in Borland Turbo Pascal for IBM */ /* flavor computers, but no unique extensions are used so it should be */ /* usable with any computer/compiler combination. */ /* */ /* Translated to C by Warren Massey, 2/22/95 */ /* */ /**************************************************************************/ #include #include #define TRUE 1 #define FALSE 0 #define DegRad 57.29587 /* Constant for converting degrees to radians. */ #define MaxErr 10.0 /* Max allowed closure error percentage */ /*======================================================================*/ /* Begin get_value function. Gets a string from the keyboard, checks */ /* to see if it contains the exit character "X" (or an error), if it */ /* does a TRUE function value is returned immediately, otherwise the */ /* input sting is converted to a floating point number and that number */ /* as well as a FALSE function value are returned. */ /*======================================================================*/ int get_value( value ) float *value; { char input[40]; if( gets( input ) == NULL ) /* if Error or EOF */ { return TRUE; } else { if( toupper( *input ) == 'X' ) /* see if an X for exit was entered */ return TRUE; else { if( sscanf( input, "%f9.0", value ) < 0 ) *value = 0.0; return FALSE; } } } /*==========================================================================*/ /* Begin SQR function. Returns the square of the (float) value passed to it.*/ /*==========================================================================*/ float SQR( r ) float r; { float s; s = r*r; return s; } /*======================================================================*/ /* Begin GEO_REDUCTION subroutine. Geodesic data reduction engine. */ /*======================================================================*/ void GEO_REDUCTION( Alt, PctClose, Baseline, Bheight, Lheight, AAzimuth, AElev, BAzimuth, BElev) float *Alt, *PctClose; /* Results of data reduction. */ float Baseline; /* Baseline length. */ float Bheight; /* Height of tracker B above A. */ float Lheight; /* Height of launcher above A. */ float AAzimuth, AElev; /* Raw tracking data from A. */ float BAzimuth, BElev; /* Raw tracking data from B. */ /***************************************************************************/ /* General purpose routine for geodesic data reduction including */ /* consideration of non-level tracking geometry. This implementation */ /* was developed directly from the derivation of the Geodesic equations */ /* obtained from their author, Mike Ady, in January 1977. */ /* The equations for the special case where the tracking system is */ /* level are included for reference but are commented out in this */ /* listing. If the height differences are zero, the general equations */ /* will reduce to the equivalent of the special level case. */ /* These equations account for height differences in the tracking */ /* system without requiring odd rotations of the coordinate system. To */ /* properly utilize these equations, the trackers >>MUST<< be */ /* independently zeroed to the horizontal. That is, align the azimuth on */ /* the opposite tracker, but align the elevation to the local horizontal */ /* with a bubble level. Aligning the trackers to a non-level plane */ /* passing through the stations and launch site is NOT appropriate. */ /* The result is that there should be no restrictions on height */ /* differences between trackers; one could be at the bottom of a well, */ /* and the other on the roof of a convenient high-rise. (You might not */ /* get many good tracks out of the well, but when you did, the equations */ /* would work fine. */ /* */ /* Tracker A is the origin of the tracking coordinate system. */ /* Tracker B is in the +X axis direction at an x distance equal to the */ /* horizontal length of the baseline. The launch site is assumed to be */ /* displaced from the baseline in the +Y, +X quadrant. Heights above */ /* tracker A are in the +Z direction. */ /* Actually this set of equations is not entirely general, since */ /* they assume that tracker A is at the origin, and that tracker B is */ /* located in the direction of the +X axis. A truly general form would */ /* allow both trackers to be located anywhere in the coordinate system. */ /* A fully general form might be useful in putting together DR for a */ /* three or four tracker system. Any comments or suggestions regarding */ /* this material would be appreciated. */ /* */ /* Gary Crowell */ /* 71461.1525@compuserve.com */ /* garyc@data-cache.com */ /* */ /***************************************************************************/ { float LOSA, LOSB; /* Line-Of-Sight lengths from trackers to */ /* intersection with geodesic line. */ float Sx, Sy, Sz; /* Geodesic vector that connects endpoints */ /* of LOS vectors A and B. */ float Mx, My, Mz; /* Vector from system origin to weighted */ /* midpoint of geodesic line. */ float Ax, Ay, Az; /* Unit vector in direction of LOSA. */ float Bx, By, Bz; /* Unit vector in direction of LOSB. */ float F, InvF; /* Intermediate terms used in calculations. */ /*Conversion of raw tracking data from degrees to radians for Pascal */ /*trig functions. */ AAzimuth = AAzimuth / DegRad; AElev = AElev / DegRad; BAzimuth = BAzimuth / DegRad; BElev = BElev / DegRad; /*Calculate unit vector of line-of-sight from tracker A; zero azimuth */ /*is aligned on tracker B with positive angles swinging towards the */ /*launch site. */ Ax = cos(AElev) * cos(AAzimuth); Ay = cos(AElev) * sin(AAzimuth); Az = sin(AElev); /*Calculate unit vector of line-of-sight from tracker B; zero azimuth */ /*is aligned on tracker A with positive angles swinging towards the */ /*launch site. */ Bx = -cos(BElev) * cos(BAzimuth); By = cos(BElev) * sin(BAzimuth); Bz = sin(BElev); /*Calculation of intermediate terms. The following equation for F */ /*is a simplification of the corresponding equation typically */ /*published for this factor. */ F = (Az * Bz) - cos(AElev) * cos(BElev) * cos(AAzimuth + BAzimuth); InvF = 1.0 / (1.0 - pow(F,2)); /*The following block of equations are the Geodesic method as */ /*presented in the September 83 Model Rocketeer (and corrected in the */ /*11/83 issue), with some minor simplification. They do not account */ /*for height differences between the tracking stations or launch site. */ /*If this block is uncommented, and the other lines below are */ /*commented out, this block will give the same results as the general */ /*equations for a level-case tracking system. */ /* LOSB = (F * Ax - Bx) * InvF; */ /* LOSA = (Ax - F * Bx) * InvF; */ /* Alt = LOSB * LOSA * (Az + Bz) / (LOSB + LOSA); */ /* PctClose = 100.0 * SQRT(InvF) * ABS(COS(BElev) * Az * */ /* SIN(BAzimuth) - COS(AElev) * Bz * SIN(AAzimuth)) / Alt; */ /* Alt = ROUND(Alt * Baseline); */ LOSB = (Baseline * (F * Ax - Bx) + Bheight * (F * Az - Bz)) * InvF; LOSA = (Baseline * (Ax - F * Bx) + Bheight * (Az - F * Bz)) * InvF; /*Note that Mx and My are not required to simply find the altitude. */ /*These terms are useful if you desire the ground position of the */ /*track. */ Mx = (LOSB * LOSA * (Ax + Bx) + LOSA * Baseline) / (LOSA + LOSB); My = (LOSB * LOSA * (Ay + By) ) / (LOSA + LOSB); Mz = (LOSB * LOSA * (Az + Bz) + LOSA * Bheight ) / (LOSA + LOSB); Sx = (LOSB * Bx + Baseline) - LOSA * Ax; Sy = (LOSB * By ) - LOSA * Ay; Sz = (LOSB * Bz + Bheight ) - LOSA * Az; *Alt = Mz - Lheight; /*A minor point of contention might be the appropriate altitude to use */ /*as the normalizing factor (the denominator) in the closure equation. */ /*Altitude above the launcher (used here)? Altitude above the origin? */ /*Something else? */ /*Note also that this equation and the use of the 'S' vector do not */ /*appear similar to the 'standard' equations. They are identical, and */ /*this form makes it easier to see that we are talking about the */ /*length of the geodesic line. */ *PctClose = fabs(100.0 * sqrt(SQR(Sx) + SQR(Sy) + SQR(Sz)) / *Alt); return; /*PROCEDURE GEO_REDUCTION*/ } /*======================================================================*/ /* Begin GEODESIC MAINLINE */ /*======================================================================*/ main() { float BASE, EA, EE, WA, WE, LH, WH, A, C; printf("\n\n\r"); printf("Program \"GEODESIC\" and procedure \"GEO_REDUCTION\" are\n\r"); printf("Copyright 1992 by Gary A. Crowell, Sr., [71461,1525].\n\r"); printf("(\"C\" version of the program by Warren Massey.)\n\r"); printf("This program is not yet approved for use by the NAR.\n\n\n\r"); printf("Enter X to Exit\n\n\r"); printf("Baseline Length? \n\r"); printf(" (horizontal distance between trackers) \n\r"); printf( " (units are arbitrary and altitudes will be in the same units)\n\r"); printf(">>> "); if( get_value( &BASE ) ) return; printf("\n\n\r"); printf("Height difference between Tracker #1 and Launch Site?\n\r"); printf(" (a positive value if Launcher is higher than Tracker #1) \n\r"); printf(" (units are the same as those used for the baseline)\n\r"); printf(">>> "); if( get_value( &LH ) ) return; printf("\n\n\r"); printf("Height difference between Tracker #1 and Tracker #2?\n\r"); printf(" (a positive value if Tracker #2 is higher than Tracker #1) \n\r"); printf(" (units are the same as those used for the baseline)\n\r"); printf(">>> "); if( get_value( &WH ) ) return; printf("\n\n\r"); printf(" Baseline length = %9.2f \n\r", BASE); printf(" Launcher height relative to Tracker #1 = %9.2f ", LH ); if( LH > 0.0 ) printf("(Launcher is higher)\n\r"); else { if( LH < 0.0 ) printf("(Launcher is lower)\n\r"); else printf("(Same elevation)\n\r"); } printf("Tracker #2 height relative to Tracker #1 = %9.2f ", WH ); if( WH > 0.0 ) printf("(Tracker #2 is higher)\n\r"); else { if( WH < 0.0 ) printf("(Tracker #2 is lower)\n\r"); else printf("(Same elevation)\n\r"); } for(;;) /*until X is pressed to exit*/ { printf("\n\n\r"); printf("Tracker #1 Azimuth? "); if( get_value( &EA ) ) break; printf("Tracker #1 Elevation? "); if( get_value( &EE ) ) break; printf("Tracker #2 Azimuth? "); if( get_value( &WA ) ) break; printf("Tracker #2 Elevation? "); if( get_value( &WE ) ) break; GEO_REDUCTION( &A, &C, BASE, WH, LH, EA, EE, WA, WE ); printf("\n\n\r"); printf(" Tracker #1 Azimuth = %.2f Tracker #2 Azimuth = %.2f \n\r", EA, WA ); printf(" Tracker #1 Elevation = %.2f Tracker #2 Elevation = %.2f \n\r", EE, WE ); /*The final altitude should be rounded off to whole meters here. */ printf("\n Altitude = %.0f Closure = %.2f%% \n\r", A, C ); if( C > MaxErr ) printf("\n\007 >>> TRACK DID NOT CLOSE <<<\n\n\r"); } return; /* GEODESIC MAINLINE */ } /*======================================================================*/ /*SAMPLE TEST DATA: */ /*BASELINE = 300; LAUNCHER HEIGHT = 0; TRACKER WEST HEIGHT = 0 */ /* EAST AZ EAST EL WEST AZ WEST EL ALT CLOSE */ /* ------- ------- ------- ------- ------ ----- */ /* 90 45 50 40 380.0 6.3 */ /* 30 45 60 45 203.3 31.2 */ /* 120 75 25 55 596.3 6.2 */ /* 30 80 40 85 1337.8 3.2 */ /*BASELINE = 300; LAUNCHER HEIGHT = 10; TRACKER WEST HEIGHT = 0 */ /* EAST AZ EAST EL WEST AZ WEST EL ALT CLOSE */ /* ------- ------- ------- ------- ------ ----- */ /* 90 45 50 40 370.0 6.5 */ /* 30 45 60 45 193.3 32.8 */ /* 120 75 25 55 586.3 6.3 */ /* 30 80 40 85 1327.8 3.2 */ /*BASELINE = 300; LAUNCHER HEIGHT = 0; TRACKER WEST HEIGHT = 10 */ /* EAST AZ EAST EL WEST AZ WEST EL ALT CLOSE */ /* ------- ------- ------- ------- ------ ----- */ /* 90 45 50 40 386.7 8.0 */ /* 30 45 60 45 208.1 27.7 */ /* 120 75 25 55 609.1 5.8 */ /* 30 80 40 85 1340.4 3.3 */ /*BASELINE = 300; LAUNCHER HEIGHT = 10; TRACKER WEST HEIGHT = 10 */ /* EAST AZ EAST EL WEST AZ WEST EL ALT CLOSE */ /* ------- ------- ------- ------- ------ ----- */ /* 90 45 50 40 376.7 8.3 */ /* 30 45 60 45 198.1 29.1 */ /* 120 75 25 55 599.1 5.8 */ /* 30 80 40 85 1330.4 3.3 */