******************************************************************************** * * * ************************* * * * * * * * 'B A T C H' - PROGRAM * * * * * * * ************************* * * * * PARAMETRIC LEAST SQUARES ADJUSTMENT * * ----------------------------------- * * OF DOPPLER SATELLITE OBSERVATIONS * * --------------------------------- * * *RANGE DIFFERENCE MODE**ORBIT KNOWN CASE* * * ----------------------------------------- * * *SINGLE TRACKING STATION**MULTI-SATELLITE PASSES* * * ------------------------------------------------- * * IN A BATCH ITERATIVE SOLUTION * * ----------------------------- * * * * * * PREPARED BY : * * 'MOHAMED M. NASSAR' / *G.S.* - (DECEMBER 1971) * * ----------------- * * * * THE PURPOSE OF THIS PROGRAM IS TO PERFORM THE BATCH ITERATIVE SOLUTION FOR * * THE PARAMETRIC LEAST SQUARES ADJUSTMENT OF THE DOPPLER SATELLITE OBSERVATION* * TAKEN FROM ONE TRACKING STATION ,UTILIZING THE 'ORBIT-KNOWN CASE'. * * * * THE OBSERVED QUANTITIES ARE THE TWO-MINUTES INTEGRATED DOPPLER COUNTS,TAKEN * * OVER DIFFERENT PERIODS OF OBSERVING MULTI-SATELLITE PASSES. THE UNKNOWN * * PARAMETERS ARE THE THREE-DIMENSIONAL CARTESIAN (XYZ) COORDINATES OF THE * * TRACKING STATION ,AND THE FREQUENCY OFFSETS (ONE PER EACH OBSERVED SATELLITE* * * * THE MAIN PROGRAM 'BATCH' IS IN (MATLAN) FOR PERFORMING THE L.S. ADJUSTMENT * * ,AND PROVIDING MANY CHECKS TO VERIFY THE CONSISTENCY OF THE COMPUTATIONS. * * * * THE DESIGN MATRICES ,NAMELY ;THE COEFFICIENT MATRIX-A OF THE UNKNOWN PARA- * * METERS AND THE MISCLOSURE VECTOR-W OF THE OBSERVATION EQUATIONS ,ARE EVAL- * * UATED VIA A (FORTRAN)-SUBROUTINE 'DESIGN' BY USING ITS ENTRY-POINT 'AANDW' * * FOR THE SUCCESSIVE ITERATIONS. * * * * THE INVERSION OF THE NORMAL EQNS. MATRIX IS DONE VIA A (MATLAN) SUBPROGRAM -* * 'SNINV' ,WHICH IN TURN CAN CALL ANOTHER (MATLAN) SUBPROGRAM -'INVERT' IF * * SUCCESSIVE PARTITIONING TECHNIQUES ARE DESIRED. * * * * THE ITERATION PROCESS IS STOPPED WHEN THE GIVEN DESIRED PRECISION LIMITS * * ARE SATISFIED ; HOWEVER, IT WILL BE TERMINATED AFTER BEING REPEATED FIVE * * TIMES ,WITH A GIVEN MESSAGE OF (NO CONVERGENCE). * * * * THE FINAL ESTIMATED VALUES OF THE (X,Y & Z) GEOCENTRIC CARTESIAN COORDINATES* * OF THE TRACKING STATION ,WILL BE CONVERTED TO ITS CORRESPONDING (LAT,LONG & * * HT) W.R.T. ANY DESIRED REFERENCE DATUM (EITHER ABSOLUTELY OR RELATIVELY * * ORIENTED) ;VIA THE (FORTRAN)-SUBROUTINE PACKAGE 'FINRES'. * * * * THE EIGEN-VALUE PROBLEM WILL BE DONE ON THE VARIANCE-COVARIANCE MATRIX OF * * THE ADJUSTED COORDINATES TO COMPUTE THE ERROR ELLIPSOID (OR (STANDARD CONFI-* * DENCE REGION) ; VIA THE FORTRAN SUBROUTINE 'EIGEN3'. * * * * THE FINAL ESTIMATED VALUES OF THE FREQUENCY OFFSETS FOR ALL THE USED PASSES * * IN THE L. S. ADJUSTMENT WILL BE TRANSFERED TO THE (FORTRAN)-SUBROUTINE * * 'RECORD' ,TO BE IMPLEMENTED WITH THE OTHER CHARACTERISTICS AND IDENTIFI- * * CATION PARAMETERS OF THESE PASSES ,INTO A BOOK-KEEPING RECORDING TABLE * * * * THE CHARACTERISTICS OF THE INPUT & OUTPUT DATA ARE NICELY PRESENTED THROUGH-* * OUT THE PROGRAM LISTINGS VIA THE EXTENSIVE COMMENT STATEMENTS ,,BUT HERE THE* * INPUT & OUTPUT ARE JUST SUMMARIZED. * * * * THE REQUIRED INPUT DATA ARE IN SEQUENCE AS FOLLOWS : * * -------------------------------------------------- * * THE INPUT DATA CAN BE DIVIDED INTO TWO CATEGORIES : I AND II * * * * I-DATA TO BE PUNCHED BY THE USER ,WHICH ARE : * * 1-HEADING : ANY DESIRABLE HEADING TO APPEAR IN THE OUTPUT UP TO 80 CHARACT.* * 2-THE REQUIRED NUMBER OF COPIES (NCOPIES) TO BE PRINTED-OUT FOR THE * * SUMMARY OF THE FINAL RESULTS AFTER THE L. S. ADJUSTMENT - REQUIRE 1 CARD* * 3-THE PROGRAM OPTIONS (ICPQLA,ICPELA) FOR COMPUTING & PRINTING THE WEIGHT- * * CEFF. AND/OR THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED OBSERVATIONS * * (DOPPLER COUNTS) - REQUIRE 2 CARDS. * * 4-THE OPTIONS & ARGUMENTS (HOWINV,C1,RC1,INC & SCALE) NECESSARY FOR THE * * INVERSION PROCESS - REQUIRE 5 CARDS. * * 5-THE DESIRED PRECISION LIMITS (PRECXYZ,PRECFREQ) NECESSARY TO STOP THE * * ITERATION PROCESS - REQUIRE 2 CARDS. * * * *II-DATA ALREADY AVAILABLE ON PUNCHED CARDS FROM THE 'PRE-BATCH' PROGRAM OUT- * * PUT ,NAMELY : * * 1-THE INFORMATIONS (NPASS,NSUMJK,NAROWS & NACOLS) NECESSARY FOR THE SPACE * * STORAGE ALLOCATION. * * 2-THE YEAR OF OBSERVATIONS (IYEAR) ,AND THE APRIORI VARIANCE FACTOR (AVF). * * 3-THE WEIGHT-COEFFICIENT MATRIX (ELB) OF THE OBSERVED DOPPLER COUNTS. * * 4-THE INITIAL APPROXIMATE VALUES (X0-VECTOR) OF THE UNKNOWN PARAMETERS. * * 5-FOR EACH PASS ; THE NUMBER OF OBSERVED SATELLITE POSITIONS (NUMPOS) ,THE * * COORDINATES OF EACH POSITION (X,Y,Z) ,AND THE CORRESPONDING DOPPLERS * * 6-FOR EACH PASS : THE IDEN. CARD ,MAX. ELEVATION ,TRACK DIRECTION ,AZIMUTH * * AND FINALLY THE QUADRANT. * * 7-THE PARAMETERS OF THE REFERENCE DATUM (SEMI-MAJOR AXIS ,1/FLAT. AND THE * * TRANSLATIONS FROM THE GEOCENTRE) ,WITH THE TRACKING STATION NAME & NO. * * 8-THE PARAMETERS OF ANY ARBITRARILY DESIRED GEOCENTRIC ELLIPSOID (IF ANY). * * * * THE MAIN-OUTPUT WILL BE AS FOLLOWS : ( NOT NECESSARILY IN SEQUENCE ) * * ---------------------------------- * * * * 1-'X'-VECTOR ; CORRECTIONS TO THE APPROXIMATE VALUES OF THE PARAMETERS. * * 2-'XA'-VECTOR ; THE ADJUSTED VALUES OF THE UNKNOWN PARAMETERS. * * 3-'QXA'-MATRIX ; THE WEIGHT COEFFICIENT MATRIX OF THE ADJUSTED PARAMETERS. * * 4-'V'-VECTOR ; THE RESIDUALS OR CORRECTIONS TO THE OBSERVED QUANTITIES. * * 5-'V'PV'-SCALAR ; THE QUADRATIC FORM OF THE WEIGHTED RESIDUALS . * * 6-'EVF'-SCALAR : THE ESTIMATED VARIANCE FACTOR. * * * * 7-'EXA'-MATRIX ; THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED PARAMETERS. * * 8-'QLA'-MATRIX ; THE WEIGHT COEFFICIENT MATRIX OF THE ADJUSTED OBSERVATIONS * * (IF REQUESTED THROUGH THE OPTION -ICPQLA) * * 9-'ELA'-MATRIX : THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED OBSERVATIONS* * (IF REQUESTED THROUGH THE OPTION -ICPELA) * * 10-THE SUMMARY OF THE INPUT AND OUTPUT RESULTS ,AFTER THE SIMULTANEOUS L. S. * * ADJUSTMENT ,WHICH WILL INCLUDE THE FOLLOWING : * * * * * * 1-THE BOOK-KEEPING RECORD OF IDENTIFICATION PARAMETERS AND MAIN CHARACTER- * * ISTICS OF THE USED PASSES WILL BE SHOWN INTO A CAREFULLY ARRANGED TABLE. * * 2-THE FINAL ADJUSTED VALUES OF THE (LAT,LONG & HT) OF THE TRACKING STATION * * W.R.T. THE GIVEN REFERENCE DATUM ,WILL BE SHOWN INTO A FANCY TABLE. * * 3-THE VARIANCE -COVARIANCE MATRIX OF THE ADJUSTED (X,Y,Z) COORDINATES OF THE* * TRACKING STATION ,WITH THE CORRESPONDING STANDARD DEVIATIONS WILL BE * * EXTRACTED AND PRINTED-OUT. * * 4-THE SEMI-AXES OF THE STANDARD ERROR ELLIPSOID (STANDARD CONFIDENCE REG- * * ION ) OF THE ADJUSTED COORDINATES OF THE TRACKING STATION. * * * * NOTE THAT THERE WILL BE SOME OTHER RESULTS WHICH APPEAR IN THE PRINTED * * OUTPUT ,BUT MOST OF THEM ARE CHECKS ON THE PERFORMANCE OF THE COMPUTATIONS * * WITHIN THE PROGRAM SEGMENTS. * * * ******************************************************************************** * MAIN BATCH TEXT 'EXECUTION OF THIS JOB STARTED AT',T ATTRIB (N,NINVERSE),DEF=1 ATTRIB (HOWINV,C1,RC1,INC),INT=1 ATTRIB (NPASS,NAROWS,NACOLS,IYEAR),INT=1 ATTRIB (NCOPIES,ICPQLA,ICPELA),INT=1 * ******************************************************************************** * * *** R E A D AND P R I N T THE I N P U T D A T A *** * HEADING * * READ-IN THE NUMBER OF REQUIRED COPIES TO BE PRINTED-OUT FOR THE SUMMARY * OF THE FINAL RESULTS AFTER THE ITERATIVE LEAST SQUARES ADJUSTMENT. * READ NCOPIES * * READ-IN THE PROGRAM OPTIONS FOR COMPUTING AND PRINTING THE WEIGHT-COEFF. * MATRIX AND/OR THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED OBSERVATIONS * NAMELY ; ICPQLA & ICPELA RESPECTIVELY * FOR ICPQLA=1 : COMPUTE AND PRINT THE WEIGHT-COEFFICIENT MATRIX (QLA) * FOR ICPELA=1 : COMPUTE AND PRINT THE VARIANCE-COVARIANCE MATRIX (ELA) * READ (ICPQLA,ICPELA) * * READ-IN THE INDICATOR 'HOWINV' ,WHICH DEFINES THE CHOSEN WAY FOR EXECUTION * OF THE INVERSION PROCESS ,NAMELY : * * * HOWINV=0 : STRAIGHT-FOREWARD INVERSION OF THE WHOLE MATRIX WITHOUT PARTITIONS * HOWINV<0 : INVERSION BY SUCCESSIVE PARTITIONING USING ONE-INVERSION-STEP * ,HOWINV IN THIS CASE CAN BE ANY -VE NUMBER. * HOWINV>0 : STRAIGHT-FOREWARD INVERSION AND/OR SUCCESSIVE PARTITIONING USING * TWO-INVERSION-STEPS. THE +VE NUMBER ASSIGNED TO HOWINV IN THIS * CASE REPRESENTS THE CHOSEN SIZE OF THE LOWER RIGHT PART OF THE * NORMAL EQNS. MATRIX TO BE INVERTED IN THE SECOND INVERSION-STEP * READ HOWINV * * READ-IN THE ARGUMENTS (C1,RC1,INC) WHICH ARE NECESSARY FOR 'CALLING' THE * SUBPROGRAM 'SNINV' ,FOR PERFORMING THE INVERSION PROCESS OF THE NORMAL EQNS. * MATRIX ; THE VALUES OF THESE ARGUMENTS INDICATE THE FOLLOWING : * C1=0 : STRAIGHT-FOREWARD INVERSION OF THE FIRST-INVERSION-STEP * C1>0 : SUCCESSIVE PARTITIONING OF THE FIRST-INVERSION-STEP ,THE +VE NUMBER * ASSIGNED TO C1 GIVES THE SIZE OF 1ST SUBMATRIX TO BE INVERTED HERE * RC1=0 : STRAIGHT-FOREWARD INVERSION OF THE SECOND-INVERSION-STEP * RC1>0 : SUCCESSIVE PARTITIONS OF THE SECOND-INVERSION-STEP ,THE +VE NUMBER * ASSIGNED TO RC1 GIVES THE SIZE OF 1ST SUBMATRIX TO BE INVERTED HERE * INC>0 : IS THE INCREMENT TO BE ADDED TO (C1 & RC1) A NUMBER OF TIMES FOR * PERFORMING THE SUCCESSIVE PARTITIONING NECESSARY TO GET THE TOTAL * SIZE OF THE MATRIX BEING INVERTED. * READ (C1,RC1,INC) * * READ-IN THE SUITABLE SCALE NEEDED TO BE MULTIPLIED BY THE NORMAL EQNS. * MATRIX BEFORE PERFORMING THE INVERSION PROCESS (IF NECESSARY) * READ SCALE DOUBLE SCALE * * READ-IN THE DESIRED PRECISION LIMITS (PRECXYZ & PRECFREQ) ,AGAIST WHICH THE * SOLUTION VECTOR (X) IS TESTED TO STOP THE ITERATION PROCESS ;WHEN ITS * VALUES ARE LESS THAN OR EQUAL TO THESE GIVEN PRECISION LIMITS * * PRECXYZ : STANDS FOR THE COORDINATES OF THE TRACKING STATION -(IN METERS) * PRECFREQ : STANDS FOR THE FREQUENCY OFFSETS - UNITS ARE IN CYCLES/SECOND. * READ (PRECXYZ,PRECFREQ) DOUBLE PRECXYZ,PRECFREQ FEED 2 TEXT 'PRINTOUT OF THE GIVEN INPUT DATA',C TEXT '================================',C FEED 2 TEXT 'REQUIRED NUMBER OF PRINTED COPIES FROM THE SUMMARY + OF THE FINAL RESULTS = #NCOPIES',5 FEED 2 TEXT 'OPTIONS FOR COMPUTING & PRINTING OF QLA AND/OR ELA + MATRIX ARE :',10 TEXT '---------------------------------------------------+ ----------',10 FEED 1 * * * TEXT 'ICPQLA = #ICPQLA : WHEN EQUALS 1 ,COMPUTE AND PRIN+ T QLA-MATRIX.',5 FEED 1 TEXT 'ICPELA = #ICPELA : WHEN EQUALS 1 ,COMPUTE AND PRIN+ T ELA-MATRIX.',5 FEED 2 TEXT 'OPTIONS AND ARGUMENTS NEEDED FOR THE INVERSION PROC+ ESS ARE :',10 TEXT '---------------------------------------------------+ --------',10 FEED 1 TEXT 'HOWINV = #HOWINV : IS THE INDICATOR WHICH DEFINES + THE CHOSEN WAY TO EXECUTE THE INVERSION PROCESS',5 FEED 1 TEXT 'C1 = #C1 : IS THE SIZE OF THE FIRST SUBMATRIX TO B+ E INVERTED IN THE FIRST-INVERSION-STEP',9 FEED 1 TEXT 'RC1 = #RC1 : IS THE SIZE OF THE FIRST SUBMATRIX TO+ BE INVERTED IN THE SECOND-INVERSION-STEP',8 FEED 1 TEXT 'INC = #INC : IS THE INCREMENT TO BE ADDED TO C1 & + RC1 TO PERFORM THE SUCCESSIVE PARTITIONING TECHNIQUE',8 FEED 1 TEXT 'SCALE : IS THE SUITABLE SCALE TO BE MULTIPLIED + BY THE NORMAL EQNS. MATRIX BEFORE ITS INVERSION',6 TEXT '(THE CHOSEN VALUE OF THIS SCALE IS AS FOLLOWS:)',18 WRITE SCALE,FORMAT=A5 * * NEWPAGE FEED 2 TEXT 'THE DESIRED PRECISION LIMITS TO STOP THE ITERATION + PROCESS ARE :',10 TEXT '---------------------------------------------------+ -----------',10 FEED 1 TEXT 'PRECXYZ : IS THE PRECISION LIMIT FOR THE COORDINAT+ ES OF THE TRACKING STATION -(IN METERS)',5 FEED 1 TEXT 'PRECFREQ : IS THE PRECISION LIMIT FOR THE FREQUENCY+ OFFSETS -(IN CYCLES/SECOND)',5 FEED 1 WRITE (PRECXYZ,PRECFREQ),FORMAT=A5 MULT 600.,PRECFREQ,PRECFREQ * ******************************************************************************** * * READ-IN THE REQUIRED DATA HERE FROM THE 'PRE-BATCH'-FORTRAN PROGRAM NAMELY ; * 1-THE TOTAL NUMBER OF ACCEPTED PASSES (NPASS) ,AND THE ACCUMULATED NUMBER * OF OBSERVED SATELLITE POSITIONS FROM THESE PASSES (NSUMJK) ,WHICH ARE * NEEDED TO ALLOCATE STORAGE SPACE IN THE FORTRAN-SUBROUTINE 'DESIGN' , * FOR THE SATELLITE COORDINATES AND THE CORRESPONDING DOPPLER COUNTS. * * * * 2-THE ROW DIMENSION OF THE DESIGN MATRICES A & W (NAROWS) ,AND THE COLUMN * DIMENSION OF THE A-MATRIX (NACOLS) ,WHICH ARE NEEDED TO ALLOCATE STORAGE * SPACE IN BOTH THE MAIN PROGRAM 'BATCH' AND THE SUBROUTINE 'DESIGN' FOR * THESE MATRICES * 3-THE YEAR OF OBSERVATIONS ,AND THE APRIORI VARIANCE-FACTOR (C/S)-SQUARED * 4-THE WEIGHT COEFFICIENT MATRIX OF THE OBSERVED ACCEPTED DOPPLERS (ELB) * ,AND THE INITIAL APPROXIMATE VALUES OF THE UNKNOWN PARAMETERS (X0). * READ (NPASS,NSUMJK) READ (NAROWS,NACOLS) READ (IYEAR,AVF) READ (ELB,X0) FEED 2 TEXT 'ALLOCATE STORAGE SPACE-DIMENSION INFORMATIONS ARE :+ ',10 TEXT '-------------------------------------------------',+ 10 FEED 1 TEXT 'NPASS = #NPASS :IS THE TOTAL NUMBER OF THE ACCEP+ TED PASSES TO BE PROCESSED HERE',5 FEED 1 TEXT 'NSUMJK = #NSUMJK :IS THE ACCUMULATED NUMBER OF THE+ OBSERVED SATELLITE POSITIONS FROM ALL PASSES',5 * FEED 1 TEXT 'NAROWS = #NAROWS :IS THE TOTAL NUMBER OF ACCEPTED + DOPPLER EQUATIONS (ROW DIMENSION OF A & W MATRICES)',5 FEED 1 TEXT 'NACOLS = #NACOLS :IS THE TOTAL NUMBER OF THE UNKNO+ WN PARAMETERS ( COLUMN DIMENSION OF THE A-MATRIX)',5 FEED 5 TEXT 'Y E A R O F O B S E R V A T I O N S W A S #IYE+ AR',10 NEWPAGE FEED 2 TEXT 'THE APRIORI VARIANCE-FACTOR (AVF) IS = #AVF (CYC+ LES/SEC)-SQUARED',10 FEED 2 TEXT 'ELB-MATRIX : WEIGHT-COEFFICIENT MATRIX OF THE OBSER+ VED QUANTITIES',C WRITE ELB,FORMAT=C5 TEXT 'INVERSION OF (ELB) MATRIX BEGAN AT',T CALL INVERT,(30,150,30,ELB,P) TEXT 'INVERSION OF (ELB) MATRIX ENDED AT',T CANCEL ELB NEWPAGE FEED 2 TEXT 'X0-VECTOR ; PRELIMINARY APPROXIMATE VALUES OF THE + UNKNOWN PARAMETERS ***FIRST INITIAL APPROXIMATION***',2 FEED 1 TEXT 'NOTE :- THE LAST 3-ELEMENTS OF THIS VECTOR ARE IN M+ ETERS AND THE REMAINING ARE IN CYCLES PER 10 MINUTES',2 WRITE X0,FORMAT=A10 * * ******************************************************************************** * * ALLOCATE STORAGE SPACE FOR THE COMING RESULTS (A , W MATRICES) FROM * THE FORTRAN-SUBROUTINE 'DESIGN' * ALLOCATE A,(NAROWS,NACOLS),2 ALLOCATE W,(NAROWS,1),2 * * CALLING THE FORTRAN-SUBROUTINE 'DESIGN' TO ALLOCATE ITS OWN STORAGE SPACE * AND INITIALIZE ZERO VALUES ;FOR THE DESIGN MATRICES (A , W) ,MOREOVER * 'DESIGN' WILL READ-IN AND STORE THE (XYZ) SATELLITE COORDINATES -WITH THE * CORRESPONDING DOPPLERS FOR ALL THE PASSES BEING PROCESSED ( WHICH WERE * COMING FROM THE PUNCHED OUTPUT OF THE SEPARATE PROGRAM 'PRE-BATCH' * CALL DESIGN,(NPASS,NAROWS,NACOLS,X0,A,W),F * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- * * ***STARTING THE LOOP OF 'I T E R A T I O N'*** * LOOP ITERATE,ICOUNT,1,5 COPY ICOUNT,ITER NEWPAGE FEED 1 TEXT '***** RESULTS AFTER I T E R A T I O N NO. #ITE+ R *****',C * * CALLING THE ENTRY-POINT 'AANDW' OF THE SUBROUTINE 'DESIGN' TO EVALUATE * (AND TRANSFERE BACK TO 'BATCH') THE COEFFICIENT MATRIX-A AND THE MISCLOSURE * VECTOR-W (USING THE PREVIOUSLY STORED SAT. COORD. WITH DOPPLERS ,AND THE * INSERTED UP-DATED APPROXIMATE VALUES OF THE VECTOR-X0) ;TO PROCEED INTO * THE ITERATIVE LEAST SQUARES ADJUSTMENT PROCESS. * CALL AANDW,(X0,A,W),F * FEED 1 TEXT 'A-MATRIX ; THE COEFFICIENTS OF THE PARAMETERS',C WRITE A,FORMAT=A9 NEWPAGE FEED 2 TEXT 'W-VECTOR : MISCLOSURE CONSTANT TERMS (UNITS ARE IN + DOPPLER COUNTS)',C WRITE W,FORMAT=A9 * ******************************************************************************** * * FORMULATION AND SOLUTION OF THE NORMAL EQUATIONS SYSTEM FOR THE 'X'-VECTOR. * NEWPAGE TEXT 'FORMATION OF NORMAL EQNS. MATRIX (N) BEGAN AT',T TEXT 'TIME SINCE EXECUTION IS',D TRANS A,AT MULT AT,P,ATP * * MULT ATP,A,N TEXT 'MATRIX (N) WAS FORMED AFTER',D FEED 2 RDIM N,JJ TEXT 'N-MATRIX : THE COEFFICIENT MATRIX OF THE NORMAL EQU+ ATIONS',C WRITE N,FORMAT=C5 NEWPAGE FEED 2 TEXT 'OVERVIEW OF THE NORMAL EQUATIONS MATRIX POPULATION'+ ,10 TEXT '--------------------------------------------------'+ ,10 PATTERN N MULT SCALE,N,N NEWPAGE FEED 2 * ******************************************************************************** * * CALLING THE SUBPROGRAM 'SNINV' TO GET THE INVERSE OF NORMAL EQNS. MATRIX * TEXT 'DETAILS OF TIME INFOMATION REGARDING THE INVERSION + PROCESS EXECUTION',5 TEXT '---------------------------------------------------+ -----------------',5 FEED 2 TEXT 'INVERSION PROCESS BEGAN AT',T FEED 2 CALL SNINV,(HOWINV,C1,RC1,INC,N,JJ,NINVERSE) FEED 2 TEXT 'INVERSION PROCESS ENDED AT',T * * A CHECK ON THE PRECISION OF THE CHOSEN AND USED **INVERSION PROCESS** * SIMPLY BY MULTIPLYING THE NORMAL EQUATIONS MATRIX BY ITS INVERSE - THE * RESULT SHOULD BE THEORITICALLY EQUAL TO 'IDENTITY' MATRIX * MULT N,NINVERSE,INVCHECK CANCEL N NEWPAGE FEED 2 TEXT 'CHECK ON THE ACCURACY OF THE INVERSION PROCESS',C FEED 1 TEXT 'INVCHECK MATRIX : = N * NINVERSE = : IDENTITY MAT+ RIX',C FEED 1 WRITE INVCHECK,FORMAT=C5 CANCEL INVCHECK * ******************************************************************************** * NEWPAGE FEED 3 * * TEXT 'PRINTOUT OF THE REQUIRED OUTPUT FROM THE PARAMETRIC+ LEAST SQUARES ADJUSTMENT',C TEXT '===================================================+ =========================',C MULT SCALE,NINVERSE,QXA CANCEL NINVERSE FEED 2 TEXT 'QXA-MATRIX ; THE WEIGHT COEFFICIENT MATRIX OF THE + ADJUSTED PARAMETERS',C TEXT '(WHICH IS THE INVERSE OF THE NORMAL EQUATIONS MATRI+ X)',C WRITE QXA,FORMAT=C5 MULT ATP,W,U NEWPAGE FEED 2 TEXT 'U-VECTOR : THE ABSOLUTE TERMS OF THE NORMAL EQS.',C WRITE U,FORMAT=A9 MULT QXA,U,X COPYC X,X NEWPAGE FEED 2 TEXT 'X-VECTOR ; CORRECTIONS FOR THE APPROXIMATE VALUES + OF THE UNKNOWN PARAMETERS',10 FEED 1 TEXT 'NOTE :- THE LAST 3-ELEMENTS OF THIS VECTOR ARE IN M+ ETERS AND THE REMAINING ARE IN CYCLES PER 10 MINUTES',2 WRITE X,FORMAT=A10 * ******************************************************************************** * * COMPUTATION OF THE ADJUSTED VALUES ( ESTIMATES ) OF THE UNKNOWN PARAMETERS * ADD X0,X,XA CANCEL X0 NEWPAGE FEED 2 TEXT 'XA-VECTOR ; THE ADJUSTED PARAMETERS',10 FEED 1 TEXT 'NOTE :- THE LAST 3-ELEMENTS OF THIS VECTOR ARE IN M+ ETERS AND THE REMAINING ARE IN CYCLES PER 10 MINUTES',2 WRITE XA,FORMAT=A10 * ******************************************************************************** * * COMPUTATION OF THE NUMBER OF DEGREES OF FREEDOM - DF * RDIM A,RDF CDIM A,CDF SUB RDF,CDF,DF * ******************************************************************************** * * CALCULATION OF THE RESIDUALS ( V'S ) * * MULT A,X,AX ADD AX,W,V NEWPAGE FEED 2 TEXT 'V-VECTOR ; RESIDUALS OR CORRECTIONS FOR THE OBSERVE+ D QUANTITIES -UNITS ARE IN DOPPLER COUNTS (CYCLES)',2 WRITE V,FORMAT=A9 * ******************************************************************************** * * A CHECK ON THE THE COMPUTATIONS IS DONE TO VERIFY THAT : A'PV = 0 * MULT ATP,V,CHECK NEWPAGE FEED 2 TEXT 'ATPV-MATRIX ; MUST EQUAL TO A NULL VECTOR AS A CHEC+ K ON THE CORRECTNESS OF OUR COMPUTATIONS',3 WRITE CHECK,FORMAT=A9 CANCEL CHECK * ******************************************************************************** * * CALCULATION OF THE QUADRATIC FORM OF THE WEIGHTED RESIDUALS ;BY -3- DIFFERENT * METHODS ( AS A CHECK ) ,WHICH ARE : * 1-EXPLICITLY BY DIRECT MULTIPLICATION ( V'PV = V' * P * V ) * 2-IMPLICITLY BY USING ANOTHER MATRICES HAVE BEEN ALLREADY COMPUTED DURING * THE ADJUSTMENT SET-UP ,NAMELY ( V'PV = X'U + W'PW ) * 3-SAME AS MENTIONED IN (2) ; BUT USING ANOTHER MATRICES ,NAMELY : * V'PV = -K'*W = V'P * W * TRANS V,VT MULT VT,P,VTP MULT VTP,V,VTPV1 TRANS X,XT MULT XT,U,XTU TRANS W,WT MULT WT,P,WTP MULT WTP,W,WTPW ADD XTU,WTPW,VTPV2 MULT VTP,W,VTPV3 * * COMPUTE THE DIFFERENCES OF (V'PV) CALCULATED THROUGH METHODS (2) & (3) * FROM THE ORIGINAL EXPLICIT METHOD (1) * SUB VTPV1,VTPV2,DIFVPV2 SUB VTPV1,VTPV3,DIFVPV3 NEWPAGE FEED 1 TEXT 'VTPV-SCALAR ; THE QUADRATIC FORM OF THE WEIGHTED RE+ SIDUALS COMPUTED BY -3- DIFFERENT METHODS FOR CHECK',2 FEED 2 TEXT '1-VTPV1 : IS COMPUTED EXPLICITLY BY DIRECT MULTIPLI+ CATIONS :- ( V''PV = V'' * P * V )',5 * * FEED 1 TEXT '2-VTPV2 : IS COMPUTED IMPLICITLY USING THE DERIVED + FORMULA :- ( V''PV = X''U + W''PW )',5 FEED 1 TEXT '3-VTPV3 : IS COMPUTED IMPLICITLY USING THE DERIVED + FORMULA :- ( V''PV = -K''W = V''P * W )',5 FEED 1 TEXT 'DIFVPV2 : IS THE DIFFERENCE OF THE (FIRST-SECOND) M+ ETHODS USED TO COMPUTE V''PV',5 FEED 1 TEXT 'DIFVPV3 : IS THE DIFFERENCE OF THE (FIRST-THIRD ) M+ ETHODS USED TO COMPUTE V''PV',5 FEED 1 WRITE (VTPV1,VTPV2,VTPV3,DIFVPV2,DIFVPV3),FORMAT=A9 FEED 2 TEXT 'DF-SCALAR : THE NUMBER OF DEGREES OF FREEDOM',C WRITE DF,FORMAT=A5 * ******************************************************************************** * * COMPUTE THE ESTIMATED VARIANCE-FACTOR EVF ; EVF = V'PV / DF * * A TEST TO AVOID THE OVERFLOW IN COMPUTING THE ESTIMATED VARIANCE FACTOR * ,IF THE NUMBER OF THE DEGREES OF FREEDOM IS EQUAL TO ZERO * FOR (DF,NE,0),SM1 FEED 2 TEXT '*** N O R E D U N D A N T - O B S E R V A T I O N+ S ***',5 GOTO ARRANGE SM1 DIV DF,VTPV1,EVF FEED 2 TEXT 'EVF-SCALAR ; THE ESTIMATED VARIANCE FACTOR',C WRITE EVF,FORMAT=A5 * ******************************************************************************** * * COMPUTE THE MAXIMUM ABSOLUTE VALUES OF THE SOLUTION VECTOR (X) ;W.R.T. THE * STATION COORDINATES ,AND THE FREQUENCY OFFSETS --SEPARATELY. * ADD NPASS,1,NCOLX MAXMOD X,(NCOLX,1),(3,1),(IXYZ,JXYZ),MAXXYZ MAXMOD X,(1,1),(NPASS,1),(IFREQ,JFREQ),MAXFREQ * * TEST THE CONDITION OF STOPPING THE ITERATION PROCESS ;BY CHECKING THE ABOVE * MAXIMUM ABSOLUTE VALUES (MAXXYZ , MAXFREQ) AGAINST THE DESIRED GIVEN * PRECISION LIMITS (PRECXYZ , PRECFREQ). * FOR (MAXXYZ,GT,PRECXYZ),MOREITER FOR (MAXFREQ,GT,PRECFREQ),MOREITER FEED 2 GOTO CONVERGE * * * * UP-DATE THE APPROXIMATE VALUES OF THE UNKNOWN PARAMETERS (X0) ,FOR * PERFORMING A NEW ITERATION. * MOREITER RENAME (XA,X0) CANCEL W ALLOCATE W,(NAROWS,1),2 ITERATE LOOPEND * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- * FEED 2 GOTO NOCONVER * * PRINT THE NUMBER AND STATUS OF THE ITERATIONS BEING DONE. * CONVERGE TEXT ' *** SOLUTION C O N V E R G E S AFTER #ITER IT+ ERATIONS ***',5 GOTO CONTINUE NOCONVER TEXT '*** SOLUTION D O E S N O T C O N V E R G E AFTE+ R #ITER ITERATIONS ***',5 * ******************************************************************************** * CONTINUE FEED 1 * * CONTINUE-ON TO COMPUTE THE VARIANCE-COVARIANCE MATRICES ,THE FREQUENCY * OFFSETS (IN CYCLES/SEC.) ,AND THE (LAT,LONG,HT) OF THE TRACKING STATION - * -USING THE RESULTS OF THE FINAL ITERATION ONLY. * NEWPAGE FEED 2 TEXT 'DEGREES OF FREEDOM = #DF',10 FEED 2 * * COMPUTE AND PRINT THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED PARAMETERS * MULT EVF,QXA,EXA TEXT 'EXA-MATRIX ; THE VARIANCE-COVARIANCE MATRIX OF THE + ADJUSTED PARAMETERS',C WRITE EXA,FORMAT=C5 * * COMPUTE & PRINT THE WEIGHT-COEFF. MATRIX PF THE ADJUSTED OBSERVATIONS * FOR (ICPQLA,EQ,1),YESQLA FOR (ICPELA,NE,1),NOELA YESQLA MULT A,QXA,AQX CANCEL A MULT AQX,AT,QLA CANCEL AT,AQX FOR (ICPQLA,NE,1),NOQLA NEWPAGE FEED 2 * * * TEXT 'QLA-MATRIX ; THE WEIGHT COEFFICIENT MATRIX OF THE + ADJUSTED OBSERVATIONS',C WRITE QLA,FORMAT=C5 NOQLA FOR (ICPELA,NE,1),NOELA * * COMPUTE & PRINT THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED OBSERVATIONS * MULT EVF,QLA,ELA NEWPAGE FEED 2 TEXT 'ELA-MATRIX ; THE VARIANCE-COVARIANCE MATRIX OF THE + ADJUSTED OBSERVATIONS',C WRITE ELA,FORMAT=C5 CANCEL QLA,ELA NOELA FEED 1 * ******************************************************************************** * * ARRANGMENT OF THE VECTOR OF THE ADJUSTED UNKNOWN PARAMETERS ;BY SPLITTING IT * INTO TWO-GROUPS -THE FIRST FOR THE COORDINATES OF THE TRACKING STATION ,AND * THE SECOND FOR THE FREQUENCY OFFSETS OF THE PROCESSED ACCEPTED PASSES. * ARRANGE FEED 1 EXSUBM XA,(NCOLX,1),(3,1),XYZVEC EXSUBM XA,(1,1),(NPASS,1),FREQP10M DIV 600.,FREQP10M,FREQPS * * EXTRACT THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED (X,Y,Z) COORD. * OF THE TRACKING STATION FROM THE WHOLE EXA-MATRIX. * EXSUBM EXA,(NCOLX,NCOLX),(3,3),EXTRACK * * COMPUTE THE STANDARD DEVIATIONS IN THE ADJUSTED (X,Y,Z) COORDINATES. * EXDIAG EXTRACK,XYZVAR1 ROWSUM XYZVAR1,XYZVAR ESQRT XYZVAR,XYZSTDEV * * CALLING THE FORTRAN-SUBROUTINE 'RECORD' TO STORE THE ADJUSTED FREQUENCY * OFFSETS ,AND TO READ & STORE THE IDENTIFICATION CARDS WITH THE MAIN * CHARATERISTICS OF THE USED PASSES , NECESSARY FOR THE BOOK-KEEPING RECORD * CALL RECORD,(IYEAR,NPASS,FREQPS),F * * CALLING THE FORTRAN-SUBROUTINE 'FINRES' ,TO CONVERT THE ADJUSTED GEOCENTRIC * (X,Y,Z) COORD. OF THE TRACKING STATION TO (LAT,LONG & HT). * DOUBLE XYZVEC CALL FINRES,(XYZVEC),F * ******************************************************************************** * * COMPUTE THE ERROR ELLIPSOID FOR THE ADJUSTED COORD. OF THE TRACKING STN. * * * COPY EXTRACK,EXYZ ALLOCATE LAM123,(3,3),2 * * CALLING THE SUBROUTINE 'EIGEN3' TO DO THE EIGEN VALUE PROBLEM ON THE * (3,3) VARIANCE-COVARIANCE MATRIX (EXTRACK) OF THE FINAL ADJUSTED (X,Y,Z) * COORD. OF THE TRACKING STATION ; RESULTING A (3,3) DIAGONAL MATRIX (LAM123) * WHOSE DIAGONAL ELEMENTS ARE THE EIGEN VALUES ,OR THE SQUARES OF THE SEMI- * AXES (MAJOR ,INTERMEDIATE ,AND MINOR) OF THE STANDARD ERROR ELLIPSOID. * CALL EIGEN3,(EXYZ,LAM123),F ROWSUM LAM123,SEMISQ ESQRT SEMISQ,SEMIAXES * ******************************************************************************** * * PRINT-OUT THE SUMMARY OF THE FINAL RESULTS AFTER THE ITERATIVE LEAST * SQUARES ADJUSTMENT . * LOOP PRINT,ICOPY,1,NCOPIES NEWPAGE FEED 15 TEXT '** C O P Y NO. #ICOPY **',C FEED 2 TEXT '** S U M M A R Y OF THE F I N A L RESULTS **',C FEED 2 TEXT '** AFTER THE S I M U L T A N E O U S LEAST SQUARE+ S ADJUSTMENT **',C FEED 3 TEXT '***************',C FEED 3 TEXT 'NOTE : THE NUMBER OF DEGREES OF FREEDOM HERE ARE = + #DF FROM #NPASS PASSES BEING PROCESSED',C * * CALLING THE ENTRY-POINT 'PRINTJ' OF THE SUBROUTINE 'RECORD' FOR PRINTING A * TABLE TO BE AS A RECORD FOR THE IDENTIFICATION AND CHARACTERISTICS OF * THE ACCEPTED PASSES ,WHICH WERE USED IN THE ADJUSTMENT TO DETERMINE * THE COORDINATES OF THE TRACKING STATION ,MOREOVER TO INCLUDE THE ADJUSTED * FREQUENCY OFFSETS IN THIS RECORD TABLE. * * CALL PRINTJ,F * ******************************************************************************** * * CALLING THE ENTRY-POINT 'PRINTI' OF THE SUBROUTINE 'FINRES' TO PRINT THE * SUMMARY OF THE FINAL ADJUSTED CARTESIAN AND CURVILINEAR COORDINATES OF * THE TRACKING STATION ,WITH THE REFERENCE DATUM PARAMETERS. * CALL PRINTI,F * * PRINT THE VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED (X,Y & Z) COORD. * OF THE TRACKING STATION. * * * NEWPAGE FEED 2 TEXT 'EXTRACK-MATRIX : THE VARIANCE-COVARIANCE MATRIX OF + THE ADJUSTED (X,Y & Z) COORD. OF THE TRACKING STATION',2 TEXT '(UNITS ARE IN : METERS SQUARED )',C WRITE EXTRACK,FORMAT=A6 * * PRINT THE STANDARD DEVIATIONS IN THE ADJUSTED (X,Y,Z) COORD. RESPECTIVELY * FEED 3 TEXT 'THE STANDARD DEVIATION IN THE (X,Y & Z) ADJUSTED CO+ ORD. ARE AS FOLLOWS : (UNITS ARE IN METERS)',5 TEXT '---------------------------------------------------+ -------------------',5 FEED 1 WRITE XYZSTDEV,NOTITLE,FORMAT=A5 * * PRINT THE SEMI-AXES (MAJOR ,INTERMEDIATE & MINOR) OF THE STANDARD ERROR * ELLIPSOID (STANDARD CONFIDENCE REGION) * FEED 3 TEXT 'SEMI-AXES OF THE STANDARD ERROR ELLIPSOIDARE AS FOL+ LOWS : MAJOR , INTERMEDIATE & MINOR RESPECTIVELY',5 TEXT '---------------------------------------------------+ ---- (UNITS ARE IN METERS)',5 FEED 1 WRITE SEMIAXES,NOTITLE,FORMAT=A5 * FEED 2 TEXT 'EXECUTION OF THIS JOB ENDED AT',T PRINT LOOPEND * ******************************************************************************** * STOP END * * ******************************************************************************** * * * SUBPROGRAM 'SNINV' : ( M. NASSAR - DEC. 1971 ) * * ------------------ * * THE PURPOSE OF THIS SUBPROGRAM IS TO PERFORM THE INVERSE * * FOR THE COEFFICIENT MATRIX OF THE NORMAL EQUATIONS SYSTEM NEEDED FOR THE * * LEAST SQUARES ADJUSTMENT. * * * * THIS SUBPROGRAM HAS BEEN DESIGNED IN A GENERAL MANNER TO EXECUTE THE INVERS-* * ION PROCESS OF THE MATRIX BY A PRE-CHOSEN WAY ACCORDING TO THE ASSIGNING OF * * THE APPROPRIATE VALUES FOR THE OPTIONS & ARGUMENTS MENTIONED AT THE BEGINNI-* * NG OF THE MAIN PROGRAM 'BATCH' ,NAMELY IT CAN HANDLE EITHER ONE OF 3-WAYS ;* * 1-STRAIGHT-FOREWARD INVERSION OF THE WHOLE MATRIX WITHOUT ANY PARTITIONING * * 2-SUCCESSIVE PARTITIONING USING ONE-INVERSION-STEP & CALLING 'INVERT' * * 3-STRAIGHT-FOREWARD AND/OR SUCCESSIVE PATITIONING (CALLING 'INVERT') ,BY * * USING TWO-INVERSION-STEPS PROCESS * * * * THE PRE-CHOSEN WAY CAN BE DETERMINED ACCORDING TO THE DIMENSIONS OF THE * * MATRIX ,THE KIND OF THE UNKNOWN PARAMETERS ,AND THE PREVIOUS EXPERIENCE FOR * * THE COMPUTER TIME TAKEN VIA EACH WAY IN A SPECIFIED CASE * * * * THIS SUBPROGRAM WORKS WITH THE ARGUMENTS EXISTING IN ITS IDENTIFICATION * * STATEMENT (SUBPRO) ,THE CORRESPONDING OF WHICH ARE FULLY EXPLAINED THROUGH * * THE MAIN PROGRAM (BATCH) BEFORE CALLING THAT 'SNINV' * * * * THE FINAL OUTPUT OF THIS SUBPROGRAM IS THE INVERSE OF THE NORMAL EQUATIONS * * MATRIX IN FULL ,I.E. WITH ITS DIAGONAL & OFF-DIAGONAL ELEMENTS. * * * ******************************************************************************** * * IDENTIFICATION CARD OF THE SUBPROGRAM 'SNINV' * SUBPRO SNINV,(NP,C1,RC1,INC,SN,JJ,NI) ATTRIB (SN,NI),DEF=1 * * NOTE-'SN' WHEN MENTIONED ANYWHERE THROUGH THIS SUBPROGRAM ,IT INDICATES * THE NORMAL EQUATIONS MATRIX NEEDED TO BE INVERTED * * PERFORMING THE INVERSE OF 'SN' BY THE STRAIGHT-FOREWARD INVERSION WITHOUT * ANY PARTITIONING ; WHEN DESIRED * FOR (NP,NE,0),S5 TEXT 'THE STRAIGHT-FOREWARD INVERSION BEGAN AT',T TEXT 'TIME SINCE EXECUTION IS ',D FEED 1 FORMS ONEINV,(JJ,JJ),(1,1),(1,1),JJ,1 DIV SN,ONEINV,NI CANCEL ONEINV TEXT 'TIME TAKEN FOR INVERSION IS ',D GOTO S6 * * PERFORMING THE INVERSE OF 'SN' USING ONE MAIN INVERSION STEP,WHICH IS HANDLED * BY SUCCESSIVE PARTITIONING THROUGH SUBPROGRAM 'INVERT' * * * S5 FOR (NP,GT,0),S7 TEXT 'THE ONE-INVERSION-STEP PROCESS BEGAN AT',T TEXT 'TIME SINCE EXECUTION',D CALL INVERT,(C1,JJ,INC,SN,NI) TEXT 'TIME TAKEN FOR INVERSION IS',D GOTO S6 * * PERFORMING THE INVERSE OF 'SN' BY SPLITING IT INTO TWO-MAIN-INVERSION-STEPS * ,EACH OF THEM CAN BE HANDLED EITHER BY ; STRAIGHT-FOREWARD INVERSION OR * SUCCESSIVE PARTITIONING THROUGH 'INVERT' ,AS DESIRED * S7 TEXT 'THE TWO-INVERSION-STEPS PROCESS BEGAN AT',T TEXT 'TIME SINCE EXECUTION IS',D SUB JJ,NP,RD ADD 1,RD,R1 EXSUBM SN,(1,1),(RD,RD),N11 EXSUBM SN,(1,R1),(RD,NP),N12 EXSUBM SN,(R1,1),(NP,RD),N21 EXSUBM SN,(R1,R1),(NP,NP),N22 * * ALLOWING STRAIGHT-FOREWARD INVERSION (I.E. WITHOUT SUCCESSIVE PARTITIONING) * OF THE FIRST INVERSION STEP : IF DESIRED * FOR (C1,NE,0),S3 FEED 1 TEXT 'FIRST INVERSION -STEP STARTED AT',T TEXT 'TIME SINCE EXECUTION IS ',D FEED 1 FORMS KLM,(RD,RD),(1,1),(1,1),RD,1 DIV N11,KLM,NI1 GOTO S4 * * PERFORMING THE FIRST INVERSION STEP BY SUCCESSIVE PARTITIONING THROUGH * THE SUBPROGRAM 'INVERT' * S3 TEXT 'FIRST INVERSION STEP STARTED AT',T TEXT 'TIME SINCE EXECUTION IS ',D FEED 1 CALL INVERT,(C1,RD,INC,N11,NI1) S4 TEXT 'TIME FOR FIRST INVERSION IS',D TEXT 'FIRST INVERSION STEP ENDED AT',T FEED 1 MULT NI1,N12,NIN CANCEL N12 MULT N21,NI1,NNI MULT N21,NIN,N2NIN1 CANCEL N21 SUB N22,N2NIN1,NEX * NOTE- THAT (NEX) IN THE ABOVE STATMENT IS THE MODIFIED (N22) SUBMATRIX , * WHICH IS ACTUALLY REQUIRED TO BE INVERTED IN THE SECOND INVERSION STEP * CANCEL N22,N2NIN1 * * * * ALLOWING STRAIGHT-FOREWARD INVERSION (I.E. WITHOUT SUCCESSIVE PARTITIONING) * OF THE SECOND INVERSION STEP : IF DESIRED * FOR (RC1,NE,0),S1 TEXT 'SECOND INVERSION STEP BEGAN AT',T TEXT 'TIME SINCE EXECUTION IS',D FEED 1 FORMS JK,(NP,NP),(1,1),(1,1),NP,1 DIV NEX,JK,NI22 GOTO S2 * * PERFORMING THE SECOND INVERSION STEP BY SUCCESSIVE PARTITIONING THROUGH * THE SUBPROGRAM 'INVERT' * S1 TEXT 'SECOND INVERSION STEP BEGAN AT',T TEXT 'TIME SINCE EXECUTION IS ',D FEED 1 CALL INVERT,(RC1,NP,INC,NEX,NI22) S2 TEXT 'TIME FOR SECOND INVERSION IS ',D TEXT 'SECOND INVERSION STEP ENDED AT',T FEED 1 CANCEL NEX * * EVALUATION OF THE SUBMATRICES FORMING THE INVERSE OF (SN)-MATRIX. * COPYC NI22,NNI2 MULT NNI2,NNI,NI21 MULT NIN,NNI2,NI12 CANCEL NN12 MULT NIN,NI21,NEN CANCEL NIN SUB NI1,NEN,NI11 CANCEL NEN * * GENERATING & BUILDING-UP THE (NI)-MATRIX ,WHICH IS THE INVERSE OF (SN)-MATRIX * NULLMAT NI,(JJ,JJ) INSUBM NI11,NI,(1,1) INSUBM NI12,NI,(1,R1) INSUBM NI21,NI,(R1,1) INSUBM NI22,NI,(R1,R1) CANCEL NI11,NI12,NI21,NI22 TEXT 'THE TWO-INVERSION-STEPS PROCESS ENDED AT',T S6 RETURN END * * * ******************************************************************************** * * * SUBPROGRAM 'INVERT' : * * ------------------- * * THIS SUBPROGRAM WAS WRITTEN BY J.N. MENSAH (APRIL 1971) * * ITS PURPOSE IS TO PERFORM THE INVERSE OF ANY MATRIX BY SUCCESSIVE PARTITIONS* * STARTING WITH THE GIVEN SIZE OF THE UPPER LEFT SUBMATRIX (C1) AND ADDING * * THE GIVEN INCREMENT (INC) UNTIL WE GET THE TOTAL SIZE OF THE MATRIX INVERTED* * * * THE SUBPROGRAM WORKS WITH THE ARGUMENTS EXISTING IN ITS IDENTIFICATION * * STATEMENT (SUBPRO) ,THE CORRESPONDING OF WHICH ARE EXPLAINED IN FULL DETAILS* * EITHER THROUGH THE MAIN PROGRAM 'BATCH' ,OR THROUGH THE SUBPROGRAM 'SNINV' * * * * THE ONLY ADDITIONS HERE TO THE ORIGINAL FORM OF 'INVERT' ARE SOME COMMENT * * CARDS TO ILLUSTRATE THE EXECUTION PROCESS. * * * ******************************************************************************** * * IDETIFICATION CARD OF THE SUBPROGRAM 'INVERT' * SUBPRO INVERT,(C1,RD,INC,NS,VI) ATTRIB (NS,VI),SYM=1 * * INVERSION OF THE INITIAL UPPER LEFT SUBMATRIX (V11) ,WHOSE GIVEN * SIZE IS (C1) * EXSUBM NS,(1,1),(C1,C1),V11 FORMS IK,(C1,C1),(1,1),(1,1),C1,1 DIV V11,IK,VI * * UP-DATING THE SIZE OF (V11) FOR THE FIRST TIME BY ADDING THE GIVEN * INCREMENT (INC). * ADD C1,INC,D1 ADD 1,C1,C2 COPY D1,II * * STARTING THE SUCCESSIVE PARTITIONING AND UP-DATING PROCESS THROUGH * A LOOP TO END-UP WITH THE MATRIX (NS) BEING INVERTED TO (VI)-MATRIX * START LOOP END1,D1,II,RD,INC EXSUBM NS,(1,1),(D1,D1),V EXSUBM V,(1,C2),(C1,INC),V12 EXSUBM V,(C2,1),(INC,C1),V21 EXSUBM V,(C2,C2),(INC,INC),V22 MULT VI,V12,VIV * * EVALUATION & INVERTING THE MODIFIED (V22) SUBMATRIX NEEDED TO GET THE * INVERSE OF THE UP-DATED (V11) * MULT V21,VI,VVI MULT V21,VIV,V2VIV1 * * * SUB V22,V2VIV1,UV FORMS IM,(INC,INC),(1,1),(1,1),INC,1 DIV UV,IM,U22 COPYC U22,NU22 * * PERFORMING THE INVERSE (VI) OF THE UP-DATED (V11) ,WHICH WILL BE FINALLY * THE INVERSE OF (NS) MATRIX ,OBVIOUSLY AFTER THE LAST CYCLE OF THE LOOP * MULT NU22,VVI,U21 MULT VIV,NU22,U12 MULT VIV,U21,VUV SUB VI,VUV,U11 NULLMAT UI,(D1,D1) INSUBM U11,UI,(1,1) INSUBM U12,UI,(1,C2) INSUBM U21,UI,(C2,1) INSUBM U22,UI,(C2,C2) RENAME (UI,VI) * * SUCCESSIVE UP-DATING OF THE SIZE OF (V11) ,BY SUCCESSIVE ADDITION OF (INC) * ADD C1,INC,C1 ADD 1,C1,C2 * * END OF THE LOOP ,RETURN BACK AGAIN TO THE LOOP-START FOR A NEW CYCLE ,TO * EVALUATE THE INVERSE OF THE UP-DATED (V11) A NUMBER OF TIMES UNTIL WE GET * THE TOTAL SIZE OF (NS)-MATRIX BEING INVERTED. * END1 LOOPEND RETURN END C C C C******************************************************************************* C* * C* SUBROUTINE 'DESIGN' : ( M. NASSAR - DEC. 1971 ) * C* ------------------- * C* * C* THE PURPOSE OF THIS SUBROUTINE IS TO EVALUATE THE * C* DESIGN MATRICES NEEDED TO START THE L. S. ADJUSTMENT PROCESS ,NAMELY : THE * C* COEFFICIENT MATRIX OF THE UNKNOWN PARAMETERS (A-MATRIX) ,AND THE MISCLOSURE * C* OR THE ABSOLUTE TERMS OF THE OBSERVATION EQUATIONS (W-VECTOR) THIS CAN BE * C* ACHIEVED VIA THE ENTRY-POINT 'AANDW' OF THIS SUBROUTINE BY CALLING IT FROM * C* THE MAIN PROGRAM 'BATCH' AT THE BEGINNING OF EACH ITERATION. * C* * C* THE REQUIRED INPUT DATA ARE : * C* --------------------------- * C* * C* I-FOR EACH ACCEPTED PASS TO BE PROCESSED ,THE FOLLOWING IS NEEDED ; * C* 1-THE NUMBER OF OBSERVED SATELLITE POSITIONS * C* 2-THE SATELLITE (X,Y,Z) COORDINATES FOR EACH OBSERVED POSITION * C* 3-THE OBSERVED DOPPLER COUNTS ,AFTER BEING CORRECTED FOR IONOSPHERIC & * C* TROPOSPHERIC REFRACTIONS. * C*II-THE INITIAL APPROXIMATE VALUES OF THE UNKNOWN PARAMETERS -(X0) * C* * C* NO. 1,2 & 3 IN (I) ,ARE OBTAINED ON PUNCHED CARDS COMING FROM THE SEPARATE * C* FORTRAN-PROGRAM 'PRE-BATCH' ,AND INSERTED WITHIN THE DATA STREAM FOR 'BATCH'* C* ACCORDING TO THE FOLLOWING FORMAT ; * C* * C* FIRST CARD : * C* ---------- * C* COLS 3 TO 5 : NUMBER OF OBSERVED SATELLITE POSITIONS .......... I3 * C* * C* SECOND CARD : * C* ----------- * C* COLS 2 TO 19 : X-COORDINATE OF THE SATELLITE POSITION (IN METERS)..D18.11* C* COLS 22 TO 39 : Y-COORDINATE OF THE SATELLITE POSITION (IN METERS)..D18.11* C* COLS 42 TO 59 : Z-COORDINATE OF THE SATELLITE POSITION (IN METERS)..D18.11* C* COLS 62 TO 79 : CORRESPONDING DOPPLER COUNTS (IN CYCLES) ....... D18.11 * C* * C* NOTE : THE SECOND CARD WILL BE REPEATED AS MANY AS THE NUMBER OF OBSERVED * C* SATELLITE POSITIONS FROM THIS PASS (AS PUNCHED ON THE FIRST CARD) * C* * C* X0 IN NO. (II) IS OBTAINED EACH ITERATION FROM THE MAIN MATLAN-PROGRAM * C* 'BATCH' THROUGH THE CALLING STATEMENT ,OF COURSE AFTER X0 BEING UP-DATED. * C* * C* SINCE IT HAS BEEN DESIRED TO DO AN ITERATIVE SOLUTION ; THIS SUBROUTINE * C* WORKS WITH THE MAIN PROGRAM 'BATCH' AS FOLLOWS : * C* 1-JUST BEFORE THE BEGINNING OF THE ITERATIONS LOOP ,'BATCH' CALLS 'DESIGN' * C* (ONLY ONCE) TO DO THE FOLLOWING ;ALLOCATE STORAGE SPACE FOR THE MATRICES * C* A,W & X0 -READ IN AND STORE FOR EACH PASS THE NUMBER OF OBSERVED SATELLITE* C* POSITIONS ,THE VALUES OF THE SAT. (X,Y,Z) COORD. & THE CORRESPONDING * C* DOPPLER COUNTS -AND FINALLY INITIALIZE ZERO VALUES FOR THE ELEMENTS OF A * C* AND W MATRICES ,THEN RETURN BACK TO 'BATCH' * C* * C* 2-FOR EACH CURRENT ITERATION ,AS A STARTING STEP OF THE ITERATIONS LOOP , * C* 'BATCH' CALLS THE ENTRY-POINT 'AANDW' ,INSERTING TO IT THE APPROXIMATE * C* (UP-DATED) VALUES OF THE UNKNOWNS (X0) ,AND GETTING BACK A & W EVALUATED. * C C C C* THE MAIN REASON OF USING THE ENTRY-POINT HERE IS TO OVERCOME A PROBLEM , * C* NAMELY : WE NEED THE SATELLITE COORD. & DOPPLERS TO BE READ-IN ONLY ONCE * C* (AND THIS IS EXACTLY WHAT 'DESIGN' DOES) , HOWEVER WE NEED THE MATRICES A * C* & W TO BE COMPUTED AS MANY AS THE NECESSARY NUMBER OF ITERATIONS ,WHICH * C* DEPENDS ON THE DESIRED PRECISION LIMITS ,(AND THIS IS ACTUALLY WHAT THE * C* ENTRY-POINT 'AANDW' DOES ,UTILIZING THE ALREADY STORED SATELLITE COORD. & * C* DOPPLERS ,PLUS THE INSERTED UP-DATED APPROX. VALUES (X0) EACH TIME) * C* * C* THE OUTPUT WILL BE THE DESIGN MATRICES : A & W ,BEING EVALUATED * C* * C******************************************************************************* C SUBROUTINE DESIGN(NPASS,NAROWS,NACOLS,X0,A,W) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X0(NACOLS),A(NAROWS,NACOLS),W(NAROWS) DIMENSION XSAT(540),YSAT(540),ZSAT(540),DOPVAC(540),NUMPOS(60) C C INSERT THE VALUES OF SOME CONSTANTS NEEDED DURING THE EXECUTION OF C THIS ROUTINE ,NAMELY : C C : IS THE VACIUM VELOCITY OF THE LIGHT -( IN METERS/SEC. ) C F0 : IS THE CONSTANT REFERENCE FREQUENCY GENERATED IN THE RECEIVER- C UNITS ARE IN CYCLES/SEC. C DT : IS THE TIME INTERVAL OVER WHICH THE DOPPLER SHIFTS ARE INTEGRATED C TO END-UP WITH THE DOPPLER COUNTS (CYCLES) BETWEEN TWO CONSECUTIVE C SATELLITE POSITIONS ,IN OUR CASE ; DT = 2 UT MINUTES ) C C=299792500.0D 00 F0=400.0D 06 DT=2.0D 00 C C SCALING THE INTEGRATION INTERVAL (DT) BY DIVIDING ITS ABOVE VALUE BY 10 C TO BE COMPATIBLE WITH THE SCALING OF THE FREQUENCY OFFSET (BEING DONE IN C THE SEPARATE (FORTRAN) PROGRAM 'PRE-BATCH') C DT=0.20D 00 C C READ-IN THE PUNCHED DATA COMINNG FROM THE SEPARATE PROGRAM-'PRE-BATCH' C J=1 K=0 DO 10 I=1,NPASS C C READ AND STORE THE NUMBER OF OBSERVED SATELLITE POSITIONS PER EACH PASS C READ(1,20)NUMPOS(I) 20 FORMAT(2X,I3) K=K+NUMPOS(I) C C READ AND STORE THE SATELLITE (XYZ) COORD. IN THE AVER. TERR. SYSTEM ,WITH C THE CORRESPONDING DOPPLER COUNTS PER EACH PASS . C DO 30 JK=J,K READ(1,40)XSAT(JK),YSAT(JK),ZSAT(JK),DOPVAC(JK) 40 FORMAT(1X,D18.11,3(2X,D18.11)) C C C 30 CONTINUE J=K+1 10 CONTINUE C C INITIALIZE ZERO VALUES FOR THE ELEMENTS OF A-MATRIX AND W-VECTOR , C ACCORDING TO THE GIVEN DIMENSIONS (NAROWS , NACOLS) C DO 50 II=1,NAROWS W(II)=0.0 DO 50 IJ=1,NACOLS A(II,IJ)=0.0 50 CONTINUE RETURN C C******************************************************************************* C ENTRY AANDW(X0,A,W) C C EVALUATE THE DESIGN MATRIX-A AND THE MISCLOSURE VECTOR-W C NOTES :1-THE APPROXIMATE VALUES OF THE UNKNOWN PARAMETERS ARE INSERTED C HERE VIA THE CALLING MAIN PROGRAM 'BATCH' C 2-ONE DOPPLER EQUATION IS SETUP PER EACH DOPPLER COUNT INTEGRATED C BETWEEN EACH TWO CONSECUTIVE SATELLITE POSITIONS (J & K) WHERE C (K = J + 1) ,EXPLICITLY ; K IS THE CURRENT SEQUENCE NUMBER OF C THE SATELLITE POSITION AND ; J IS THE LATEST PRECEDING ONE , C W.R.T. THE TRACKING STATION. C C C COMPUTE THE CURRENT COLUMN NUMBER OF THE A-MATRIX PERTAINING TO THE C X , Y & Z COORDINATES OF THE TRACKING STATION. C NCOLX=NPASS+1 NCOLY=NPASS+2 NCOLZ=NPASS+3 M=1 N=0 NROW=0 DO 60 L=1,NPASS C C COMPUTE THE CURRENT COLUMN NUMBER OF THE A-MATRIX PERTAINING TO THE C FREQUENCY OFFSET OF THE PASS # L - (NCOLF) C NCOLF=L N=N+NUMPOS(L) DO 70 MN=M,N C C COMPUTE THE RANGES BETWEEN THE TRACKING STATION (I) AND THE CURRENT C SATELLITE POSITION (K) C DXIK=XSAT(MN)-X0(NCOLX) DYIK=YSAT(MN)-X0(NCOLY) DZIK=ZSAT(MN)-X0(NCOLZ) DRIK=DSQRT(DXIK*DXIK+DYIK*DYIK+DZIK*DZIK) C C COMPUTE THE CONTRIBUTION OF THE CURRENT POSITION 'K' TO THE COEFFICIENTS C OF THE A-MATRIX W.R.T. THE TRACKING STATION COORDINATES. C C C AXIK=DXIK/DRIK AYIK=DYIK/DRIK AZIK=DZIK/DRIK C C CHECK TO VERIFY THE EXISTENCE OF A DOPPLER COUNT ; IF NO ,I.E. WHEN C DOPVAC(MN)=0 - NEGLECT THIS OBSERVATION EQUATION ,AND JUMP TO FORM C ANOTHER ONE C IF(DOPVAC(MN).EQ.0) GO TO 80 C C UP-DATE THE CURRENT NUMBER OF THE DOPPLER OBSERVATION EQUATION ; C (WHICH IS THE SAME ROW NUMBER OF THE A-MATRIX AT THIS POINT) C NROW=NROW+1 C C COMPUTE THE NON-ZERO ELEMENTS OF THE A-MATRIX ; NOTE THAT THE FIRST THREE C COLUMNS PERTAIN TO THE COORD. OF TRACKING STATION ,AND THE REMAINING C COLUMNS PERTAIN TO THE FREQUENCY OFFSETS (AS MANY AS THE PROCESSED PASSES) C A(NROW,NCOLX)=(F0/C)*(AXIJ-AXIK) A(NROW,NCOLY)=(F0/C)*(AYIJ-AYIK) A(NROW,NCOLZ)=(F0/C)*(AZIJ-AZIK) A(NROW,NCOLF)=DT C C COMPUTE THE ELEMENTS OF THE MISCLOSURE VECTOR-W ; W IS THE DIFFERENCE C BETWEEN THE APPROXIMATE COMPUTED DOPPLER COUNT (FROM THE APPROX.VALUES OF C THE UNKNOWN PARAMETERS) ,AND THE OBSERVED COUNT (AFTER BEING CORRECTED C FOR THE IONOSPHERIC & TROPOSPHERIC REFRACTIONS) ; W = APROXN - DOPVAC C DRKIJ=DRIK-DRIJ APROXN=X0(NCOLF)*DT+(F0/C)*DRKIJ W(NROW)=APROXN-DOPVAC(MN) C C INITIALIZE CONDITIONS FOR FORMING A NEW DOPPLER EQUATION BY SETTING THE C DATA & INFORMATIONS ALREADY COMPUTED AT THE SATELLITE POSITION 'K' IN C THE LAST DOPPLER EQUATION ,TO BE THE SAME VALUES FOR THE POSITION 'J' C IN THE NEW EQUATION ,THEN PROCEED THE COMPUTATIONS AT A NEW POSITION 'K' C 80 AXIJ=AXIK AYIJ=AYIK AZIJ=AZIK DRIJ=DRIK 70 CONTINUE M=N+1 60 CONTINUE RETURN END C C******************************************************************************* C* * C* SUBROUTINE 'RECORD' : ( M. NASSAR - JAN. 1972 ) * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO ARRANGE AND * C* PRINT-OUT A COMPLETE TABLE AS A BOOK-KEEPING RECORD FOR THE IDENTIFICATION * C* AND THE MAIN CHARACTERISTICS (W.R.T. THE TRACKING STATION) OF ALL THE * C* ACCEPTED PASSES ,WHICH WERE USED IN THE LEAST SQUARES ADJUSTMENT PROCESS * C* TO DETERMINE THE COORDINATES OF THE TRACKING STATION. * C* * C* THIS WILL BE ACHIEVED VIA THE ENTRY-POINT 'PRINTJ' OF THIS SUBROUTINE * C* ,WHICH ENABLES US TO HAVE AS MANY COPIES OF THIS TABLE AS DESIRED * C* (ACCORDING TO THE ASSIGNED VALUE FOR (NCOPIES) IN THE INPUT DATA), * C* * C* THE REQUIRED INPUT DATA ARE : * C* --------------------------- * C* * C* 1-FOR EACH PASS ; THE IDEN. CARD ,THE MAX. ELEVATION ,TRACK DIRECTION ,THE * C* AZIMUTH ,AND THE QUADRANT. * C* THIS CAN BE EASILY OBTAINED ON PUNCHED CARDS COMING FROM THE OUTPUT OF * C* THE 'PRE-BATCH' PROGRAM ,AND INSERTED WITHIN THE INPUT DATA STREAM OF * C* 'BATCH' PROGRAM ; ACCORDING TO THE FOLLOWING FORMAT :- * C* * C* COLS 11 TO 22 : THE PASS IDENTIFICATION CARD ......... I7 AND I5 * C* COLS 28 TO 32 : THE MAXIMUM ELEVATION (IN DEGREES) .......... F5.2 * C* COLS 38 TO 40 : THE DIRECTION OF THE TRACK ................ A3 * C* COLS 46 TO 49 : THE AZIMUTH (EITHER EAST OR WEST) .......... A4 * C* COLS 55 TO 57 : THE QUADRANT OF THE PASS .............. A3 * C* * C* 2-THE ESTIMATED VALUES OF THE FREQUENCY OFFSETS FROM THE RESULTS OF THE * C* LEAST SQUARES ADJUSTMENT ,FOR ALL THE ACCEPTED PASSES. * C* THIS CAN BE EASILY OBTAINED FROM THE CALLING MAIN-PROGRAM 'BATCH' ,VIA * C* THE CALL STATEMENT. * C* * C* THE OUTPUT WILL BE THE BOOK-KEEPING RECORDING TABLE CONTAINING : SEQUENCE * C* # ,DAY # ,LOCKON TIME (HOUR & MIN) ,PASS # ,SATELLITE # ,MAX. ELEVATION , * C* TRACK DIRECTION ,AZIMUTH ,QUADRANT ,# OF ACCEPTED DOPPLER EQUATIONS ,AND * C* FINALLY THE ADJUSTED VALUES OF THE FREQUENCY OFFSETS (IN CYCLES/SEC) * C* * C* THIS TABLE WILL BE SPLITTED INTO GROUPS ,EACH OF FOUR PASSES ,TO IMPROVE * C* THE READABILITY. * C* * C******************************************************************************* C SUBROUTINE RECORD(IYEAR,NPASS,FREQPS) REAL*8 FREQPS(NPASS) REAL VAMAX(100) DIMENSION IDL(100),ISPE(100),IDIR(100),IAZ(100),IQUADR(100) DIMENSION IDAY(100),ILHR(100),ILMIN(100),IPASS(100),ISAT(100) DIMENSION IEQNS(100) C C READ-IN THE PASS IDENTIFICATION CARD AND ITS CHARACTERISTICS W.R.T. THE C TRACKING STATION ,FROM THE PUNCHED CARDS COMING FROM THE SEPARATE - C PROGRAM 'PRE-BATCH' C DO 50 I=1,NPASS READ(1,51)IDL(I),ISPE(I),VAMAX(I),IDIR(I),IAZ(I),IQUADR(I) C C 51 FORMAT(10X,I7,I5,5X,F5.2,5X,A3,5X,A4,5X,A3) C C SPLITTING THE PASS IDENTIFICATION CARD TO ITS ORIGINAL COMPONENTS ,WHICH C ARE : DAY #,LOCKON TIME (HR & MIN) ,SAT # ,PASS # ,AND # OF DOPPLERS C IDAY(I)=IDL(I)/10000 ILHM=IDL(I)-IDAY(I)*10000 ILHR(I)=ILHM/100 ILMIN(I)=ILHM-ILHR(I)*100 ISAT(I)=ISPE(I)/1000 IPE=ISPE(I)-ISAT(I)*1000 IPASS(I)=IPE/10 IEQNS(I)=IPE-IPASS(I)*10 50 CONTINUE RETURN C C******************************************************************************* C C *** PRINT-OUT THE PASS IDENTIFICATION BOOK-KEEPING RECORD *** C ENTRY PRINTJ C WRITE(3,10)NPASS 10 FORMAT('1'//10X,'TOTAL NUMBER OF THE ACCEPTED PASSES IS =',I5/) WRITE(3,70)IYEAR 70 FORMAT(25X,'YEAR OF OBSERVATIONS WAS : ',I4/) WRITE(3,20) 20 FORMAT(20X,'IDENTIFICATION AND MAIN CHARACTERISTICS RECORD OF THE 1ACCEPTED PASSES'/20X,'============================================ 2========================='/) WRITE(3,30) 30 FORMAT(2X,'SEQUENCE',2X,'DAY',3X,'LOCKON-TIME',2X,'PASS',2X,'SATEL 1LITE',1X,'MAX.-ELEVATION',3X,'TRACK',3X,'AZIMUTH',1X,'QUADRANT',1X 2,'DOPPLER',2X,'FREQ.-OFFSET'/3X,'NUMBER',2X,'NUMBER',2X,'HOUR',2X, 3'MIN',2X,'NUMBER',2X,'NUMBER',4X,'(IN DEGREES)',2X,'DIRECTION',19X 4,'EQNS.',3X,'(CYCLES/SEC)'/) C C C PRINT-OUT FOR EACH ACCEPTED PASS ; ITS IDENTIFICATION AND THE MAIN C CHARACTERISTICS INCLUDING THE ADJUSTED FREQUENCY OFFSET ,,W.R.T. THE C TRACKING STATION - AND ARRANGE THESE PASSES IN GROUPS OF FOUR. C JCOUNT=0 DO 60 I=1,NPASS WRITE(3,52)I,IDAY(I),ILHR(I),ILMIN(I),IPASS(I),ISAT(I),VAMAX(I),ID 1IR(I),IAZ(I),IQUADR(I),IEQNS(I),FREQPS(I) 52 FORMAT(4X,I3,5X,I3,5X,I2,4X,I2,4X,I2,6X,I2,9X,F5.2,9X,A3,6X,A4,5X, 1A3,6X,I1,5X,F12.6) JCOUNT=JCOUNT+1 IF(JCOUNT.NE.4) GO TO 60 WRITE(3,54) 54 FORMAT(/) JCOUNT=0 60 CONTINUE RETURN END C C C******************************************************************************* C* * C* * C* SUBROUTINE 'FINRES' : ( M. NASSAR - DEC. 1971 ) * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO CONVERT THE FINAL * C* ADJUSTED GEOCENTRIC CARTESIAN COORDINATES (X,Y,Z) OF THE TRACKING STATION * C* OBTAINED FROM THE ITERATIVE BATCH PARAMETRIC L. S. ADJUSTMENT ; TO ITS * C* CORRESPONDING GEOCENTRIC GEOGRAPHICAL COORDINATES (LAT,LONG,HT) WITH RESPECT* C* TO THE GIVEN REFERENCE ELLIPSOID. * C* * C* ON THE OTHER HAND ,IF THE DATUM TRANSLATIONS OF THAT REFERENCE ELLIPSOID * C* FROM THE GEOCENTRE WERE AVAILABLE AND INPUTTED WITH THE DATA INPUT STREAM * C* THIS SUBROUTINE WOULD CALCULATE THE GEODETIC CARTESIAN COORDINATES (X,Y,Z) * C* BY APPLYING THE DATUM TRANSLATIONS ,AND AFTER THAT CONVERT THESE (XYZ) TO * C* THEIR CORRESPONDING GEODETIC COORDINATES (LAT,LONG,HT). * C* * C* MOREOVER ,IF THE ELLIPSOIDAL PARAMETERS OF ANY DESIRED ARBITRARILY CHOSEN * C* GEOCENTRIC ELLIPSOID ARE INPUTTED TO THIS SUBROUTINE ; THE CORRESPONDING * C* (LAT,LONG,HT) OF THE FINAL ADJUSTED (X,Y,Z) COORDINATES OF THE TRACKING * C* STATION WILL BE EVALUATED ,W.R.T. THIS GEOCENTRIC ELLIPSOID. * C* * C* FINALLY ; BOTH THE GIVEN AND COMPUTED PARAMETERS OF THE REFERENCE ELLIPSOID * C* ,THE DATUM TRANSLATIONS ,AND ALL THE ABOVE MENTIONED SETS OF COORDINATES * C* WILL BE SHOWN ON THE PRINTED OUTPUT -INTO NICELY ARRANGED TABLES WITH THE * C* NECESSARY DESCRIPTION (WHEN NEEDED) FOR AN EASIER READOUT. THIS WILL BE * C* ACHIEVED VIA THE ENTRY-POINT 'PRINTI' OF THIS SUBROUTINE. THE USE OF THE * C* ENTRY-POINT HERE WILL PERMIT US TO PRINT THESE RESULTS AS MANY COPIES AS * C* REQUIRED (ACCORDING TO THE ASSIGNED VALUE FOR (NCOPIES) IN THE INPUT DATA). * C* * C* THE REQUIRED INPUT DATA ARE : * C* ---------------------------- * C* * C* 1-THE VECTOR OF FINAL ADJUSTED GEOCENTRIC CARTESIAN COORDINATES (XYZ) OF * C* THE TRACKING STATION ,COMING FROM THE CALLING MAIN PROGRAM 'BATCH'. * C* 2-THREE OR FOUR PUNCHED CARDS ,COMING FROM THE SEPARATE FORTRAN PROGRAM * C* 'PRE-BATCH' ,THESE CARDS SHOULD CONTAIN THE FOLLOWING : * C* * C* FIRST CARD : * C* ---------- * C* COLS 3 TO 5 : GEODETIC ELLIPSOID NAME .................... A3 * C* COLS 8 TO 10 : GEOCENTRIC ELLIPSOID NAME ................... A3 * C* * C* SECOND CARD : * C* ----------- * C* COLS 2 TO 11 : SEMI MAJOR AXIS OF THE GEODETIC REFERENCE ELLIPSOID * C* -(IN METERS) ............................ F10.2 * C* COLS 14 TO 23 : RECIPROCAL OF THE FLATTENNING -(UNITLESS) ...... F10.6 * C* COLS 26 TO 32 : X-COMPONENT OF TRANSLATIONS -IN (METERS) ...... F7.2 * C* COLS 35 TO 41 : Y-COMPONENT OF TRANSLATIONS -IN (METERS) ...... F7.2 * C* COLS 44 TO 50 : Z-COMPONENT OF TRANSLATIONS -IN (METERS) ...... F7.2 * C* COL 55 : IDENTIFICATION OF DATUM TRANSLATIONS EXISTANCE .... I1 * C* ,WHICH SHOULD BE 1 ; IF ANYONE OF THE THREE TRANSLATION* C* COMPONENTS (XTRANS,YTRANS OR ZTRANS) IS NOT EQUAL ZERO.* C* NOTE : IF ALL THE TRANSLATION COMPONENTS ARE EQUAL TO ZERO ,THE FIRST * C* CARD FROM COL. 21 ON ,WILL BE AS IT WAS LEFT BLANK. * C C C* * C* THIRD CARD : ( IF ANY ) * C* ---------- * C* COLS 2 TO 11 : SEMI-MAJOR AXIS OF THE GEOCENTRIC ELLIPSOID ..... F10.2 * C* COLS 14 TO 23 : GEOCENTRIC-ELLIPSOID ,RECIPROCAL OF FLATTENNING...F10.6 * C* * C* FOURTH CARD : * C* ----------- * C* COLS 2 TO 5 : TRACKING STATION ABBREVIATED NAME ...... A4 * C* COLS 8 TO 9 : TRACKING STATION NUMBER ................... I2 * C* * C* THE OUTPUT WILL BE : * C* ------------------ * C* * C* 1-THE GEOCENTRIC GEOGRAPHICAL COORDINATES ,W.R.T. THE GEODETIC ELLIPSOID ;* C* LATITUDE .... IN DEG , MIN , SEC. * C* LONGITUDE.... IN DEG , MIN , SEC. * C* GEOMETRICAL HEIGHT ABOVE THE REFERENCE ELLIPSOID ....... IN METERS. * C* * C* 2-IF DESIRED ; WE CAN APLLY THE DATUM TRANSLATIONS AS EXPLAINED EARLIER , * C* TO GET - THE GEODETIC ( LAT , LONG & HT ) IN THE SAME UNITS AS NUMBER 1 * C* * C* 3-IF REQUIRED ,ANOTHER SET OF GEOCENTRIC (LAT,LONG,HT) W.R.T. ANY CHOSEN * C* GEOCENTRIC ELLIPSOID ,CAN BE OBTAINED IN THE SAME UNITS AS ABOVE. * C* * C* THE USED SUBROUTINES WITH THIS ROUTINE ARE : * C* ------------------------------------------ * C* 1-'ABOESA' ;TO CONVERT (X,Y,Z) TO (LAT,LONG)-IN RADIANS & (HT)-IN METERS * C* WHERE THE LONG. IS COMPUTED VIA THE SUBROUTINE 'MADAD' * C* 2-'SIMSIM' :TO CONVERT THE (LAT,LONG) FROM RADIANS TO ; DEG , MIN & SEC. * C* * C******************************************************************************* C SUBROUTINE FINRES(XAVEC) IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION DATAN,DSIN,DCOS,DABS,DSQRT DIMENSION XAVEC(3) INTEGER ELNAME/'NIL'/,CGELIP PI=3.141592653589793 RHO=180.*3600./PI X=XAVEC(1) Y=XAVEC(2) Z=XAVEC(3) C C READ-IN THE NAMES OF THE GEODETIC AND THE GEOCENTRIC ELLIPSOIDS RESPECTIVELY. C READ(1,70)GELLIP,CGELIP 70 FORMAT(2(2X,A3)) C C READ-IN THE PARAMETERS OF THE GEODETIC ELLIPSOID ,NAMELY ; THE SEMI-MAJOR C AXIS (A)-IN METERS & THE RECIPROCAL OF THE FLATTENNING (1/F)-UNITLESS , THE C DATUM TRANSLATIONS OF THIS REFERENCE ELLIPSOID FROM THE GEOCENTRE (XTRANS , C YTRANS , ZTRANS)-IN METERS , AND THE TRANSLATIONS INDICATOR (IDTRAN) WHICH C SHOULD BE PUNCHED AS 1 ;IF ANYONE OF THE 3-TRANSLATION COMPONENTS IS NON-ZERO C C C C C READ(1,20)A,FINV,XTRANS,YTRANS,ZTRANS,IDTRAN 20 FORMAT(1X,F10.2,2X,F10.6,3(2X,F7.2),4X,I1) C C COMPUTATION OF THE OTHER PARAMETERS OF THE GEODETIC ELLIPSOID (F , B & E) C F=1./FINV B=A*(1.-F) E=DSQRT(2.*F-F*F) IF(CGELIP.EQ.ELNAME) GO TO 71 C C READ-IN THE SEMI-MAJOR AXIS AND THE RECIPROCAL OF FLATTENNING ,FOR THE C ARBITRARILY CHOSEN GEOCENTRIC ELLIPSOID ( IF ANY ). C READ(1,72)ACG,FINVCG 72 FORMAT(1X,F10.2,2X,F10.6) C C COMPUTE THE OTHER PARAMETERS OF THE GEOCENTRIC ELLIPSOID (IF ANY) C FCG=1./FINVCG BCG=ACG*(1.-FCG) ECG=DSQRT(2.*FCG-FCG*FCG) 71 CONTINUE C C READ-IN THE TRACKING STATION NAME AND NUMBER C READ(1,10)NAME,ISTNUM 10 FORMAT(1X,A4,2X,I2) IF(IDTRAN.NE.1) GO TO 1000 C C COMPUTE THE GEODETIC (XG,YG & ZG) CARTESIAN COORDINATES OF THE TRACKING C STATION BY APPLYING THE DATUM TRANSLATIONS TO THE GEOCENTRIC (X,Y,Z) COORD. C XG=X-XTRANS YG=Y-YTRANS ZG=Z-ZTRANS C C CALLING THE SUBROUTINE 'ABOESA' TO CONVERT THE GEOCENTRIC (X,Y,Z) TO THE C CORRESPONDING GEOCENTRIC (LAT,LONG,HT) ,USING THE PARAMETERS OF THE C GEODETIC ELLIPSOID ;I.E. THE GEODETIC ELLIPSOID IS ASSUMED AS GEOCENTERED. C 1000 CALL ABOESA(A,FINV,X,Y,Z,PHI,ALAMBA,H,ITER) IF(ALAMBA.GT.PI) ALAMBA=ALAMBA-2.*PI C C CALLING THE SUBROUTINE 'SIMSIM' TO CONVERT THE GEOCENTRIC (LAT,LONG) FROM C RADIANS TO DEG,MIN & SEC. C CALL SIMSIM(PHI,IDEG1,MIN1,SEC1) CALL SIMSIM(ALAMBA,IDEG2,MIN2,SEC2) IF(IDTRAN.NE.1) GO TO 2000 C C CALLING THE SUBROUTINE ABOESA TO TRANSFORM THE GEODETIC (X,Y,Z) TO ITS C CORRESPONDING GEODETIC (LAT,LONG,HT) ,W.R.T. THE GEODETIC ELLIPSOID. C CALL ABOESA(A,FINV,XG,YG,ZG,PHI,ALAMBA,HG,ITERG) IF(ALAMBA.GT.PI) ALAMBA=ALAMBA-2.*PI C C C C CALLING THE SUBROUTINE 'SIMSIM' TO CONVERT THE GEODETIC (LAT,LONG) FROM C RADIANS TO DEG,MIN & SEC. C CALL SIMSIM(PHI,IDEG3,MIN3,SEC3) CALL SIMSIM(ALAMBA,IDEG4,MIN4,SEC4) 2000 IF(CGELIP.EQ.ELNAME) GO TO 75 C C CALLING THE SUBROUTINE 'ABOESA' TO TRANSFORM THE GEOCENTRIC (X,Y,Z) TO ITS C CORRESPONDING GEOCENTRIC (LAT,LONG,HT) ,W.R.T. THE GEOCENTRIC ELLIPSOID. C CALL ABOESA(ACG,FINVCG,X,Y,Z,PHI,ALAMBA,HCG,ITERCG) IF(ALAMBA.GT.PI) ALAMBA=ALAMBA-2.*PI C C CALLING THE SUBROUTINE 'SIMSIM' TO CONVERT THE GEOCENTRIC LAT & LONG (W.R.T. C THE GEOCENTRIC ELLIPSOID) ,FROM RADIANS TO ; DEG , MIN & SEC. C CALL SIMSIM(PHI,IDEG5,MIN5,SEC5) CALL SIMSIM(ALAMBA,IDEG6,MIN6,SEC6) C 75 RETURN C C******************************************************************************* C C *** PRINT-OUT THE SUMMARY OF THE INPUT AND OUTPUT RESULTS *** C ENTRY PRINTI C WRITE(3,591) 591 FORMAT('1') WRITE(3,592) 592 FORMAT(//) C C PRINT THE NAME ,PARAMETERS ,AND TRANSLATIONS OF THE GEODETIC ELLIPSOID. C WRITE(3,5003) 5003 FORMAT(10X,'PARAMETERS OF THE GEODETIC REFERENCE ELLIPSOID :'/10X, 1'----------------------------------------------') WRITE(3,30)GELLIP,A,B,FINV,E 30 FORMAT(62X,'NAME = ',A3/65X,'A = ',F10.2,3X,'METERS'/65X,'B = ',F1 10.2,3X,'METERS'/65X,'F = 1/',F10.6/65X,'E = ',F10.8) WRITE(3,592) WRITE(3,40)XTRANS,YTRANS,ZTRANS 40 FORMAT(10X,'GEODETIC DATUM TRANSLATIONS FROM THE GEOCENTRE :'/10X, 1'----------------------------------------------'/60X,'XTRANS = ',F 210.2,3X,'METERS'/60X,'YTRANS = ',F10.2,3X,'METERS'/60X,'ZTRANS = ' 3,F10.2,3X,'METERS') WRITE(3,592) IF(CGELIP.EQ.ELNAME) GO TO 76 C C PRINT THE NAME AND PARAMETERS OF THE CHOSEN GEOCENTRIC ELLIPSOID (IF ANY). C WRITE(3,5004) 5004 FORMAT(10X,'PARAMETERS OF THE CHOZEN GEOCENTRIC ELLIPSOID :'/10X,' 1---------------------------------------------') WRITE(3,30)CGELIP,ACG,BCG,FINVCG,ECG C C C WRITE(3,592) 76 CONTINUE C C PRINT THE FINAL ADJUSTED VALUES OF THE GEOCENTRIC (AND GEODETIC-IF DESIRED) C CARTESIAN COORDINATES OF THE TRACKING STATION. C WRITE(3,600) 600 FORMAT(//17X,'TABLE OF FINAL ADJUSTED CARTESIAN COORDINATES OF THE 1 TRACKING STATION'/17X,'========================================== 2==========================='/) WRITE(3,700) 700 FORMAT(5X,'STATION NAME',2X,'ST. NO.',6X,'X -(METERS)',6X,'Y -(MET 1ERS)',6X,'Z -(METERS)',10X,'DESCRIPTION'/) WRITE(3,800)NAME,ISTNUM,X,Y,Z 800 FORMAT(9X,A4,8X,I2,3X,3(5X,F12.3)) WRITE(3,900) 900 FORMAT('+',84X,'***GEOCENTRIC***') IF(IDTRAN.NE.1) GO TO 3000 WRITE(3,800)NAME,ISTNUM,XG,YG,ZG WRITE(3,860)GELLIP 860 FORMAT('+',82X,'*** GEODETIC-',A3,' ***') C C PRINT THE FINAL ADJUSTED VALUES OF THE GEOCENTRIC (AND GEODETIC-IF DESIRED) C GEOGRAPHICAL COORDINATES (LAT,LONG,HT) OF THE TRACKING STATION C 3000 WRITE(3,592) WRITE(3,200) 200 FORMAT(//15X,'TABLE OF FINAL ADJUSTED GEOGRAPHICAL COORDINATES OF 1THE TRACKING STATION'/15X,'======================================= 2================================='/) WRITE(3,300) 300 FORMAT(3X,'STATION',2X,'STATION',5X,'LATITUDE',10X,'LONGITUDE',11X 1,'HEIGHT',4X,'NUMBER OF ',5X,'DESCRIPTION'/5X,'NAME',4X,'NUMBER',3 2X,'(+VE NORTH)',9X,'(+VE EAST)',9X,'(METERS)',3X,'ITERATIONS'/21X, 3'DEG',2X,'MIN',2X,'SEC',6X,'DEG',2X,'MIN',2X,'SEC') WRITE(3,400)NAME,ISTNUM,IDEG1,MIN1,SEC1,IDEG2,MIN2,SEC2,H,ITER 400 FORMAT(5X,A4,5X,I2,5X,I3,2X,I2,2X,F7.4,3X,I3,2X,I2,2X,F7.4,3X,F10. 15,5X,I2) WRITE(3,500)GELLIP 500 FORMAT('+',82X,'***GEOCENTRIC-',A3,' ***') IF(IDTRAN.NE.1) GO TO 4000 WRITE(3,400)NAME,ISTNUM,IDEG3,MIN3,SEC3,IDEG4,MIN4,SEC4,HG,ITERG WRITE(3,860)GELLIP 4000 IF(CGELIP.EQ.ELNAME) GO TO 81 C C PRINT THE GEOCENTRIC (LAT,LONG,HT) W.R.T. THE GEOCENTRIC ELLIPSOID. C WRITE(3,400)NAME,ISTNUM,IDEG5,MIN5,SEC5,IDEG6,MIN6,SEC6,HCG,ITERCG WRITE(3,500)CGELIP 81 RETURN END C C C C C******************************************************************************* C* * C* * C* SUBROUTINE 'ABOESA' : ( M. NASSAR - DEC. 1971 ) * C* =================== * C* THE PURPOSE OF THIS SUBROUTINE IS TO SOLVE THE INVERSE * C* CASE OF TRANSFORMATION BETWEEN THE GEODETIC AND THE RECTANGULAR COORDINATES * C* WITH RESPECT TO A GIVEN REFERENCE ELLIPSOID. ,NAMELY COMPUTING (LATITUDE , * C* LONGITUDE & ELLIPSOIDAL HEIGHT) ; FROM (X , Y & Z) CARTESIAN COORDINATES. * C* * C* ATTENTION : THIS SUBROUTINE IS VALID ONLY FOR ONE TRANSFORMATION PER EACH * C* --------- CALLING STATEMENT. * C* * C* THE REQUIRED INPUT DATA ARE : * C* --------------------------- * C* 1-THE PARAMETERS OF THE REFERENCE ELLIPSOID ; NAMELY THE SEMI-MAJOR AXIS * C* (A)-IN METERS ,AND THE RECIPROCAL OF THE FLATTENNING (1/F)-UNITLESS * C* 2-THE CARTESIAN RECTANGULAR COORDINATES (X , Y & Z)-IN METERS. * C* * C* THE OUTPUT WILL BE : * C* ------------------ * C* 1-THE GEODETIC LATITUDE ; IN RADIANS. * C* 2-THE GEODETIC LONGITUDE ; IN RADIANS. * C* 3-THE GEOMETRICAL HEIGHT ABOVE THE REFERENCE ELLIPSOID ; IN METERS. * C* * C* THIS SUBROUTINT USES ANOTHER SUBROUTINE CALLED 'MADAD' TO COMPUTE LONGITUDE.* C* * C******************************************************************************* C SUBROUTINE ABOESA(A,FI,X,Y,Z,PHI,ALAMBA,H,J) IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION DATAN,DSIN,DCOS,DSQRT,DABS PI=3.141592653589793 RHO=180.*3600./PI C C COMPUTATION OF THE OTHER PARAMETERS OF THE REFERENCE ELLIPSOID. C F=1/FI B=A*(1.-F) E=DSQRT(2.*F-F*F) C C COMPUTATION OF THE GEODETIC LONGITUDE FROM ; TAN (LAMBDA) = Y/X C BY CALLING THE FORTRAN-SUBROUTINE 'MADAD' C CALL MADAD(Y,X,ALAMBA) C C COMPUTATION OF THE GEODETIC LATITUDE AND ELLIPSOIDAL HEIGHT ; ( ITERATIVELY ) C J=1 P=DSQRT(X*X+Y*Y) C C COMPUTE AN APPROXIMATE INITIAL VALUE FOR THE LATITUDE (PHIA) FROM C WHICH WE CAN GET VALUES FOR N & H C C C C PHIA=DATAN(Z/(P*(1.-E*E))) C C START THE ITERATION PROCESS TO GET PHI & H C 20 RCNA=A*A/DSQRT(A*A*DCOS(PHIA)**2+B*B*DSIN(PHIA)**2) HA=P/DCOS(PHIA)-RCNA C C COMPUTE THE LATEST VALUE OF PHI USING THE RIGOROUS FORMULA ,AND TESTING C THE DIFFERENCE BETWEEN THIS VALUE (PHIG) AND THE LATEST APPROXIMATE C ONE (PHIA) AGAINST THE PRECISION LIMIT OF 0.0001 ARCSECOND TO STOP THE C ITERATION PROCESS C PHIG=DATAN(Z*(RCNA+HA)/(P*(RCNA*(1.-E*E)+HA))) IF(DABS(PHIG-PHIA).LT.(0.0001/RHO)) GO TO 30 C C UP-DATING THE APPROX. VALUE OF PHI (PHIA) TO PERFORM ANOTHER ITERATION C PHIA=PHIG J=J+1 GO TO 20 C C TRANSFERE THE FINAL ACCEPTED VALUES OF (PHI & H) AFTER STOPPING THE C ITERATION ,BACK TO THE CALLING ROUTINE C 30 PHI=PHIG H=HA RETURN END C C C C C******************************************************************************* C* * C* SUBROUTINE 'MADAD' : ( M. NASSAR - DEC. 1971 ) * C* -------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO DETERMINE THE * C* ANGLE (IN-RADIANS) FROM ITS TANGENT WHICH IS DEFINED BY THE OPPOSITE AND * C* THE ADJACENT SIDES OF THE REQUIRED ANGLE. MOREOVER THIS SUBROUTINE DETERMINE* C* THE CORRECT QUADRANT OF THIS ANGLE. * C* * C* THE REQUIRED INPUT ARE THE OPPOSITE & ADJACENT SIDES OF THE ANGLE (Y & X) * C* RESPECTIVELY ,BOTH HAVE TO BE IN THE SAME UNITS (E.G. METERS OR KILOMETERS) * C* * C* THE OUTPUT WILL BE THE ANGLE (IN-RADIANS) IN THE CORRECT QUADRANT ; I.E. ITS* C* VALUE WILL LIE BETWEEN 0 AND 2*PI ,DEPENDING ON THE QUADRANT NUMBER. * C* * C* NOTE : THIS SUBROUTINE IS VALID ONLY FOR COMPUTING ONE ANGLE PER EACH * C* ---- CALLING STATEMENT. * C* * C******************************************************************************* C SUBROUTINE MADAD(Y,X,F) IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION DATAN,DABS PI=3.141592653589793 IF(X)27,25,26 25 IF(Y)11,12,13 11 F=3.*PI/2. GO TO 98 12 PRINT33 33 FORMAT(//20X,'CHECK YOUR INPUT DATA') GO TO 98 13 F=PI/2. GO TO 98 26 IF(Y)14,15,16 14 Z=DABS(Y/X) F=2.*PI-DATAN(Z) GO TO 98 15 F=0. GO TO 98 16 F=DATAN(Y/X) GO TO 98 27 IF(Y)17,18,19 17 F=PI+DATAN(Y/X) GO TO 98 18 F=PI GO TO 98 19 Z=DABS(Y/X) F=PI-DATAN(Z) 98 RETURN END C C C C C******************************************************************************* C* * C* * C* SUBROUTINE 'SIMSIM' : ( M. NASSAR - DEC. 1971 ) * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO CONVERT THE * C* GIVEN ANGLE FROM RADIANS TO : DEGREES ,MINUTES & SECONDS * C* * C* NOTE HERE THAT ,IF THE GIVEN ANGLE IS -VE ,THE DEGREES ONLY WILL BE -VE * C* AFTER TRANSFORMATION ,BUT THE MINUTES & SECONDS WILL BE +VE : THIS HAS BEEN * C* DESIGNED TO GET A BETTER APPEARANCE OF THE PRINTOUT BEING DONE THROUGH THE * C* CALLING ROUTINE. AND THE -VE SIGN IN FRONT OF THE DEGREES WILL INDICATE * C* THAT THE WHOLE ANGLE IS -VE. * C* * C* NOTE : THIS SUBROUTINE CAN CONVERT ONLY ONE ANGLE PER EACH CALLING STATEMENT* C* * C******************************************************************************* C SUBROUTINE SIMSIM(AB,ID,M,S) IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION DABS PI=3.141592653589793 RHO=180.*3600./PI SIGN=AB AB=DABS(AB) S=AB*RHO+0.005 ID=S/3600. D=ID M=S/60.-D*60. RM=M IF(SIGN.LT.0.) ID=-ID S=S-RM*60.-D*3600.-0.005 RM=DABS(RM) M=RM S=DABS(S) RETURN END C C C C******************************************************************************* C* * C* SUBROUTINE 'EIGEN3' : * C* ------------------- * C* * C* (REFERENCE ; 'APPLIED NUMERICAL METHODS' BY CARNAHAN * C* PAGES 250-260) * C* THE PURPOSE OF THIS SUBROUTINT IS TO COMPUTE THE EIGENVALUES AND THE * C* EIGENVECTORS OF A (3,3) REAL SYMMETRIC MATRIX-A ,USING 'JACOBI-METHOD' , * C* WHOSE ALGORITHM ITERATIVELY REDUCES THE INPUTTED MATRIX-A TO A NEARLY * C* DIAGONAL MATRIX-B ; BY SUCCESSIVELY ANNIHILATING OFF-DIAGONAL ELEMENTS * C* THROUGH A SERIES OF ROTATIONS. * C* * C* THE FINAL PRODUCT OF THESE ROTATIONS WILL BE THE ORTHOGONAL MATRIX-R * C* WHOSE COLUMNS ARE THE THREE EIGENVECTORS OF THE ORIGINAL INPUTTED MATRIX-A * C* * C* THE DIAGONAL ELEMENTS OF THE RESULTING DIAGONAL MATRIX-B ,WILL BE THE * C* EIGENVALUES OF THE MATRIX-A ,WHICH IN TURN ARE THE SQUARES OF THE SEMI- * C* AXES (MAJOR , INTERMEDIATE & MINOR) OF THE STANDARD ERROR ELLIPSOID. * C* * C* THE INPUT IS : THE (3,3) REAL SYMMETRIC MATRIX-A * C* THE OUTPUT IS : THE (3,3) DIAGONAL MATRIX-B * C* * C******************************************************************************* C SUBROUTINE EIGEN3(A,B) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3,3),B(3,3),R(3,3),Y(3) C C SPECIFY THE TESTING AIDS NECESSARY FOR THE ITERATION AND DIAGONALIZATION C PROCESSES ,NAMELY ; C ITMAX : IS THE MAXIMUM NUMBER OF ITERATIONS TO BE DONE C EPS1 : IS A TEST FOR THE ROTATION ANGLE (IF ALL DIAGONAL ELEMENTS ARE C EQUAL ,ROTATE BY 45 DEGREES) C EPS2 : IS THE PRECISION LIMIT BELOW WHICH NEGLECT THE DIAGONAL ELEMENTS C EPS3 : IS THE PRECISION LIMIT BELOW WHICH THE ITERATION PROCESS IS STOPPED C ,WHEN THE SUM OF SQUARES OF THE DIAGONAL ELEMENTS WILL BE THE SAME C FOR TWO CONSECUTIVE ITERATIONS. C ITMAX=50 EPS1=0.1D-12 EPS2=0.1D-12 EPS3=0.1D-10 N=3 NM1=N-1 C C SET MATRIX-R EQUALS IDENTITY , AND UPPER TRIANGULATE MATRIX-A C DO 16 I=1,N DO 16 J=1,N R(I,J)=0. IF(I.EQ.J) R(I,J)=1. IF(I.GT.J) A(I,J)=0. 16 CONTINUE C C C C COMPUTE SIGMA1 AND S C SIGMA1=0. OFFDSQ=0. DO 17 I=1,N SIGMA1=SIGMA1+A(I,I)*A(I,I) IP1=I+1 IF(I.GE.N) GOTO 18 DO 17 J=IP1,N 17 OFFDSQ=OFFDSQ+A(I,J)*A(I,J) 18 S=2.0*OFFDSQ+SIGMA1 C C BEGIN THE JACOBI ITERATION C DO 29 ITER=1,ITMAX DO 27 I=1,NM1 IP1=I+1 DO 27 J=IP1,N Q=DABS(A(I,I)-A(J,J)) C C COMPUTE THE SINE AND COSINE OF THE ROTATION ANGLE C IF(Q.LE.EPS1) GOTO 19 IF(DABS(A(I,J)).LE.EPS2) GOTO 27 P=2.0*A(I,J)*Q/(A(I,I)-A(J,J)) SPQ=DSQRT(P*P+Q*Q) CSA=DSQRT((1.0+Q/SPQ)/2.0) SNA=P/(2.0*CSA*SPQ) GOTO 20 19 CSA=1.0/DSQRT(2.0D0) SNA=CSA 20 CONTINUE C C UP-DATE COLUMNS I & J OF MATRIX-R ; EQUIVALENT TO MULTIPLICATION BY THE C ANNIHILATION MATRIX C DO 21 K=1,N HOLDKI=R(K,I) R(K,I)=HOLDKI*CSA+R(K,J)*SNA 21 R(K,J)=HOLDKI*SNA-R(K,J)*CSA C C COMPUTE THE NEW ELEMENTS OF THE MATRIX-A IN ROWS I & J C DO 24 K=1,N IF(K.GT.J) GOTO 23 Y(K)=A(I,K) A(I,K)=CSA*Y(K)+SNA*A(K,J) IF(K.NE.J) GOTO 22 A(J,K)=SNA*Y(K)-CSA*A(J,K) 22 GOTO 24 23 HOLDIK=A(I,K) A(I,K)=CSA*HOLDIK+SNA*A(J,K) A(J,K)=SNA*HOLDIK-CSA*A(J,K) 24 CONTINUE C C COMPUTE THE NEW ELEMENTS OF MATRIX-A IN THE COLUMNS I & J C C C Y(J)=SNA*Y(I)-CSA*Y(J) C C WHEN K IS LARGER THAN R C DO 26 K=1,J IF(K.LE.I) GOTO 25 A(K,J)=SNA*Y(K)-CSA*A(K,J) GOTO 26 25 HOLDKI=A(K,I) A(K,I)=CSA*HOLDKI+SNA*A(K,J) A(K,J)=SNA*HOLDKI-CSA*A(K,J) 26 CONTINUE 27 A(I,J)=0. C C FIND SIGMA2 FOR THE TRANSFORMED-A , AND TEST FOR CONVERGENCE. C SIGMA2=0. DO 28 I=1,N 28 SIGMA2=SIGMA2+A(I,I)*A(I,I) IF(1.0-SIGMA1/SIGMA2.LT.EPS3) GO TO 800 29 SIGMA1=SIGMA2 WRITE(3,500)ITER 500 FORMAT(//5X,'THE EIGEN-VALUE PROBLEM DOES NOT CONVERGE AFTER',I4,2 1X,'ITERATIONS'/) 800 CONTINUE DO 700 I=1,N DO 700 J=1,N 700 B(I,J)=A(I,J) RETURN END