C*************************************************************************** C* * C* * C* *************************** * C* * * * C* * 'D O P S A T S'-PROGRAM * * C* * * * C* *************************** * C* * C* PARAMETRIC LEAST SQUARES ADJUSTMENT * C* ----------------------------------- * C* OF DOPPLER SATELLITE OBSERVATIONS * C* --------------------------------- * C* *RANGE DIFFERENCE MODE**ORBIT KNOWN CASE* * C* ----------------------------------------- * C* *SINGLE TRACKING STATION**MULTI-SATELLITE PASSES* * C* ------------------------------------------------- * C* IN A PHASED (OR SEQUENTIAL)-TYPE SOLUTION ,USING ONE PASS AT A TIME * C* ------------------------------------------------------------------- * C* VIA UPDATING AND WEIGHTING OF THE UNKNOWN PARAMETERS AFTER EACH PHASE * C* --------------------------------------------------------------------- * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* UNB - SURVEYING ENGINEERING PROGRAM LIBRARY * C* PROGRAM DOCUMENTATION - PROGRAM NUMBER * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* 'DOPSATS'-PROGRAM IS THE UNB MODIFIED VERSION FROM THE ORIGINAL * C* PROGRAM : 'SIMULTANEOUS SATELLITE NAVIGATION PROGRAM' WHICH WAS * C* DEVELOPED BY :-SHELL CANADA DEVELOPMENT PROGRAM - AREA SERVICES , * C* SHELL CANADA LTD. ,EDMONTON ,ALBERTA - DATED ; DECEMBER 1970. * C* * C* UNB MODIFICATIONS : JULY 1971 D.WELLS * C* NOV. 1971 M.NASSAR * C* DEC. 1971 P.KIRKHAM * C* JUNE 1972 M.NASSAR * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* *NASR C* THE PURPOSE OF THIS PROGRAM IS TO PERFORM THE PHASED-TYPE SOLUTION *NASR C* FOR THE PARAMETRIC LEAST SQUARES ADJUSTMENT OF DOPPLER SATELLITE *NASR C* OBSERVATIONS ,TAKEN FROM ONE CHOSEN MOTIONLESS TRACKING STATION *NASR C* (ICHOOZ) ,AND COLLECTED OVER DIFFERENT PERIODS OF OBSERVING MULTI- *NASR C* SATELLITE-PASSES . *NASR C* *NASR C* THE PROGRAM ACCEPTS MANY DOPPLER DATA FROM STATIONS OTHER THAN *NASR C* ICHOOZ-STATION (UP TO 7 STATIONS) ,FOR DOPPLER TRIMMING PURPOSES *NASR C* WHEN DEALING WITH TRANSLOCATION AND/OR SIMULTANEOUS MODES. *NASR C* *NASR C* THE SOLUTION WILL BE PERFORMED FOR THE THREE-DIMENSIONAL CARTESIAN *NASR C* COORDINATES (X,Y,Z) OF THE ICHOOZ-STATION IN THE AVERAGE TERRESTR- *NASR C* IAL SYSTEM ,ALONG WITH THE ASSOCIATED FREQUENCY OFFSET UNKNOWNS *NASR C* (ONE PER EACH SATELLITE-PASS). *NASR C* *NASR C* THE ALGORITHM BEING USED IN THIS PROGRAM FOR THE PHASED-TYPE ADJUS-*NASR C* TMENT SOLVES FOR ONE PASS ONLY AT A TIME ,UTILIZING THE TECHNIQUES *NASR C* OF WEIGHTING AND UPDATING THE INVOLVED UNKNOWN PARAMETERS AFTER *NASR C* EACH PASS (PHASE) BEING PROCESSED. *NASR C* *NASR C* IN ADDITION THE PROGRAM PERFORMS A CHI-SQUARED STATISTICAL TEST , *NASR C* FOR REJECTING THE POOR PASSES ,WITHOUT CARRYING-OUT THEIR CONTRIB- *NASR C* UTION TO THE SUBSEQUENT PHASES.. *NASR C* *NASR C* FINALLY ,AFTER PROCESSING ALL INPUTTED PASSES ,THE PROGRAM WILL *NASR C* PRINT-OUT A SUMMARY OF THE FINAL RESULTS AND ACCURACY ESTIMATES *NASR C* ,THIS WILL APPEAR AT THE END OF THE WHOLE PRINTED OUTPUT. *NASR C* *NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (IN THE FOLLOWING SEQUENCE) *NASR C* --------------------------- ( F R E E -FORMAT IS REQUIRED) *NASR C* *NASR C* I-THE PROGRAM OPTIONS : ITROP,HOR,ITRIM,IPUNCH,SIGMR,TPCORR, *NASR C* INTEIG ,AND IDNT. (ALL ON ONE CARD) *NASR C* *NASR C* II-THE PARAMETERS OF THE REFERENCE ELLIPSOID ;NAMELY : SEMI-MAJOR *NASR C* AXIS IN METERS ,RECIPROCAL OF FLATTENNING (DIMENSIONLESS) ,AND *NASR C* THE THREE TRANSLATION COMPONENTS FROM THE GEOCENTRE -IN METERS *NASR C* (ALL ON ONE CARD) *NASR C* *NASR C*III-THE TOTAL NUMBER OF INPUTTED GROUND STATIONS TO BE PROCESSED IN *NASR C* THIS RUN (UP TO 7) , AND THE PARTICULAR CHOSEN STATION-NUMBER *NASR C* 'ICHOOZ' ,FOR WHICH ONLY THE PHASED-SOLUTION WILL BE PERFORMED *NASR C* (BOTH NUMBERS ON ONE CARD) *NASR C* *NASR C* IV-A NUMBER OF CARDS EQUAL TO THE TOTAL NUMBER OF GROUND STATIONS *NASR C* (UP TO 7) ,EACH OF WHICH SHOULD CONTAIN :- STATION NUMBER ,STATN.*NASR C* APPROXIMATE LATITUDE IN (DEG,MIN & SEC) +VE NORTH ,STATION APPRO-*NASR C* XIMATE LONGITUDE IN (DEG,MIN & SEC) +VE EAST , AND THE STATION *NASR C* APPROXIMATE ELLIPSOIDAL HEIGHT -IN METERS. *NASR C* *NASR C* V-THE MAJORITY VOTED DATA FOR ALL COLLECTED PASSES FROM ALL DESIRED*NASR C* STATIONS ,ON PUNCHED CARDS ,AS PREPARED BY THE MAJORITY VOTE *NASR C* PROGRAM ,OR MANUALLY BY THE USER. *NASR C* FOR EACH PASS ,THE FOLLOWING INPUT IS REQUIRED :- *NASR C* 1-THE SATELLITE IDENTIFICATION CARD (DAY # ,PASS # ,SAT # ,AND *NASR C* LOCKON TIME (HOUR & MIN)) ............. REQUIRE ONE CARD.*NASR C* 2-THE UNSCALED ELEVEN FIXED ORBITAL PARAMETERS. ...ON ONE CARD. *NASR C* 3-THE SETS OF UNSCALED VARIABLE ORBITAL PARAMETERS (UP TO 17 SETS*NASR C* ,EACH CONTAINS 3 PARAMETERS CAN BE READ-IN) --REQUIRE 2 CARDS. *NASR C* 4-A DELIMETER 999 ,INDICATING THE END OF SATELLITE ORBITAL *NASR C* PARAMETERS (BOTH THE FIXED AND THE VARIABLE). *NASR C* 5-THREE CARDS PER EACH INPUTTED TRACKING STATION ,CONTAIN : THE *NASR C* STATION NUMBER (FIRST CARD) ,AND THE PAIRS OF DOPPLER DATA *NASR C* (400 MHZ & 150 MHZ) TAKEN AT THE STATION FROM THIS PASS -NOTE *NASR C* THAT THE NUMBER OF THESE DOPPLER PAIRS SHOULD BE LESS BY 3 THAN*NASR C* THE NUMBER OF VARIABLE PARAMETERS SETS , IN OTHER WORDS ,IF *NASR C* THERE ARE NO ENOUGH DOPPLERS ,COMPLETE THEIR APPROPRIATE NUMBER*NASR C* BY PUNCHING ZEROS. *NASR C* 6-A DELIMITER 99 ,INDICATING THE END OF DOPPLER DATA FROM ALL *NASR C* INPUTTED STATIONS TO THIS PARTICULAR PASS. *NASR C* *NASR C* VI-AFTER REPEATING -V- FOR EACH PASS ,A DELIMITER 9999 ,IS REQUIRED*NASR C* TO INDICATE THE END OF ALL DATA CARDS INPUTTED FOR THIS RUN. *NASR C* *NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C* *NASR C* THE PRINTED OUTPUT WILL INCLUDE : *NASR C* ------------------------------- *NASR C* *NASR C* I-THE DESIRED OPTIONS FOR THIS RUN. *NASR C* II-THE PARAMETERS OF THE GIVEN REFERENCE DATUM. *NASR C*III-THE APPROXIMATE GEOGRAPHICAL AND APPROXIMATE CARTESIAN COORDINAT-*NASR C* ES OF INPUTTED TRACKING STATIONS. *NASR C* *NASR C* IV-FOR EACH PASS BEING PROCESSED ,THE PROGRAM PRINTS THE FOLLOWING:-*NASR C* 1-THE GENERAL HEADING ,AND THE PASS IDENTIFICATION CARD. *NASR C* 2-THE UNSCALED AND SCALED FIXED ORBITAL PARAMETERS. *NASR C* 3-THE UNSCALED AND SCALED VARIABLE ORBITAL PARAMETERS. *NASR C* 4-THE SATELLITE CARTESIAN COORDINATES (X,Y,Z) IN THE AVERAGE *NASR C* TERRESTFIAL SYSTEM OF COORDINATES (GEOCENTRIC COORDINATES). *NASR C* 5-DOPPLER DATA FROM THE ICHOOZ-STATION ,BEFORE AND AFTER THE *NASR C* IONOSPHERIC AND TROPOSPHERIC REFRACTION CORRECTIONS. *NASR C* 6-DETAILS ABOUT THE REJECTED DOPPLERS UNDER DIFFERENT CRITERIA , *NASR C* AND THE REMAINING USABLE DOPPLERS. *NASR C* *NASR C* V-REGARDING THE PHASED-ADJUSTMENT RESULTS AFTER PROCESSING EACH *NASR C* PASS ,THE PROGRAM PRINTS THE FOLLOWING :- *NASR C* 1-THE DESIGN MATRIX-A ,AND THE MISCLOSURE VECTOR-W. *NASR C* 2-THE RELATIVE VARIANCE-COVARIANCE MATRIX (IN CASE OF CORRELATED *NASR C* WEIGHTING ONLY) ,AND THE WEIGHT MATRIX-P OF THE USED DOPPLERS. *NASR C* 3-THE INVERSE OF THE NORMAL EQUATIONS MATRIX (WEIGHT COEFFICIENT *NASR C* MATRIX OF THE ESTIMATED PARAMETERS) ,AND THE ABSOLUTE TERMS OF *NASR C* THE NORMAL EQUATIONS (U-VECTOR).. *NASR C* 4-THE SOLUTION VECTOR-X. *NASR C* 5-THE ORIGINAL APPROXIMATE AND THE LATEST ESTIMATED COORDINATES *NASR C* OF THE ICHOOZ-STATION ,ALONG WITH THE DIFFERENCES BETWEEN THEM.*NASR C* 6-THE RESIDUAL VECTOR-V ,IN THE UNITS OF DOPPLER COUNTS (CYCLES)*NASR C* 7-THE ACCUMULATED SQUARE SUM OF THE WEIGHTED RESIDUALS V'PV , *NASR C* COMPUTED BY TWO DIFFERENT METHODS. ALSO THE CONTRIBUTION OF *NASR C* THE CURRENT PASS ONLY TO V'PV ,WILL BE PRINTED-OUT. *NASR C* 8-THE ACCUMULATED NUMBER OF DEGREES OF FREEDOM ,THE ESTIMATED *NASR C* VARIANCE-FACTOR ,AND THE STANDARD DEVIATION OF UNIT WEIGHT - *NASR C* (IN CYCLES). *NASR C* 9-THE ESTIMATED VARIANCE-COVARIANCE MATRIX OF THE ESTIMATED PARA-*NASR C* METERS , AND THE STANDARD DEVIATIONS OF THE (X,Y,Z) COORDINATES*NASR C* OF THE ICHOOZ-STATION -(IN METERS). *NASR C* 10-IF ASKED FOR ; THE THREE SEMI-AXES OF THE STANDARD ERROR ELLIP-*NASR C* SOID FOR THE ESTIMATED (X,Y,Z) COORDINATES WILL BE COMPUTED *NASR C* AND PRINTED IN METERS. *NASR C* *NASR C* VI-AT THE END OF THE PRINTED OUTPUT FOR ALL INPUTTED PASSES ,THE *NASR C* PROGRAM PRINTS THE SUMMARY OF THE FINAL RESULTS ,WHICH WILL BE :-*NASR C* 1-GENERAL AND SPECIFIC HEADINGS. AND SOME DETAILS ABOUT THE ACTU-*NASR C* AL NUMBER OF THE ACCEPTED PASSES. *NASR C* 2-PARAMETERS OF THE REFERENCE ELLIPSOID (GIVEN GEODETIC DATUM) *NASR C* 3-FINAL ADJUSTED COORDINATES OF THE ICHOOZ-STATION ,BOTH GEOGRA- *NASR C* PHICAL AND CARTESIAN COORDINATES W.R.T. THE GIVEN REFERENCE *NASR C* DATUM WILL BE GIVEN. *NASR C* 4-ACCUMULATED NUMBER OF DEGREES OF FREEDOM. AND THE STANDARD DEV-*NASR C* IATION OF UNIT WEIGHT. *NASR C* 5-VARIANCE-COVARIANCE MATRIX AND STANDARD DEVIATIONS OF THE *NASR C* ESTIMATED (X,Y,Z) COORDINATES OF THE ICHOOZ-STATION. *NASR C* 6-VARIANCE-COVARIANCE MATRIX AND STANDARD DEVIATIONS OF THE CORR-*NASR C* ESPONDING (LAT,LONG & HT) OF THE ICHOOZ-STATION-IN METER UNITS.*NASR C* 7-THE THREE SEMI-AXES OF THE STANDARD ERROR ELLIPSOID FOR THE *NASR C* FINAL ESTIMATED COORDINATES ,AFTER THE LAST PASS -IN METERS. *NASR C* *NASR C* THE PROGRAM WILL PUNCH-OUT THE ESTIMATED COORDINATES OF THE ICHOOZ *NASR C* -STATION AND THEIR ESTIMATED VARIANCE-COVARIANCE MATRIX ,AFTER *NASR C* EVERY SPECIFIED NUMBER OF PASSES ; IF ASKED FOR ,OTHERWISE NOPUNCH *NASR C* *NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C* *NASR C* NOTE :- THERE IS NO LIMITATION ON THE NUMBER OF INPUTTED PASSES *NASR C* TO BE PROCESSED BY THIS PROGRAM. HOWEVER ,IT IS RESTRIC-*NASR C* TED ONLY FOR A SINGLE MOTIONLESS UNKNOWN TRACKING STATION*NASR C* TO BE PHASED-ADJUSTED. *NASR C* *NASR C***********************************************************************NASR C NASR IMPLICIT REAL*8 (A-H,O-Z) 0300 DIMENSION X(7),Y(7),Z(7),NDESC(18),NG(7),XS(9),YS(9),ZS(9),DETA(17 0310 1),DAK(17),DEK(17),DOP(7,14),COUNT(14),KDPSI(7,14),ROT(3,3),CR(3 0320 2),LR(3),NTEST(7),ALAT(7),ALOG(7),ELMAX(7),SK(7,9),XN(7),YN(7),ZN(7 0330 3),FDI(7),C(8,24),A(8,4),W(8),U(4),PX(4,4),ANORM(4,4),VAR(6,6),AROT 0340 4(6,6),DAR(6,6),XX(04),YY(24),AX(3),IDENT(07),P(3),Q(3) NASR0350 DIMENSION EX(3,3),EVECT(3,3),STORE(50,3),ORGAPR(7,3) NASR DIMENSION PT(9,3),PTOP(9,3),QT(3),LRT(3),XWT(4),V(8) NASR0370 REAL DIR(7),SPEED(7),ELEV(9),DIRCT,CHS(18),VA(9) 0390 REAL*8 LO 0410 INTEGER R(3) 0380 C NASR C SPECIFY A LABELLED COMMON BLOCK /DOPADJ/ ,FOR THE INPUT DATA TO ,AND NASR C THE OUTPUT DATA FROM THE SUBROUTINE 'ADJUST' ,TO SHARE THIS AREA NASR C WITH THE MAIN PROGRAM 'DOPSATS' ; THIS IS TO SIMPLIFY THE CALLING NASR C NASR COMMON /DOPADJ/ D1,ANORM,PX,A,W,XX,U,V,VPV,VPV1,MACOL,MROW, NASR 1 NGRPTS,KDPSI,IDNT NASR C NASR C***********************************************************************NASR C NASR C SETTING-UP A GENERAL HEADING TO BE PRINTED-OUT ON THE TOP OF EACH NASR C DESIRED OUTPUT PAGE ,WHEN ASKED FOR. NASR C NASR DATA NDESC /'DOPS','ATS ','SIMU','LTAN','EOUS',' FOR', 2940 1 'MAT ','- SI','NGLE',' STA','TION',' SOL','UTIO','N US', 2950 2 'ING ','MULT','I PA','SSES'/ 2960 C NASR C INSERTING THE CHI-SQUARE TABULATED VALUES FOR THE UPPER LIMIT ONLY NASR C ACCORDING TO DIFFERENT DEGREES OF FREEDOM ,WITH 95% PROBABILITY NASR C THESE VALUES WILL BE USED FOR THE STATISTICAL TEST BUILT-IN THE NASR C PROGRAM FOR REJECTING THE POOR PASSES. NASR C NASR DATA CHS/5.02,7.38,9.35,11.1,12.8,14.4,16.,17.5,19.,20.5,21.9,23.3 0420 1,24.7,26.1,27.5,28.8,30.2,31.5/ 0430 C NASR C***********************************************************************NASR C NASR C ASSIGN VALUES FOR THE PROGRAM CONSTANTS AS FOLLOWS : NASR C NASR C INPUT/OUTPUT UNITS NASR C NASR C INPT- INPUT UNIT NUMBER 0510 C IOUT-OUTPUT UNIT NUMBER 0520 C NASR INPT=5 0530 IOUT=6 0540 C NASR C PROGRAM ORGANIZING CODES - FOR READING THE INPUT DATA CARDS. NASR C NASR C EDS : IS THE END OF SETS OF VARIABLE PARAMETERS OV THE SATELLITE NASR0570 C ORBIT NASR0580 C EDP : IS THE END OF THE SATELLITE PASS , IT MEANS THE END OF ALL NASR0590 C DOPPLER DATA FROM ALL OBSERVING STATIONS TO THIS PASS NASR0600 C EDD : IS THE END OF DATA FOR ALL PASSES NASR0610 C NASR EOS=999. 0630 EOP=99. 0640 EOD=9999. 0650 C 0660 C PHYSICAL CONSTANTS 0670 C 0620 PI=3.141592653589793D0 0680 RHO=3600.D0*180.D0/PI 0690 RADG=PI/180.D0 0700 WE=4.3752695D-03 0710 XVEL=299792500.D0 0720 C 0730 C NNSS CONSTANTS 0740 C 0550 XFREQN=400.D6 0750 LO=XVEL/XFREQN 0760 WAVENO=1.0D0/LO 0770 FO=1.920000D+06 0780 C NASR C TROPOSPHERIC REFRACTION CORRECTION CONSTANTS - ACCORDING TO HOPFIELD NASR C MODEL (HOPFIELD 1969). NASR C NASR XKD=2.31D0 0820 XKW=0.20D0 0830 THD=(2.5*PI/180.)**2 0840 THW=(1.5*PI/180.)**2 0850 C NASR C***********************************************************************NASR C NASR C INITIALIZING VALUES FOR THE ARRAY INDICIES. NASR C NASR MROW=8 0900 MACOL=4 0910 MCCOL=24 0920 C NASR C INITIALIZING ZERO VALUES FOR SOME SCALARS. NASR C NASR ISEQ=0 0960 NPTS=0 0970 NGRPTS=0 0980 VPV1=0. 0990 VPV=0. 1000 MDFM=0 1010 C NASR C INIALIZING ZERO VALUES FOR SOME ARRAYS. NASR C NASR DO 100 IA=1,7 1050 FDI(IA)=FO 1060 X(IA)=0. 1070 Y(IA)=0. 1080 Z(IA)=0. 1090 NTEST(IA)=0 1100 ALAT(IA)=0. 1110 ALOG(IA)=0. 1120 NG(IA)=0 1130 SPEED(IA)=0. 1140 DIR(IA)=0. 1150 ELMAX(IA)=0. 1160 ELEV(IA+2)=0. 1170 DO 6400 IB=1,9 1180 6400 SK(IA,IB)=0. 1190 DO 100 IB=1,14 1210 100 DOP(IA,IB)=0. 1220 DO 6331 IPZ=1,14 1230 6331 KDPSI(1,IPZ)=0 1240 DO 590 IB=1,4 1250 DO 590 IA=1,4 1260 590 PX(IA,IB)=0. 1270 C NASR C***********************************************************************NASR C NASR C READ-IN THE PROGRAM OPTIONS FOR THIS RUN. NASR C ALL WILL BE ON ONE CARD (FREE FORMAT) IN THE FOLLOWING SEQUENCE : NASR C ITROP : IS A SWITCH FOR THE TROPOSPHERIC CORRECTION ; NASR C ITROP = 1 ,MEANS APPLY THE CORRECTION. NASR C HOR : IS THE MAXIMUM ELEVATION ANGLE (IN DEGREES) ,BELOW WHICH NASR C REJECT OBSERVED DOPPLERS AT THAT LOW ELEVATION ,WHICH IS NASR C KNOWN AS THE HORIZON CUT-OFF CRITERION. NASR C ITRIM : IS A SWITCH FOR TRIMMING DOPPLERS FROM ALL PROCESSSED STNS. NASR C TO KEEP ONLY THE SIMULTANEOUS DOPPLERS ,WHICH ARE ESSENTIAL NASR C FOR THE TRANSLOCATION AND/OR THE SIMULTANEOUS MODES . NASR C ITRIM = 1 ,MEANS TRIM DOPPLERS. NASR C IPUNCH: IS A SWITCH FOR PUNCHING-OUT THE ESTIMATED COORDINATES OF NASR C THE TRACKING STATION ,ALONG WITH THEIR VARIANCE-COV. MATRIX NASR C AFTER SPECIFIED NUMBER OF PASSES (FOR FURTHER INVESTIGATIONS)NASR C IPUNCH=1 ,PUNCH WITH FORMAT : 3E13.5 NASR C IPUNCH=2 ,PUNCH WITH FORMAT : 39X,3E13.5 NASR C OTHERWISE ,IPUNCH SHOULD BE INPUTTED AS 0 (ZERO) NASR C SIGMR : IS THE APRIORI VARIANCE FACTOR - IN (CYCLES)-SQUARED. NASR C TPCORR: IS THE NECESSARY CORRECTION FOR THE TIME OF PERIGEE PASSAGE NASR C (IF ANY) - UNITS ARE IN MINUTES OF UNIVERSAL TIME. NASR C INTEIG: IS THE PASS INTERVAL AFTER WHICH THE ERROR ELLIPSOID WILL NASR C BE COMPUTED FOR THE ESTIMATED COORDINATES ,AND AFTER WHICH NASR C THE VAR.-COV. MATRIX WILL BE PUNCHED (IF IPUNCH = 1 OR 2) NASR C IDNT : IS A SWITCH FOR ESTABLISHING THE REQUIRED WEIGHT COEFFICIENT NASR C MATRIX OF THE USED DOPPLERS ; EITHER IDENTITY OR WITH OFF- NASR C DIAGONALS (I.E. WITH CORRELATION) , NASR C IDNT = 1 ,MEANS USE IDENTITY WEIGHT-COEFF. MATRIX NASR ITROP=FREAD(1,KX) 1430 HOR=FREAD(1,KX) 1440 ITRIM=FREAD(1,KX) 1450 IPUNCH=FREAD(1,KX) 1460 SIGMR=FREAD(1,KX) 1470 TPCORR=FREAD(1,KK) NASR INTEIG=FREAD(1,KX) 1490 AINT=INTEIG 1500 IDNT=FREAD(1,KX) 1510 C NASR C PRINT-OUT THE DESIRED OPTIONS FOR THIS RUN. NASR C NASR PRINT 55551,(NDESC(I),I=1,18) NASR PRINT 35046 NASR 35046 FORMAT(/,6X,'THE INPUTTED OPTIONS FOR THIS RUN ARE :'/6X,'--------NASR 1-----------------------------') NASR PRINT 6101,ITROP,HOR,ITRIM,IPUNCH 1550 6101 FORMAT(' ',/10X,'TROP SWITCH= ', NASR1560 1 I6,5X,'(1=APPLY TROPO CORRECTION)',//10X,'HORIZON CUTOFF= ', NASR1570 2 F5.1,3X,'DEGREES'//,10X,'TRIM SWITCH= ',I6,5X,'(1=TRIM DOPS)', NASR1580 3 //10X,'PUNCH SWITCH= ',I5,5X,'(1=PUNCH COVARIANCE MATRIX ACCORDI 1590 4NG TO FORMAT (3E13.5);',/,35X,'2= ACCORDING TO FORMAT (39X,3E13.5)NASR 5; 1 OR 2= PUNCH STN XYZ)',/) 1610 PRINT87654,SIGMR NASR 87654 FORMAT(10X,'A PRIORI VARIANCE FACTOR =',F10.5,3X,'(CYCLES)-SQUAREDNASR 1') NASR PRINT6102,TPCORR,INTEIG,IDNT NASR 6102 FORMAT(/10X,'TP CORRECTION =',D11.2,2X,'(MINUTES OF UNIVERSAL TIMENASR 1)'/' '/10X,'PASS INTERVAL=',I5,5X,'(EIGEN-VALUESNASR 2 COMPUTED AT THIS INTERVAL;',/,35X,'VAR-COV AND STN XYZ PUNCHED AT 1650 3 THIS INTERVAL)'//10X,'IDNT=',I14,5X,'(1=USE AN IDENTITY AS WEIGHT 1660 4 MATRIX;'/35X,'0=USE A CORRELATED WEIGHT MATRIX)') NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE PARAMETERS OF THE REFERENCE DATUM ,NAMELY :- NASR C SEMI-MAJOR AXIS IN METERS , RECIPROCAL OF FLATTENNING (1/F) DIMENSIO-NASR C NLESS , AND THE TRANSLATION COMPONENTS FROM GEOCENTRE IN METERS. NASR C NASR A1=FREAD(1,KX) 1720 F=FREAD(1,KX) 1730 FLAT=F NASR XTRANS=FREAD(1,KX) 1740 YTRANS=FREAD(1,KX) 1750 ZTRANS=FREAD(1,KX) 1760 C NASR C COMPUTE THE SEMI-MINOR AXIS AND THE SQUARE OF FIRST ECCENTRICITY NASR C OF THE REFERENCE ELLIPSOID. NASR C NASR AM=A1*(F-1.)/F 1780 ESQ=1./F*(2.-1./F) 1790 C NASR C PRINT-OUT THE PARAMETERS AND TRANSLATIONS FROM GEOCENTRE OF THE NASR C GIVEN REFERENCE ELLIPSOID. NASR C NASR PRINT11223 NASR 11223 FORMAT(//10X,'PARAMETERS OF THE REFERENCE DATUM ARE :'/10X,'------NASR 1-------------------------------') NASR PRINT32211,A1,AM,F,XTRANS,YTRANS,ZTRANS NASR 32211 FORMAT(65X,'A = ',F10.2,3X,'METERS'/65X,'B = ',F10.2,3X,'METERS'/,NASR 165X,'F = 1/',F10.5/,60X,'XTRANS = ',F10.2,3X,'METERS'/,60X,'YTRANSNASR 2 = ',F10.2,3X,'METERS'/,60X,'ZTRANS = ',F10.2,3X,'METERS'/) NASR C NASR C***********************************************************************NASR C NASR C READ-IN THE TOTAL NUMBER OF TRACKING STATIONS TO BE PROCESSED IN THISNASR C RUN ,WHICH SIMULTANEOUSLY OBSERVED THE INPUTTED SATELLITE PASSES ; NASR C ( UP TO 7 STATIONS ARE ACCEPTED ) NASR C NASR NSTA=FREAD(1,KX) 1860 C NASR C READ-IN THE STATION-NUMBER OF THE PARTICULAR CHOSEN STATION 'ICHOOZ' NASR C ,FOR ONLY WHICH THE PHASED-ADJUSTMENT AND SOLUTION WILL BE PERFORMED.NASR C NASR ICHOOZ=FREAD(1,KX) 1870 C NASR C READ NUMBER AND APPROXIMATE COORDINATES (LAT.,LONG.&HT.) FOR EATCH NASR1900 C TRACKING STATION ,COMPUTE THE APPROXIMATE (X,Y&Z) ,AND PRINT-OUT THE NASR C BOTH SETS OF APPROXIMATE COORDINATES. NASR C NASR PRINT5001 NASR1940 5001 FORMAT(/,/,/5X,'GIVEN APPROXIMATE COORDINATES OF THE OBSERVING') NASR1950 PRINT 6007 1960 6007 FORMAT(' ',4X,'GROUND STATIONS',/,/,4X,'NO',8X,'LAT',9X, NASR 1 'LONG',8X,'HEIGHT',11X,'X',15X,'Y',15X,'Z',5X,'(IN METERS)'/) NASR IDENT(1)=ICHOOZ 1990 K=1 2000 DO 6300 I=1,NSTA 2010 ISTN=FREAD(1,KX) 2020 IF(ISTN.EQ.ICHOOZ) GOTO 6302 2030 K=K+1 2040 IDENT(K)=ISTN 2050 6302 CONTINUE 2060 M1=FREAD(1,KX) 2070 MPK1=M1 2080 L1=IABS(M1) 2090 N1=FREAD(1,KX) 2100 N1=IABS(N1) 2110 S1=FREAD(1,KX) 2120 S1=DABS(S1) 2130 M2=FREAD(1,KX) 2140 MPK2=M2 2150 L2=IABS(M2) 2160 N2=FREAD(1,KX) 2170 N2=IABS(N2) 2180 S2=FREAD(1,KX) 2190 S2=DABS(S2) 2200 H=FREAD(1,KX) 2210 J=K 2220 IF(ISTN.EQ.ICHOOZ) J=1 2230 ALAT(J)=M1/L1*(L1*3600.D0+N1*60.D0+S1)/RHO 2240 ALOG(J)=M2/L2*(L2*3600.D0+N2*60.D0+S2)/RHO 2250 C NASR C CALLING THE SUBROUTINE 'DGEOCA' TO CONVERT THE APPROXIMATE (LAT,LONG NASR C & HT) OF GIVEN STATIONS TO THEIR CORRESPONDING (X,Y & Z) COORDINATES.NASR C NASR CALL DGEOCA(D,E,F,A1,ESQ,ALAT(J),ALOG(J),H) 2260 C 1930 X(J)=D 2270 Y(J)=E 2290 Z(J)=F 2300 C NASR C PRINT BOTH GEOGRAPHICAL & CARTESIAN COORDINATES OF EACH GIVEN STATN. NASR C NASR PRINT 6006,ISTN,M1,N1,S1,M2,N2,S2,H,D,E,F 2310 6006 FORMAT('0',2X,I4,2X,2(I3,I3,F7.3,1X),1X,F8.2,4X,3(F12.2,3X)) 2320 IF(ISTN.EQ.ICHOOZ) PRINT 6008 2330 6008 FORMAT('+ *') 2340 C NASR2350 C IMPROVING THE INITIAL APPROXIMATION FOR THE RECEIVER GEOCENTRIC NASR2360 C COORDINATES BY APPLYING THE DATUM TRANSLATIONS FROM THE GEOCENTER. NASR2370 C NASR2380 D=D+XTRANS NASR2390 E=E+YTRANS NASR2410 F=F+ZTRANS NASR2420 X(J)=D NASR2430 Y(J)=E NASR2440 Z(J)=F NASR2450 PRINT5002,ISTN,X(J),Y(J),Z(J) NASR2460 5002 FORMAT('0',3X,'ST. NO. =',I4,1X,'APPROXIMATE GEOCENTRIC COORD.',2XNASR24 1,'X=',F11.2,2X,'Y=',F11.2,2X,'Z=',F11.2) NASR2480 IF(ISTN.EQ.ICHOOZ) PRINT 6008 NASR2490 C NASR C RESERVING THE ORIGINAL APPROXIMATE COORDINATES OF THE RECEIVER POSI- NASR C TION FOR FURTHER NEEDS (PRINTING THEM ALONG WITH THEIR CHANGE AFTER NASR C EACH PHASE ADJUSTMENT) NASR C NASR ORGAPR(J,1)=D NASR ORGAPR(J,2)=E NASR ORGAPR(J,3)=F NASR 6300 CONTINUE 2500 PRINT5095 NASR2510 5095 FORMAT('0',15X,'INITIAL APPROXIMATE VALUE FOR THE FREQUENCY OFFSETNASR2520 1 = 32000. CYCLES/SECOND') NASR2530 C NASR C***********************************************************************NASR C NASR C INITIALIZING OF : NPASS ,IWRONG ,IDOPL2 ,INSING ,ICHI95 ,AND KPASS , NASR C WHERE :- NASR C NPASS : IS THE TOTAL NUMBER OF INPUTTED PROCESSED PASSES. NASR C IWRONG : NUMBER OF REJECTED PASSES ,DUE TO WRONG INPUTTED STATION NASR C NUMBERS ,OR DUPLICATED DATA FROM THE SAME CHOSEN STATION. NASR C IDOPL2 : IS THE TOTAL NUMBER OF REJECTED PASSES ,WHICH HAVE LESS NASR C THAN TWO USABLE DOPPLERS. NASR C INSING : NUMBER OF REJECTED PASSES ,DUE TO THE SINGULARITY OF THE NASR C NORMAL EQUATION M ATRIX (E.G. INSUFFICIENT DOPPLER EQUATIONSNASR C ICHI95 : IS THE TOTAL NUMBER OF REJECTED PASSES ,WHICH FAILED THE NASR C CHI-SQUARED TEST AT 95% PROBABILITY. NASR C KPASS : IS THE TOTAL NUMBER OF ACCEPTED PASSES ,AFTER ALL REJECTIONSNASR C NASR NPASS=0 NASR IWRONG=0 NASR IDOPL2=0 NASR INSING=0 NASR ICHI95=0 NASR KPASS=0 NASR C NASR C***********************************************************************NASR C NASR C RETURN HERE TO READ-IN AND PROCESS THE DATA OF ANOTHER SATELLITE PASSNASR C NASR 591 CONTINUE NASR 104 CONTINUE NASR 1047 CONTINUE NASR C NASR DO 12345 IA=1,7 DO 12345 IB=1,14 12345 DOP(IA,IB)=0. C***********************************************************************NASR C NASR C READ-IN AND PRINT-OUT THE IDENTIFICATION CARD OF THE PASS (DAY # , NASR C PASS # , STELLITE # , AND LOCKON TIME -HOUR & MIN-) NASR C 2880 IDAY=FREAD(1,KX) 2890 C 2900 C CHECK FOR END OF DATA 2910 C 2920 IF(IDAY.EQ.EOD)GO TO 9549 2930 NPASS=NPASS+1 NASR IPASS=FREAD(1,KX) 3000 ISAT=FREAD(1,KX) 3010 IH=FREAD(1,KX) 3020 M=FREAD(1,KX) 3030 PRINT 203,(NDESC(I),I=1,18),NPASS 2970 203 FORMAT('1',/,10X,18A4,10X,'PASS NO. ',I5,/) 2980 PRINT 940,IDAY,IPASS,ISAT,IH,M 3040 940 FORMAT(10X,'DAY= ',I5,5X,'PASS= ',I5,5X,'SAT= ',I5,5X,'HOUR= ', 3050 1 I5,5X,'MIN= ',I5/) NASR C NASR C READ-IN AND PRINT-OUT THE UNSCALED ELEVEN FIXED SATELLITE ORBITAL NASR C PARAMETERS ,THEN PERFORM THE NECESSARY SCALING TO THEM. NASR C NASR PRINT5003 NASR3270 5003 FORMAT(5X,'UNSCALED SATELLITE MEAN ORBIT PARAMETERS'/) NASR3280 506 TP=FREAD(2,KX)/1D4 3290 C NASR C CORRECTING THE TIME OF PERIGEE PASSAGE (IF REQUIRED). NASR C NASR TP=TP+TPCORR NASR AN=(FREAD(2,KX)/1D7+3.D0)*RADG 3310 OMEGA=FREAD(2,KX)*RADG/1D4 3320 DOMEGA=FREAD(2,KX)*RADG/1D7 3330 EXC=FREAD(2,KX)/1D6 3340 AO=FREAD(2,KX)*1D1 3350 COMEGA=FREAD(2,KX)*RADG/1D4 3360 DCOMG=FREAD(2,KX)*RADG/1D7 3370 CI=FREAD(2,KX)/1D6 3380 GLAM=FREAD(2,KX)*RADG/1D4 3390 SI=FREAD(2,KX)/1D6 3400 C NASR C READ-IN AND PRINT-OUT THE UNSCALED VARIABLE ORBITAL PARAMETERS NASR C ( UP TO 17 SETS ,EACH OF WHICH INCLUDES 3 PARAMETERS CAN BE ACCEPTED)NASR C NASR PRINT5004 NASR3440 5004 FORMAT(/15X,'UNSCALED VARIABLE PARAMETERS'/) NASR3450 NUM=0 3460 DO 502 I=1,17 3470 DEK(I)=FREAD(1,KX) 3480 IF(DEK(I).EQ.EOS) GOTO 2502 3490 NUM=NUM+1 3500 DAK(I)=FREAD(1,KX) 3510 DETA(I)=FREAD(1,KX) 3520 PRINT 1502,DEK(I),DAK(I),DETA(I) 3530 1502 FORMAT(' ',10X,3F10.0) 3540 C NASR C SCALING THE SATELLITE VARIABLE ORBITAL PARAMETERS. NASR C NASR DEK(I)=DEK(I)*RADG/1D4 3550 DAK(I)=DAK(I)*1D1 3560 DETA(I)=DETA(I)*1D1 3570 502 CONTINUE 3580 2502 CONTINUE 3590 C NASR C***********************************************************************NASR C NASR C READ-IN FOR EACH INPUTTED STATION WHICH OBSERVED THE SATELLIE PASS NASR C THE FOLLOWING :- ITS NUMBER ,AND THE OBSERVED DOPPLERS TO THIS PASS NASR C ,THEN STORE THE DOPPLER DATA AT THE 'ICHOOZ'-STATION ONLY (NOTE THAT NASR C EITHER 'MAGNAVOX' OR 'ITT' FORMAT IS ACCEPTABLE) NASR C NASR NSKIP=0 3630 LOOP=0 3640 NUMD=NUM-3 3650 6310 N=FREAD(1,KX) 3660 IF(N.EQ.EOP) GOTO 6360 3670 LOOP=LOOP+1 3680 IF(LOOP.GT.7)GO TO 6355 3690 I=0 3700 DO 6320 K=1,NSTA 3710 IF(N.EQ.IDENT(K)) I=K 3720 6320 CONTINUE 3730 IF(I.EQ.1) GOTO 6340 3740 DO 6330 J=1,NUMD 3750 DUM1=FREAD(1,KX) 3760 DUM2=FREAD(1,KX) 3770 IF(I.EQ.0) GOTO 6330 3780 KDPSI(I,J)=0 3790 IF(DUM1.NE.0.AND.DUM2.NE.0) KDPSI(I,J)=J 3800 6330 CONTINUE 3810 GOTO 6310 3820 6340 DO 6350 J=1,NUMD 3830 DOP(1,J)=FREAD(1,KX) 3840 COUNT(J)=FREAD(1,KX) 3850 KDPSI(1,J)=0 3860 IF(DOP(1,J).NE.0.AND.COUNT(J).NE.0) KDPSI(1,J)=J 3870 6350 CONTINUE 3880 NSKIP=NSKIP+1 3890 GOTO 6310 3900 C NASR C DOPPLER DATA FROM MORE THAN 7 STATIONS - BAD NEWS (INPUT ERROR) NASR C ,NEGLECT THIS PASS COMPLETELY AND GO TO READ AND PROCESS ANOTHER PASSNASR C NASR 6355 PRINT 6356,LOOP 3940 6356 FORMAT(////8X,'THIS PASS IS REJECTED ,BECAUSE THERE ARE DOPPLER DANASR 1TA FROM',I4,3X,'TRACKING STATIONS TO THIS PASS') NASR 6357 DUM=FREAD(1,KX) 3960 IF(DUM.EQ.EOP) GO TO 10801 NASR GO TO 6357 3980 10801 IWRONG=IWRONG+1 NASR GO TO 1047 NASR 6360 IF(NSKIP.EQ.1)GO TO 6365 4000 C 4010 C ZERO OR MORE THAN ONE SET OF DATA FOR STATION ICHOOZ- BAD NEWS 4020 C ,NEGLECT THIS PASS COMPLETELY AND GO TO READ AND PROCESS ANOTHER PASSNASR C 4030 PRINT 6362,NSKIP 4040 6362 FORMAT(////8X,'THIS PASS IS REJECTED ,BECAUSE THERE ARE',I4,3X,'SENASR 1TS OF DOPPLER DATA AT THE SINGLE STATION-ICHOOZ') NASR IWRONG=IWRONG+1 NASR GO TO 1047 4060 C 4070 C REDUCE THE DOPPLER ARRAY FOR THE 'ICHOOZ'-STATION TO 8 ELEMENTS. NASR C 4090 6365 CONTINUE 4100 NRED=0 4110 K=NUMD-8 4120 IF(K.LE.0) GOTO 6380 4130 DO 6370 J=1,K 4140 IF(KDPSI(1,J).NE.0) GOTO 6380 4150 6370 NRED=NRED+1 4160 6380 DO 6390 J=1,8 4170 DOP(1,J)=DOP(1,J+NRED) 4180 COUNT(J)=COUNT(J+NRED) 4190 DO 6390 K=1,NSTA 4200 6390 KDPSI(K,J)=KDPSI(K,J+NRED) 4210 C NASR C REDUCE THE VARIABLE ORBITAL PARAMETERS ARRAY TO 11 ELEMENTS. NASR C NASR DO 6395 J=1,11 4220 DEK(J)=DEK(J+NRED) 4230 DAK(J)=DAK(J+NRED) 4240 6395 DETA(J)=DETA(J+NRED) 4250 C NASR C COMPUTE AND PRINT THE TIME IN UT MINUTES THAT EXCEEDS THE INPUTTED NASR C LOCKON-TIME ,IN CASE THAT THE VARIABLE PARAMETERS ARRAY EXCEEDS 11 NASR C ELEMENTS ; TO COMPUTE THE ACTUAL LOCKON-TIME. NASR C NASR MM=2*NRED 4260 PRINT 6397,MM 4270 6397 FORMAT('0',10X,'LOCKON NOW',I5,' MINUTES LATER THAN SHOWN ABOVE') 4280 C 4290 C************************************************************************* 4300 C 4310 C COMPUTE DTP=TIME FROM PERIGEE TO LOCKON 4330 C 4340 KTC=IH*60+M+2*NRED+0.1 4350 KTCX=KTC-4*(KTC/4) 4360 DTP=KTC-TP 4370 IF(DTP.LE.0) DTP=DTP+1440. 4380 IF(DTP+2.*PI/AN.GE.1440.) DTP=DTP-1440. 4390 C 4320 C INTERPOLATE THE OUT-OF-PLANE COMPONENTS (LAGRANGE INTERPOLATION) 4400 C KTCX IS ZERO IF KTC IS DIVISIBLE BY FOUR 4410 C 4420 M1=3 4430 M2=9 4440 IF(KTCX.EQ.0)GO TO 503 4450 M1=2 4460 M2=10 4470 503 DO 504 I=M1,M2,2 4480 DETA(I)=(DETA(I-1)+DETA(I+1))/2 4490 III=I-1 4500 IF(I.GE.9)III=I-3 4510 504 DETA(I)=DETA(I)-(DETA(III)+DETA(III+4)-2*DETA(III+2))/8 4520 C 4530 C REDUCE THE VARIABLE PARAMETERS ARRAYS FROM 11 TO 9 ELEMENTS. 4540 C 4550 DO 510 I=1,9 4560 DEK(I)=DEK(I+1) 4570 DAK(I)=DAK(I+1) 4580 510 DETA(I)=DETA(I+1) 4590 C NASR4600 C PRINT-OUT THE SCALED (FIXED & VARIABLE) SATELLITE ORBIT PARAMETERS NASR4610 C NASR4620 PRINT5005,TP,AN,OMEGA,DOMEGA,EXC,AO,COMEGA,DCOMG,CI,GLAM,SI NASR4630 5005 FORMAT(/5X,'SCALED FIXED PARAMETERS'/,11(/5X,F20.10)) NASR4640 PRINT5006 NASR4650 5006 FORMAT(/17X,'SCALED VARIABLE PARAMETERS'/) NASR4660 DO 5007 I=1,9 NASR4670 5007 PRINT5008,DEK(I),DAK(I),DETA(I) NASR4680 5008 FORMAT(10X,F10.8,2(5X,F10.4)) NASR4690 C NASR C***********************************************************************NASR C NASR C COMPUTATION OF THE SATELLITE GEOCENTRIC COORDINATES (XYZ) REFERRED NASR4720 C TO THE AVERAGE TERRESTRIAL SYSTEM - IN METERS. NASR4730 C 4740 DO 7105 K=1,9 4750 DTK=DTP+2.*(K-1) 4760 FMK=AN*DTK 4770 EK=FMK+EXC*DSIN(FMK)+DEK(K) 4780 AK=AO+DAK(K) 4790 C NASR C COMPUTE THE SATELLITE COORDINATES IN ITS ORBIT-PLANE SYSTEM. NASR C NASR P(1)=AK*(DCOS(EK)-EXC) 4800 P(2)=AK*DSIN(EK) 4810 P(3)=DETA(K) 4820 Q(1)=-(OMEGA-DABS(DOMEGA)*DTK) 4830 Q(3)=-(COMEGA-GLAM) -(DCOMG -WE)*DTK 4840 Q(2)=-PI/2.+DATAN(CI/SI) 4850 R(1)=3 4860 R(2)=1 4870 R(3)=3 4880 NO=3 4890 C NASR C CALLING THE SUBROUTINE 'DROTAT' TO COMPUTE THE PRODUCT ROTATION MATR-NASR C IX REQUIRED FOR SATELLITE COORDINATE TRANSFORMATION TO THE AVERAGE NASR C TERRESTRIAL SYSTEM. NASR C NASR CALL DROTAT(Q,R,NO,ROT) 4900 C NASR XS(K)=0. NASR4910 YS(K)=0. NASR4920 ZS(K)=0. NASR4930 DO 7106 I=1,3 4940 XS(K)=XS(K)+ROT(1,I)*P(I) 4950 YS(K)=YS(K)+ROT(2,I)*P(I) 4960 7106 ZS(K)=ZS(K)+ROT(3,I)*P(I) 4970 7105 CONTINUE 4980 PRINT203,(NDESC(I),I=1,18),NPASS NASR4990 C 203 FORMAT('1',/,20X,18A4,10X,'PASS NO. ',I5,/) NASR5000 C NASR C PRINT-OUT THE SATELLITE AVERAGE TERRESTRIAL CARTESIAN COORDINATES NASR C IN METERS. NASR C NASR PRINT 586 5010 586 FORMAT(/5X,'SATELLITE POSITIONS (XYZ-CARTESIAN GEOCENTRIC COORD.)NASR5020 1-IN THE AVER. TEERES. SYSTEM'//13X,'X',19X,'Y',19X,'Z',10X,'(UNITSNASR 2 ARE IN METERS)'/) NASR WRITE(IOUT,1600)(XS(I),YS(I),ZS(I),I=1,9) 5040 1600 FORMAT(' ',3F20.4) C NASR C***********************************************************************NASR5130 C NASR C INCREASE THE NUMBER OF THE DEGREES OF FREEDOM BY "3" ,FOR EACH PHASENASR C EXCEPT THE FIRST ONE ; WHICH IS DUE TO THE WEIGHTING OF THE (X,Y,Z) NASR C UNKNOWN PARAMETERS (COORDINATES OF THE TRACKING STATION "ICHOOZ" NASR C NASR I=1 5080 MDF=0 5090 526 IF(I.GT.NPTS) GO TO 528 5100 IF(NTEST(I).EQ.10) MDF=MDF+3 5110 GO TO 211 5120 C NASR C***********************************************************************NASR C NASR C OPTION CHARACTERISTICS FOR UNKNOWN STATION 5150 C (IF SOLUTION TURNS TO BE SINGULAR DUE TO NUIS. PAR. Y 5160 C PROGRAM SWITCHES TO ORBITAL SOLUTION ) 5170 C NOTE :- THESE POSSIBILITIES MIGHT BE USED FOR FUTURE VERSIONS ONLY. NASR C 5180 C NTEST(I)=10- UNKNOWN STA-NO SOLUTION FOR SPEED AND DIR. 5190 C NTEST(I)=12 SOLUTION FOR BOTH SPEED AND DIR. 5200 C NTEST(I)=13 SOLUTION FOR SPEED ONLY 5210 C NTEST(1)=14 SOLUITON FOR DIR. ONLY 5220 C NTEST(I)=15 CONTINUATION OF A PHASE ADJUSTMENT (I.E. SOLUTION OF 5230 C A PREVIOUS PHASE :- (VPV,MDFM,ISEQ,PX) ARE FEED IN) NASR C 5250 528 NTEST(1)=10 5260 IF(NTEST(I).NE.15)GO TO 107 5350 C NASR C***********************************************************************NASR C NASR C NOTE : THE FOLLOWING POSSIBILITY MIGHT BE USED IN FUTURE VERSIONS NASR C ONLY (IF REALLY NEEDED) ,BUT IT IS NEGLECTED IN THIS VERSION. NASR C READ-IN THE RESULTS (VPV,MDFM,ISEQ & PX) ,AS OBTAINED FROM A PHASED- NASR C ADJUSTMENT IN A PREVIOUS RUN OF THIS PROGRAM ; IF DESIRED TO CONTINUENASR C ON THIS PHASE-TYPE ADJUSTMENT. NASR C MDFM-TOTAL NO. OF DEGREES OF FREEDOM 5300 C VPV - " SQUARE SUM OF RESIDUALS 5310 C ISEQ-POSITIONAL RANK NO.(-SPECIFIES POSITION WITHIN PX (OR XX VECTOR)5320 C PX - THE WEIGHT MATRIX OF THE ESTIMATED RECEIVER'S COORDINATES NASR C (THE UPPER TRIANGLE MATRIX ONLY (BY ROWS) SINCE PX MATRIX NASR C IS A SYMMETRIC MATRIX. NASR C 5340 NTEST(I)=10 5360 MDF=MDF+3 5370 VPV=FREAD(1,KK) NASR MDFM=FREAD(1,KK) NASR ISEQ=FREAD(1,KK) NASR WRITE(IOUT,109)VPV,MDFM,ISEQ 5400 109 FORMAT('0',' PHASE VPV=',F20.4,//,' MDFM=',I4,//,' IS 5410 1EQ=',I4,/) 5420 VPV1=VPV 5430 JK=ISEQ+2 5440 JL=ISEQ+1 5450 C NASR C READ-IN THE UPPER TRIANGLE MATRIX OF PX-MATRIX NASR C NASR DO 112 J=ISEQ,JK NASR 112 PX(ISEQ,J)=FREAD(1,KK) NASR PX(JL,JL)=FREAD(1,KK) NASR PX(JL,JK)=FREAD(1,KK) NASR PX(JK,JK)=FREAD(1,KK) NASR C NASR C FILL-IN THE LOWER TRIANGLE MATRIX OF PX-MATRIX. NASR C NASR DO 110 L=ISEQ,JK 5490 J=L 5500 DO 110 K=J,JK 5510 110 PX(K,L)=PX(L,K) 5520 WRITE(IOUT,113)ISEQ,((PX(L,J),J=ISEQ,JK),L=ISEQ,JK) 5530 113 FORMAT('0',' PHASE WEIGT. COEF. MATRIX ISEQ=',I4,//,3E20.5,/,3E20 5540 1.5,/,3E20.5,//) 5550 C NASR C***********************************************************************NASR C NASR 107 CONTINUE 5560 211 CONTINUE 5600 DVPV=VPV 5610 DVPV1=VPV1 5620 C NASR C***********************************************************************NASR C NASR5640 C COMPUTE RANGES FROM THE OBSERVER TO THE SATELLITE POSITIONS NASR5650 C AS AN APPROXIMATE VALUE NEEDED IN THE ADJUTMENT TO EVALUATE THE NASR5660 C DESIGN MATRICES NASR5670 C NASR5680 DO 6050 J=1,9 NASR5690 PT(J,1)=XS(J)-X(I) NASR5700 PT(J,2)=YS(J)-Y(I) NASR5710 PT(J,3)=ZS(J)-Z(I) NASR5720 SK(I,J)=DSQRT(PT(J,1)**2+PT(J,2)**2+PT(J,3)**2) NASR5730 C NASR5740 C ROTATION OF THE RANGE VECTOR BETWEEN OBSERVER AND SATELLITE INTO THE NASR5750 C TOPOCENTRIC COORD. SYSTEM FOR EVALUATION OF ELEVATION ANGLES NASR5760 C NASR5770 QT(1)=ALOG(I)-PI NASR5780 QT(2)=ALAT(I)-PI/2. NASR5790 LRT(1)=3 NASR5800 LRT(2)=2 NASR5810 C NASR C CALLING THE SUBROUTINE 'DROTAT' TO COMPUTE THE PRODUCT ROTATION MATR-NASR C IX NECESSARY FOC COORDINATE TRANSFORMATION FROM AVERAGE TERRESTRIAL NASR C TO THE TOPOCENTRIC SYSTEM. NASR C NASR CALL DROTAT(QT,LRT,2,ROT) NASR5820 C NASR PTOP(J,1)=0. NASR5830 PTOP(J,2)=0. NASR5840 PTOP(J,3)=0. NASR5850 DO 5009 LV=1,3 NASR5860 DO 5009 KV=1,3 NASR5870 5009 PTOP(J,LV)=PTOP(J,LV)+PT(J,KV)*ROT(LV,KV) NASR5880 C NASR5890 C REFLECT THE SECONDARY AXIS - NAMELY PTOP(J,2) AXIS NASR5900 C NASR5910 PTOP(J,2)=-PTOP(J,2) NASR5920 C NASR5930 C COMPUTE ELEVATION ANGLES FROM OBSERVER TO SATELLITE POSITIONS NASR5940 C NASR5950 PTANG=DSQRT(PTOP(J,1)**2+PTOP(J,2)**2) NASR5960 VA(J)=DATAN2(PTOP(J,3),PTANG) NASR5970 ELEV(J)=VA(J)/RADG NASR5980 6050 CONTINUE NASR5990 C WELLS C***********************************************************************WELLS C WELLS C ...FIND TIME AND ELEVATION AT CLOSEST APPRAOCH FOR STN 'ICHOOZ' WELLS C ...FIRST FIND TCA WELLS C WELLS DO 7200 J=1,8 WELLS IF(ELEV(J+1).GE.ELEV(J)) GOTO 7200 WELLS IF(J.EQ.1) GOTO 7299 WELLS GOTO 7210 WELLS 7200 CONTINUE WELLS C ...TCA FOUND. COMPUTE IT. WELLS 7210 DTCA=(ELEV(J+1)-ELEV(J-1))/(ELEV(J)-AMIN1(ELEV(J+1),ELEV(J-1))) WELLS TCA=KTC+2*(J-1)+DTCA WELLS ITCAHR=TCA/60. WELLS TCAMIN=TCA-ITCAHR*60 WELLS JPM=J+1 WELLS IF(DTCA.LT.0.) JPM=J-1 WELLS C ...NOW FIND ELEVATION AT TCA. BEGIN BY INTERPOLATING VARIABLES WELLS EK=DEK(J)+(DTCA/2.)*(DEK(JPM)-DEK(J)) WELLS AK=DAK(J)+(DTCA/2.)*(DAK(JPM)-DAK(J)) WELLS ETA=DETA(J)+(DTCA/2.)*(DETA(JPM)-DETA(J)) WELLS C ...COMPUTE RANGE VECTOR IN AVERAGE TERRESTRIAL SYSTEM WELLS DTK=TCA-TP WELLS IF(DTK.LE.0.) DTK=DTK+1440. WELLS IF(DTK+2*PI/AN.GE.1440.) DTK=DTK-1440. WELLS FMK=AN*DTK WELLS EK=FMK+EXC*DSIN(FMK)+EK WELLS AK=AO+AK WELLS P(1)=AK*(DCOS(EK)-EXC) WELLS P(2)=AK*DSIN(EK) WELLS P(3)=ETA WELLS Q(1)=-(OMEGA-DABS(DOMEGA)*DTK) WELLS Q(2)=-PI/2.+DATAN(CI/SI) WELLS Q(3)=-(COMEGA-GLAM)-(DCOMG-WE)*DTK WELLS R(1)=3 WELLS R(2)=1 WELLS R(3)=3 WELLS CALL DROTAT(Q,R,3,ROT) WELLS XTCA=-X(I) WELLS YTCA=-Y(I) WELLS ZTCA=-Z(I) WELLS DO 7220 K=1,3 WELLS XTCA=XTCA+ROT(1,K)*P(K) WELLS YTCA=YTCA+ROT(2,K)*P(K) WELLS 7220 ZTCA=ZTCA+ROT(3,K)*P(K) WELLS C ...ROTATE RANGE VECTOR INTO TOPOCENTRIC SYSTEM WELLS P(1)=XTCA WELLS P(2)=YTCA WELLS P(3)=ZTCA WELLS Q(1)=ALOG(I)-PI WELLS Q(2)=ALAT(I)-PI/2. WELLS R(1)=3 WELLS R(2)=2 WELLS CALL DROTAT(Q,R,2,ROT) WELLS XTCA=0. WELLS YTCA=0. WELLS ZTCA=0. WELLS DO 7230 K=1,3 WELLS XTCA=XTCA+ROT(1,K)*P(K) WELLS YTCA=YTCA+ROT(2,K)*P(K) WELLS 7230 ZTCA=ZTCA+ROT(3,K)*P(K) WELLS C ...COMPUTE ELEVATION ANGLE 'ECA' WELLS PTANG=DSQRT(XTCA*XTCA+YTCA*YTCA) WELLS ECA=DATAN2(ZTCA,PTANG)/RADG WELLS PRINT 7240,ITCAHR,TCAMIN,ECA WELLS 7240 FORMAT(75X,'CLOSEST APPROACH AT',I5,F8.3,/,75X,'MAXIMUM ELEVATION' 1 , F8.3) 7299 CONTINUE WELLS C WELLS C***********************************************************************WELLS C WELLS C NASR C FIND THE MAXIMUM ELEVATION (IN DEGREES) OF THIS PASS WITH RESPECT TO NASR C THE TRACKING STATION 'ICHOOZ' NASR C NASR IB=1 6030 DO 533 J=2,9 6040 IF(ELEV(J).GE.ELEV(J-1).OR.IB.GT.1) GOTO 533 6050 ELMAX(I)=(ELEV(J)-ELEV(J-2))**2/(8.0*(2.0*ELEV(J-1)-ELEV(J)-ELEV(J 6060 1-2)))+ELEV(J-1) 6070 IB=IB+1 6080 533 CONTINUE 6090 C NASR C***********************************************************************NASR C NASR C CORRECT THE OBSERVED DOPPLERS AT ICHOOZ-STATION ; FOR BOTH IONOSPHER-NASR C IC AND TROPOSPHERIC REFRACTIONS. NASR C NASR PRINT 6150 NASR6150 6150 FORMAT('-',27X,'DOPPLER DATA'//6X,'ELEV',8X,'HI DOP',8X,'LO DOP',8NASR6160 1X,'IONOS',8X,'TROPO') NASR6170 I=1 6190 IA=1 6200 NG(I)=0 6210 DO 6051 J=1,8 6220 DOP1=DOP(I,J) NASR DOP2=0.D0 6240 KDPSI(I,J)=0 6250 IF(DOP(I,J).EQ.0.OR.COUNT(J).EQ.0) GOTO 524 6260 KDPSI(I,J)=J 6270 C 6290 C IONOSPHERIC REFRACTION CORRECTION 6300 C 6310 IF(COUNT(J).LT.10000.)GO TO 514 6320 C 6330 C MAGNAVOX FORMAT 6340 C 6350 DOP(I,J)=DOP(I,J)+(DOP(I,J)-COUNT(J))*9./55. 6360 GO TO 525 6370 C 6380 C ITT FORMAT NASR C 6400 514 DOP(I,J)=DOP(I,J)-(COUNT(J)-2000.)*24./55. 6410 525 IA=IA+1 6420 DOP2=DOP(I,J) 6430 C 6440 C TROPOSPHERIC REFRACTION CORRECTION ; USING 'HOPFIELD'-MODEL (1969). NASR C 6460 IF(ITROP.NE.1) GOTO 524 6470 ELEVD1=DSQRT(VA(J)**2+THD) 6480 ELEVD2=DSQRT(VA(J+1)**2+THD) 6490 ELEVW1=DSQRT(VA(J)**2+THW) 6500 ELEVW2=DSQRT(VA(J+1)**2+THW) 6510 TROP2=XKD/DSIN(ELEVD2)+XKW/DSIN(ELEVW2) 6520 TROP1=XKD/DSIN(ELEVD1)+XKW/DSIN(ELEVW1) 6530 DOP(I,J)=DOP(I,J)-WAVENO*(TROP2-TROP1) 6540 524 CONTINUE 6550 IF(DOP(I,J).EQ.0.OR.COUNT(J).EQ.0) DOP(I,J)=0.0D 0 NASR C NASR C PRINT-OUT THE DOPPLER DATA BEFORE AND AFTER APPLYING THE REFRACTION NASR C CORRECTIONS. NASR C NASR PRINT 6053,ELEV(J),KDPSI(I,J),DOP1,COUNT(J),DOP2,DOP(I,J) 6560 6053 FORMAT(' ',F10.2/10X,I2,4F14.2) NASR6570 6051 CONTINUE 6580 PRINT 6053,ELEV(9) 6590 C NASR C***********************************************************************NASR6600 C NASR C FIND THE NUMBER (NG) OF THE USABLE NON-ZERO DOPPLERS AT THE ICHOOZ- NASR C STATION ; AFTER THE HORIZON AND TRIMMING REJECTIONS. NASR C 6630 NG(I)=IA-1 6640 C 6650 C REJECT HORIZON DOPPLERS AND COUNT THE REMAINING USEFUL DOPPLERS NASR6670 C NASR6680 NREJ=0 6690 DO 6080 J=1,8 6710 IF(KDPSI(I,J).EQ.0) GOTO 6080 6720 IF(ELEV(J).GT.HOR) GOTO 6073 6730 GOTO 6074 6740 6073 IF(ELEV(J+1).GT.HOR) GOTO 6080 6750 6074 KDPSI(I,J)=0 6760 NREJ=NREJ+1 6770 6080 CONTINUE 6780 NG(I)=NG(I)-NREJ 6790 PRINT 6081,NREJ,HOR 6800 6081 FORMAT('0',10X,I10,10X,'DOPPLERS REJECTED ON ',F5.1,5X, 6810 1'DEGREES ELEVATION (HORIZON CUT-OFF CRITERION)') NASR C NASR C TRIM TO ONLY THE SIMULTANEOUS DOPPLERS THAT OCCURED FROM ALL INPUTTEDNASR C TRACKING STATIONS THAT OBSERVED THIS PASS. (IF REQUIRED) NASR C 6850 NTRIM=0 6860 IF(ITRIM.NE.1) GOTO 6085 6870 DO 6084 K=1,NSTA 6880 DO 6084 J=1,8 6890 IF(KDPSI(K,J).NE.0) GOTO 6084 6900 KDPSI(I,J)=0 6910 6084 CONTINUE 6920 DO 6087 J=1,8 6930 IF(KDPSI(I,J).NE.0) NTRIM=NTRIM+1 6940 6087 CONTINUE 6950 NTRIM=NG(I)-NTRIM 6960 6085 NG(I)=NG(I)-NTRIM 6970 PRINT 6086,NTRIM,ITRIM 6980 6086 FORMAT('0',10X,I10,10X,'DOPPLERS REJECTED TO TRIM TO ONLY', 6990 1 ' SIMULTANEOUS DOPPLERS WITH ITRIM= ',I5) 7000 C NASR7010 C REJECT PASSES HAVING LESS THAN TWO DOPPLERS NASR7020 C NASR7030 PRINT5014,IDAY,IPASS,ISAT NASR7040 5014 FORMAT('0',5X,'PASS ON DAY ',I5,5X,'PASS NUMBER',I5,5X,'SATELLITE NASR7050 1NUMBER',I5) NASR7060 IF(NG(I).GE.2) GO TO 5015 NASR7070 PRINT 5016 NASR7080 5016 FORMAT('+',70X,'***REJECTED-HAS LESS THAN 2 DOPPLERS***') NASR IDOPL2=IDOPL2+1 NASR GO TO 104 NASR7110 5015 PRINT 5017 NASR7120 5017 FORMAT('+',70X,'***ACCEPTED-HAS MORE THAN 2 DOPPLERS***') NASR PRINT 5018 NASR7140 5018 FORMAT('0',5X,'ST. NO.',5X,'MAX. ELEV.',8X,'DOPPLERS TO BE USED',1NASR7150 10X,'AVAILABLE DOPPLER EQUATIONS'/) NASR7160 NGRPTS=1 NASR DO 5019 JK=1,NGRPTS NASR7170 5019 PRINT5020,ISTN,ELMAX(JK),(KDPSI(JK,JL),JL=1,8),NG(JK) NASR7180 5020 FORMAT(6X,I4,10X,F6.3,5X,8I4,10X,I2) NASR7190 C NASR C***********************************************************************NASR C NASR C CHECK ; TO BE SURE THAT THE FIRST ACCEPTED PASS ,UNTIL THIS STAGE , NASR C HAS AT LEAST 4 USABLE DOPPLERS ,TO FIND A SOLUTION ; OTHERWISE , NASR C THIS PASS SHOULD BE REJECTED ,AFTER PRINTING-OUT A MESSAGE. NASR C NASR IF(KPASS.EQ.0.AND.NG(JK).LT.4) GO TO 45454 NASR GO TO 42526 NASR 45454 PRINT 203,(NDESC(I),I=1,18),NPASS NASR C 203 FORMAT('1',/,20X,18A4,10X,'PASS NO. ',I5,/) NASR PRINT 64719 NASR 64719 FORMAT(' '///,5X,'THIS PASS IS REJECTED ,BECAUSE THE NUMBER OF DOPNASR 1PLER EQUATIONS IS LESS THAN THE NUMBER OF UNKNOWNS'/) NASR INSING=INSING+1 NASR GO TO 1047 NASR 42526 CONTINUE NASR C NASR C***********************************************************************NASR7380 C NASR7390 C REINIALIZING ZERO VALUES FOR SOME OF THE INVOLVED ARRAYS IN THE NASR C PHASED ADJUSTMENT PROCESS. NASR C NASR DO 102 IZERO=1,MROW NASR W(IZERO)=0. NASR V(IZERO)=0. NASR DO 103 JZERO=1,MCCOL NASR 103 C(IZERO,JZERO)=0. NASR DO 102 KZERO=1,MACOL NASR 102 A(IZERO,KZERO)=0. NASR DO 106 LZERO=1,MACOL NASR U(LZERO)=0. NASR XX(LZERO)=0. NASR DO 106 MZERO=1,MACOL NASR 106 ANORM(LZERO,MZERO)=0. NASR C NASR C***********************************************************************NASR C NASR C FORMING THE DESIGN MATRICES ; COEFFICIENT MATRIX-A AND MISCLOSURE NASR C VECTOR-W ,AS A PRE-REQUISITE FOR THE PHASED LEAST SQUARES ADJUSTMENT NASR C NASR7410 I=I+1 7240 527 NGRPTS=I-1 ITER=0 7300 SPEEDM=0. 7310 GO TO 572 7320 573 DO 571 J=1,NGRPTS 7340 DO 571 I=1,9 7350 571 SK(J,I)=DSQRT((XS(I)-X(J))**2+ (YS(I)-Y(J))**2+(ZS(I)-Z(J))**2) 7360 572 IJI=1 7420 KK=0 7430 IC=1 7440 JDI=1 7450 KL=0 7460 IF(NGRPTS.GT.3) PRINT 1601 7280 1601 FORMAT('0',//,' M A T R I X C',//,' ROW COL ',//) 7290 DO 512 II=1,NGRPTS 7480 IF(NTEST(II).GE.10) KK=KK+1 7490 XN(II)=X(II) 7500 YN(II)=Y(II) 7510 ZN(II)=Z(II) 7520 IDI=0 7530 NPTS=NGRPTS 7580 DO 512 J=1,8 7540 IF(KDPSI(II,J).EQ.0)GO TO 519 7550 IB=IB+1 7560 IF(SPEED(II).NE.0.) GO TO 535 7570 C 7590 C FORMING THE W-VECTOR : ABSOLUTE TERMS OF THE OBSERVATION EQUATIONS NASR7600 C 7610 521 W(IC)=-(SK(II,J)-SK(II,J+1))/LO+2.*FDI(II) -DOP(II,J) 7620 IF(SPEED(II).NE.0.) GO TO 522 7630 C 7640 C FORMING MATRIX-C ; COEFFICIENT MATRIX OF THE UNKNOWN SATELLITE COORD-NASR C INATES (IN CASE OF ORBIT UNKNOWN MODE) ,OTHERWISE IT WILL BE USED IN NASR C COMPUTING THE ELEMENTS OF THE MATRIX-A. NASR C 7660 535 JD=J*3-2 7670 C IF(JD.GT.JDI) IDI=JD -JDI 7680 C JD=JD-IDI 7690 IF(JD.GT.19) JD=19 7700 DIS=SK(II,J)*LO 7710 C(IC,JD)= (X(II)-XS(J))/DIS 7720 C(IC,JD+1)= (Y(II)-YS(J))/DIS 7730 C(IC,JD+2)= (Z(II)-ZS(J))/DIS 7740 IF(SPEED(II).NE.0.) GO TO 519 7750 522 DIST=SK(II,J+1)*LO 7760 C(IC,JD+3)=-(X(II)-XS(J+1))/DIST 7770 C(IC,JD+4)=-(Y(II)-YS(J+1))/DIST 7780 C(IC,JD+5)=-(Z(II)-ZS(J+1))/DIST 7790 JDC=JD+5 7800 IF(NGRPTS.GT.3)WRITE(IOUT,1610) IC,JD,(C(IC,JS),JS=JD,JDC) 7810 1610 FORMAT(' ',I4,I4,2(4X,3F8.4)) 7820 IF(NTEST(II).LT.10)GO TO 523 7830 C 7840 C FORMING MATRIX-A ; COEFFICIENT MATRIX OF THE DESIRED PARAMETERS (I.E.NASR C THE COORDINATES OF THE GROUND TRACKING STATION ) NASR C 7860 101 KI=KL+3*KK-2 7870 A(IC,KI)= -C(IC,JD) +(X(II)-XS(J+1))/DIST 7880 A(IC,KI+1)= -C(IC,JD+1)+(Y(II)-YS(J+1))/DIST 7890 A(IC,KI+2)= -C(IC,JD+2)+(Z(II)-ZS(J+1))/DIST 7900 IF(SPEED(II).NE.0.) GO TO 537 7910 ISEQ=KI 7920 GO TO 523 7930 C NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C NASR C COMPUTE THE (X,Y,Z) INCREMENTS DUE TO THE SPEED OF THE RECEIVER (FOR NASR C THE MOVING TRACKING STATION - IF ANY) NASR C NOTE :- THIS PART IS NEGLECTED IN OUR CASE OF PHASED-SOLUTION FOR NASR C A STATIONARY TRACKING STATION 'ICHOOZ' NASR C NASR 519 IF(SPEED(II).EQ.0.) GO TO 512 7980 SPEEDM=SPEED(II)*1852./60. 7990 DIRCT=DIR(II)*RADG 8000 AX(1)=2.*SPEEDM*COS(DIRCT) 8010 AX(2)=-2.*SPEEDM*SIN(DIRCT) 8020 AX(3)=0. 8030 WF=DSQRT(1.-ESQ*DSIN(ALAT(II))**2) 8050 AXXK2= AX(2)*(J*2-1)/2 8060 AXXL1= AX(1)*(J*2-1)/2 8070 CR(2)=PI-(ALOG(II)-AXXK2*WF/(A1*DCOS(ALAT(II)))) 8080 CR(1)=PI/2.-(ALAT(II)+AXXL1*WF**3/(A1*(1.-ESQ))) 8090 LR(1)=2 8100 LR(2)=3 8110 C 7970 C ROTATION MATRIX 8120 C 8130 CALL DROTAT(CR,LR,2,ROT) 8140 C 7950 IF(J.NE.1) GO TO 529 8150 NPTS=NGRPTS-1 8160 529 DO 520 IA=1,3 8170 X(II)=X(II)+ROT(1,IA)*AX(IA) 8180 Y(II)=Y(II)+ROT(2,IA)*AX(IA) 8190 520 Z(II)=Z(II)+ROT(3,IA)*AX(IA) 8200 IA=J+1 8210 SK(II,IA)=DSQRT((X(II)-XS(IA))**2+(Y(II)-YS(IA))**2+(Z(II)-ZS(IA)) 8220 1**2) 8230 IF(KDPSI(II,J).EQ.0)GO TO 512 8240 GO TO 521 8250 537 IF(NTEST(II).EQ.10)GO TO 523 8260 KE=JD+2 8270 DO 551 JAN=JD,KE 8280 JN=JAN-JD+1 8290 RAT=RAT-2.*(C(IC,JD)*(J-1)+C(IC,KE+JN)*J)*ROT(JN,1) 8300 551 RAB=RAB-2.*(C(IC,JD)*(J-1)+C(IC,KE+JN)*J)*ROT(JN,2) 8310 KZ=KI+3 8320 IF(NTEST(II).EQ.14)GO TO 550 8330 A(IC,KZ)= (RAT*COS(DIRCT)-RAB*SIN(DIRCT))*60./1852. 8340 KL=KL+1 8350 IF(NTEST(II).EQ.13)GO TO 523 8360 KZ=KZ+1 8370 550 A(IC,KZ)= -(RAT*SIN(DIRCT)*SPEEDM+RAB*COS(DIRCT)*SPEEDM)/RA 8380 1DG 8390 KL=KL+1 8400 C NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C NASR 523 JDI=JD+6 8430 IC=IC+1 8440 IF(I.GT.IJI)IJI=I 8450 512 CONTINUE 8460 C NASR C***********************************************************************NASR C NASR C DIMENSIONS OF THE A & C MATRICES ; (BEFORE INCLUDING THE COLUMNS NASR C PERTAINING TO THE FREQUENCY OFFSET UNKNOWNS. NASR C NASR MCCOL=JDI-1 8500 MACOL=(KK)*3+KL 8510 MROW=IC-1 8520 C 8530 C INCLUDE THE COEFFICIENTS OF THE FREQUENCY OFFSET UNKNOWNS AS COLUMNS NASR C IN THE A-MATRIX ,CONSIDERING THEM AS DESIRED PARAMETERS TO BE SOLVED NASR C 8550 IFF=0 8560 IC=1 8570 DO 536 IA=1,NGRPTS 8580 IF(NG(IA).EQ.0) IFF=IFF+1 8590 DO 536 IB=1,8 8600 IF(KDPSI(IA,IB).EQ.0)GO TO 536 8610 IE=MACOL+IA-IFF 8620 A(IC,IE)=0.20D 00 NASR8630 IC=IC+1 8640 536 CONTINUE 8650 MACLM=MACOL+NGRPTS 8660 MACOL=MACLM-IFF 8670 C NASR8730 C PRINT-OUT THE DESIGN MATRICES 'A' , 'W' NASR8740 C NASR8750 PRINT 203,(NDESC(I),I=1,18),NPASS NASR8760 C 203 FORMAT('1',/,20X,18A4,10X,'PASS NO. ',I5,/) NASR8770 WRITE(IOUT,1602) 8780 1602 FORMAT(/,035X,'M A T R I X - A',38X,'M A T R I X - W'/) NASR879 C 8800 DO 8505 JK=1,MROW 8810 WRITE(IOUT,538) (A(JK,JL),JL=1,MACOL) 8820 8505 WRITE(6,8506) W(JK) 8830 538 FORMAT(15D20.10) 8840 8506 FORMAT('+',85X,D20.10) 8850 C C RESET VPV & VPV1 C VPV=DVPV 8710 VPV1=DVPV1 8720 C NASR C***********************************************************************NASR8870 C***********************************************************************NASR8880 C NASR C***** P H A S E A D J U S T M E N T C O M P U T A T I O N ******** NASR C NASR C CALLING THE SUBROUTINE 'A D J U S T' TO PERFORM THE PARAMETRIC LEAST NASR C SQUARES ADJUSTMENT OF THE DOPPLER SATELLITE OBSERVATIONS FROM THIS NASR C PASS ,INCORPORATING THE WEIGHTING TECHNIQUE OF THE UNKNOWN PARAMETERSNASR C OF THE 'ICHOOZ'-STATION ONLY. NASR C NASR CALL ADJUST 8920 C NASR C CHECK ; IF THERE IS NO AVAILABLE SOLUTION FOR THIS PASS ,DUE TO THE NASR C SINGULARITY OF THE NORMAL EQUATIONS MATRIX ? - PRINT A MESSAGE , NASR C REJECT THIS PASS , AND GO BACK TO READ AND PROCESS ANOTHER NEW PASS. NASR C NASR IF(D1.EQ.0.) GO TO 71117 NASR GO TO 92229 NASR 71117 PRINT 81118,NG(1) NASR 81118 FORMAT(//,5X,'THE USED DOPPLER EQUATIONS ARE =',I5,3X,'CHECK YOUR NASR 1DATA INPUT '/) NASR INSING=INSING+1 NASR GO TO 1047 NASR 92229 CONTINUE NASR C NASR C***********************************************************************NASR C NASR C UPDATING THE APPROXIMATE VALUES OF THE PARAMETERS (WHICH WERE THE NASR C ESTIMATED VALUES FROM THE PREVIOUS PHASE SOLUTION) ,BY ADDING THE NASR C SOLUTION VECTOR OBTAINED THROUGH THIS PHASE. NASR NROW=0 2810 NCOL=4 2820 MAX=0 8970 ITER=ITER+1 8980 IFF=0 8990 JROW=1 9000 DO 556 JK=1,NGRPTS 9010 JT=JK-IFF 9020 IF(NG(JK).NE.0) GO TO 701 9030 IFF=IFF+1 9040 FDI(JK)=FO 9050 GO TO 556 9060 C NASR C UPDATING THE FREQENCY OFFSET UNKNOWN (IN CYCLES/MIN.) NASR C NASR 701 FDI(JK)=FDI(JK)+XX(MACLM-NGRPTS+JT)/10.00000 9070 IF(NTEST(JK).LT.10) GO TO 556 9080 IF(NROW.EQ.1) GO TO 556 9090 IF(SPEED(JK).EQ.0.)GO TO 581 9100 C NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C NASR IF(DABS(XX(JROW)).GT.50.)MAX=50 9110 IF(DABS(XX(JROW+1)).GT.50.)MAX=50 9120 IF(DABS(XX(JROW+2)).GT.50.)MAX=50 9130 IF(DABS(XX(MACLM-NGRPTS+JT)).GT.50.)MAX=50 9140 C NASR C UPDATING THE APPROXIMATE COORD. OF THE TRACKING STATION ,IF IT WAS NASR C IN MOTION. NASR C NASR X(JK)=XN(JK)+XX(JROW) 9150 Y(JK)=YN(JK)+XX(JROW+1) 9160 Z(JK)=ZN(JK)+XX(JROW+2) 9170 C NASR C UPDATING THE SPEED ,IF IT WSS REQUIRED TO BE SOLVED FOR. NASR C NASR JROW=JROW+3 9180 IF(NTEST(JK).EQ.10) GO TO 556 9190 IF(NTEST(JK).EQ.14) GO TO 557 9200 SPEED(JK)=SPEED(JK)+XX(JROW) 9210 C NASR C UPDATING THE DIRECTION (OR COURSE) OF THE MOVING STATION ,IF IT NASR C WAS REQUIRED TO BE SOLVED FOR. NASR C NASR JROW=JROW+1 9220 IF(NTEST(JK).EQ.13) GO TO 556 9230 557 DIR(JK)=DIR(JK)+XX(JROW) 9240 DIR(JK)=DIR(JK)+XX(JROW) 9250 JROW=JROW+1 9260 GO TO 556 9270 C NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C NASR 581 JROW=JROW+3 9280 556 CONTINUE 9290 IF(NROW.EQ.1) MAX=0 9310 IF(ITER.EQ.10) GO TO 591 9320 IF(MAX.GE.50) GO TO 573 9330 JROW=1 9370 DO 584 JK=1,NGRPTS 9430 IF(NTEST(JK).LT.10) GO TO 584 9470 C NASR C COMPUTE THE UPDATED FREQUENCY OFFSET IN CYCLES/SECOND. ,FOR PRINTOUT NASR C NASR FOFF=FDI(JK)/60. 9440 IF(NROW.EQ.1) GO TO 584 9480 IF(SPEED(JK).NE.0.)GO TO 585 9490 C NASR C UPDATING THE TRACKING STATION COORD. ,IF THE STATION IS MOTIONLESS NASR C NASR X(JK)=XN(JK)+XX(JROW) 9530 Y(JK)=YN(JK)+XX(JROW+1) 9540 Z(JK)=ZN(JK)+XX(JROW+2) 9550 585 JROW=JROW+3 9560 IF(NTEST(JK).EQ.12) JROW=JROW+2 9570 IF(NTEST(JK).LE.12) GO TO 584 9580 JROW=JROW+1 9590 584 CONTINUE 9600 C NASR C***********************************************************************NASR C NASR C PRINT-OUT THE WEIGHT COEFFICIENT MATRIX (THE INVERSE OF THE NORMAL NASR C EQUATIONS) AND THE ABSOLUTE TERM OF THE NORMAL EQUATIONS (U-VECTOR) NASR C NASR WRITE(IOUT,154) NASR8330 154 FORMAT(//6X,'WEIGHT COEFFICIENT MATRIX OF THE ADJUSTED PARAMETERS NASR 1') NASR PRINT5022 NASR 5022 FORMAT('+',69X,'U - VECTOR'/) NASR MCC=MACOL 8350 IF(MACOL.GT.9) MCC=MACOL-NGRPTS 8360 DO 163 I=1,MACOL NASR8370 WRITE(IOUT,150) (ANORM(I,J),J=1,MCC) NASR8380 163 PRINT5023,U(I) NASR 5023 FORMAT('+',62X,D20.10) NASR 150 FORMAT(' ',3(3X,3E13.5)) NASR IF(MACOL.LE.9) GO TO 1612 8400 WRITE(IOUT,1611) 8410 1611 FORMAT('0',//,' MATRIX CONTINUED',//) 8420 MCD=MCC+1 8430 DO 1613 I=1,MACOL 8440 1613 WRITE(IOUT,150) (ANORM(I,J),J=MCD,MACOL) 8450 C150 FORMAT(' ',3(3X,3E13.5)) NASR C 9610 C PRINT SOLUTION VECTOR X 9620 C 9630 1612 WRITE(IOUT,156) 9640 156 FORMAT('0'//5X,'S O L U T I O N V E C T O R X (UNITS ARE IN METENASR 1RS EXCEPT LAST ELEMENT IN CYCLES/10 MIN.)'/) NASR WRITE(IOUT,157) (XX(I),I=1,MACOL) 9660 157 FORMAT(' ',F15.3) 9670 C 9680 C PRINT-OUT THE ORIGINAL APPROX. AND THE LATEST ADJUSTED VALUES OF NASR9690 C THE RECEIVER GEOCENTRIC COORDINATES NASR9700 C NASR9710 582 FORMAT(///10X,'THIS PASS IS REJECTED AT 95% PROBABILITY LEVEL') NASR ICHI95=ICHI95+1 NASR IF(SPEEDM.EQ.0.) GO TO 104 0890 598 IF(ISEQ.EQ.0) GO TO 594 0900 C NASR C***********************************************************************NASR C NASR0950 C COMPUTATION OF THE WEIGHT MATRIX OF THE LATEST ESTIMATE OF THE NASR0960 C RECEIVER COORDINATES TO BE USED IN THE FOLLOWING PHASE NASR0970 C THIS WEIGHT MATRIX IS CALLED THE 'PHASE COEFFICIENT MATRIX' NASR0980 C (ONLY ONE STATIONARY TRACKING STATION IS ALLOWED TO BE ADJUSTED BY NASR C THE PHASED-SOLUTION.) NASR C 1110 DO 592 ICK=1,3 1010 DO 592 JCK=1,3 1020 592 PX(ICK,JCK)=ANORM(ICK,JCK) 1030 C NASR C CALLING THE SUBROUTINE 'DINVPV' TO INVERT THE WEIGHT COEFFICIENT MAT-NASR C RIX OF THE LATEST ESTIMATE OF THE RECEIVER COORD. ,TO OBTAIN THEIR NASR C WEIGHT MATRIX TO BE USED IN THE NEXT PHASE. NASR C NASR CALL DINVPV(PX,3,DC,4) 1040 C 594 IF(MDFM.EQ.0) PRINT 79053 NASR 79053 FORMAT(////7X,'NOTE :- THAT THERE ARE NO ACCURACY ESTIMATES TO BE NASR 1COMPUTED FOR THIS PASS,'/,15X,'BECAUSE THE ACCUMULATED NUMBER OF NASR 2DEGREES OF FREEDOM EQUALS ZERO FOR THIS PHASE') NASR IF(MDFM.EQ.0) GO TO 508 NASR C NASR C PRINT-OUT THE ACCUMULATED NUMBER OF DEGREES OF FREEDOM (MDFM) NASR C NASR WRITE(IOUT,561) MDFM 1150 561 FORMAT(' '/5X,'ACCUMULATED NO. OF DEGREES OF FREEDOM =',I4//05X,'ENASR1160 1STIMATED VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED PARAMETERS'/) NASR C NASR C COMPUTING THE ESTIMATED VARIANCE FACTOR NASR1120 C 1130 SIGM=VPV/MDFM 1180 C NASR1190 C ESTIMATING AND PRINTING-OUT THE VARIANCE-COVARIANCE MATRIX OF THE NASR C ADJUSTED PARAMETERS. NASR C NASR1210 DO 164 JK=1,MACOL 1220 DO 560 IK=1,MACOL 1230 560 ANORM(JK,IK)=SIGM*ANORM(JK,IK) 1240 164 WRITE(IOUT,150) (ANORM(JK,IK),IK=1,MACOL) 1250 C150 FORMAT(' ',3(3X,3E13.5)) NASR IF(MACOL.LE.9) GO TO 1614 1420 WRITE(IOUT,1611) 1430 C1611 FORMAT('0',//,' MATRIX CONTINUED',//) 1440 DO 1615 I=1,MACOL 1450 1615 WRITE(IOUT,150) (ANORM(I,J),J=MCD,MACOL) 1460 C150 FORMAT(' ',3(3X,3E13.5)) NASR 1614 CONTINUE 1480 C NASR C PRINT-OUT THE ESTIMATED VARIANCE FACTOR. NASR C NASR WRITE(IOUT,562) SIGM 1490 562 FORMAT(//,5X,'ESTIMATED VARIANCE FACTOR =',F16.3) NASR C NASR C COMPUTE AND PRINT THE STANDARD DEVIATION OF UNIT WEIGHT (SQUARE ROOT NASR C OF THE ESTIMATED VARIANCE FACTOR) NASR C NASR SIGM=DSQRT(SIGM) NASR PRINT 66619,SIGM NASR 66619 FORMAT(5X,'STANDARD DEVIATION OF UNIT WEIGHT =',F8.3,3X,'(CYCLES)'NASR 1/) NASR C NASR1510 C COMPUTE AND PRINT THE STANDARD DEVIATIONS OF THE ESTIMATED COORDINAT-NASR C ES OF THE TRACKING STATION , IN METERS UNITS. NASR C NASR1530 SIGX=DSQRT(ANORM(1,1)) 1540 SIGY=DSQRT(ANORM(2,2)) 1550 SIGZ=DSQRT(ANORM(3,3)) 1560 PRINT 6100,SIGX,SIGY,SIGZ NASR 6100 FORMAT(/5X,'STANDARD DEVIATIONS OF (X,Y & Z) COORD. (UNITS ARE IN NASR 1METERS)'//8X,'FOR X',F15.3/8X,'FOR Y',F15.3/8X,'FOR Z',F15.3) NASR C NASR C***********************************************************************NASR C NASR C COMPUTE THE ACTUAL NUMBER OF THE ACCEPTED PASSES AFTER ALL REJECTIONSNASR C NASR KPASS=NPASS-IWRONG-IDOPL2-INSING-ICHI95 NASR C NASR C STORE THE LATEST ESTIMATED VALUES OF THE RECEIVER COORDINATES AFTER NASR C EVERY (INTEIG) PASSES ,TO BE PUNCHED-OUT WITH THEIR ESTIMATED VAR.- NASR C COV. MATRIX IF DESIRED (WHEN IPUNCH = 1) NASR C NASR IOK=0 9920 KWW=KPASS/INTEIG NASR RWW=KPASS NASR AWW=RWW/AINT 9960 TEST=AWW-KWW 9970 IF(TEST.GT.1.0D-08)GO TO 508 NASR IOK=1 9990 IF(IPUNCH.EQ.0)GO TO 9549 NASR JK=1 NASR STORE(KWW,1)=X(JK) 0000 STORE(KWW,2)=Y(JK) 0010 STORE(KWW,3)=Z(JK) 0020 C 1270 C PUNCH COVARIANCE MATRIX FOR THE ESTIMATED RECEIVER COORDINATES AFTER NASR1280 C EVERY INTEIG PASS (IF REQUIRED) 1290 C 1300 IF(IPUNCH.EQ.1) WRITE(7,7000) ((ANORM(JC,IK),IK=1,3),JC=1,3) 1340 7000 FORMAT(3E13.5) 1350 IF(IPUNCH.EQ.2) WRITE(7,7002) ((ANORM(JC,IK),IK=1,3),JC=1,3) 1360 7002 FORMAT(39X,3E13.5) 1370 C NASR C EIGEN-IZE THE (3,3) VARIANCE-COVARIANCE MATRIX OF THE RECEIVER COORD.NASR1630 C AFTER EVERY (INTEIG) PASS AND AFTER THE LAST PASS ,TO COMPUTE THE NASR C STANDARD ERROR ELLIPSOID FOR THE RECEIVER POSITION. NASR C 1660 9549 CONTINUE 1680 IF(IDAY.EQ.EOD.AND.IOK.EQ.1)GO TO 507 1700 IF(KPASS.EQ.0) GO TO 507 NASR DO 9500 IPK=1,3 1710 DO 9500 JPK=1,3 1720 9500 EX(IPK,JPK)=ANORM(IPK,JPK) 1730 C NASR C CALL THE SUBROUTINE 'MEIGEN' FOR PERFORMING THE EIGEN VALUE PROBLEM NASR C ON EX-MATRIX. NASR C NASR CALL MEIGEN(EX,EVECT,3,3,3) 1740 C NASR WRITE(6,9502) INTEIG 1750 DO 9501 IPK=1,3 1760 C NASR C COMPUTE AND PRINT THE SEMI-AXES OF THE ERROR ELLIPSOID ,IN METERS. NASR C NASR SEMI=DSQRT(EX(IPK,IPK)) 1770 9501 WRITE(6,9503) SEMI 1780 9502 FORMAT(//,5X,'SEMI-AXES OF THE STANDARD ERROR ELLIPSOID : IN METERNASR 1S'/,5X,'PRINTED-OUT AFTER EVERY',I5,2X,'ACCEPTED PASSES AND AFTER NASR 2THE LAST ACCEPTED PASS'/) NASR 9503 FORMAT(13X,F15.3) NASR C 1830 C END OF EIGEN VALUE PROBLEM 1840 C 1850 IF(IDAY.EQ.EOD)GO TO 507 1820 508 IF(SPEEDM.EQ.0)GO TO 104 1860 GO TO 547 1870 C NASR C***********************************************************************NASR1890 C***********************************************************************NASR1900 C NASR C JUMP TO HERE TO PRINT-OUT THE SUMMARY OF THE PHASE ADJUSTMENT AFTER NASR1920 C THE LAST PASS NASR1930 C NASR 507 IF(SPEEDM.NE.0.) GO TO 509 NASR 547 CONTINUE 1960 PRINT55551,(NDESC(I),I=1,18) NASR1970 55551 FORMAT('1'//20X,18A4/) NASR1980 WRITE(IOUT,565) NASR 565 FORMAT(45X,'ORBIT-KNOWN MODE'/) NASR PRINT 559 1990 559 FORMAT(//30X,'P H A S E S O L U T I O N - S U M M A R Y ',//) NASR C NASR C PRINT-OUT THE DETAILS OF PROCESSED ,REJECTED ,AND ACCEPTED PASSES. NASR C NASR PRINT87778,NPASS NASR 87778 FORMAT(/5X,'DETAIL INFORMATION CONCERNING THE USED NUMBER OF PASSENASR 1S ARE :'/,5X,'----------------------------------------------------NASR 2-------'//10X,'TOTAL NUMBER OF PASSES BEING PROCESSED IN THIS RUN NASR 3ARE',10X,'=',I5/) NASR PRINT22944,IWRONG,INSING NASR 22944 FORMAT(10X,'NUMBER OF REJECTED PASSES ,DUE TO INCONSISTANT STATIONNASR 1 NUMBERS =',I5,/10X,'NUMBER OF REJECTED PASSES ,WHICH DO NOT HAVENASR 2 ANY SOLUTION ARE =',I5) NASR PRINT78887,IDOPL2,ICHI95,KPASS NASR 78887 FORMAT(10X,'NUMBER OF REJECTED PASSES ,WHICH HAVE LESS THAN 2 DOPPNASR 1LERS ARE =',I5/,10X,'NUMBER OF REJECTED PASSES ,AT 95% PROBABILITNASR 2Y-LEVEL ARE',9X,'=',I5//,10X,'TOTAL NUMBER OF ACCEPTED PASSES AFTENASR 3R THE ABOVE REJECTIONS ARE =',I5//) NASR IF(KPASS.EQ.0) GO TO 544 NASR C NASR C PRINT THE PARAMETERS AND TRANSLATIONS OF THE GIVEN REFERENCE DATUM. NASR C NASR PRINT 11223 NASR PRINT 32211,A1,AM,FLAT,XTRANS,YTRANS,ZTRANS NASR C NASR C COMPUTATION OF THE ADJUSTED GEOD. CURVILINEAR COORD. (LAT.,LONG.,HT.)NASR2070 C OF THE RECEIVER POSITION FROM THE FINAL ESTIMATE OF (XYZ) GEOCENTRIC NASR2080 C CARTESIAN COORDINATES. NASR C NASR2090 JROW=0 2100 DO 509 IN=1,NGRPTS 2110 FOFF=FDI(IN)/60. 2120 IF(NTEST(IN).LT.10) GO TO 509 2150 C NASR2160 C COMPUTE ADJUSTED (XYZ) FOR THE RECEIVER POSITION REFERRED TO THE NASR2170 C USED GEODETIC DATUM NASR2180 C NASR2190 XRNAD=X(IN)-XTRANS NASR2200 YRNAD=Y(IN)-YTRANS NASR2210 ZRNAD=Z(IN)-ZTRANS NASR2220 C NASR2230 C COMPUTE (LAT.,LONG.,HT.) FROM (XYZ) - W.R.T. THE USED GEODETIC DATUM NASR2240 C NASR2250 CALL CARGEO(XRNAD,YRNAD,ZRNAD,A1,ESQ,BLAT,BLOG,H,RAN) NASR2260 C NASR2270 C COMPUTE POSITION NCOL MIN. AFTER LOCK-ON TIME FOR THE MOVING SHIP NASR2280 C NASR INISH=0 2290 445 DAN=NCOL*SPEED(IN)*1852./(60.*RAN) 2300 DIRCT=DIR(IN)*RADG 2310 DLAT= BLAT+DAN*COS(DIRCT)*(1.-ESQ*DSIN(BLAT))/(1.-ESQ) 2320 DLOG= BLOG+DAN*SIN(DIRCT)/(DCOS(BLAT)) 2330 IF(INISH.EQ.1) CALL DGEOCA(X(IN),Y(IN),Z(IN),A1,ESQ,DLAT,DLOG,H) 2340 DLOG=DLOG/RADG 2350 DLAT=DLAT/RADG 2360 SGN=DLAT/DABS(DLAT) 2370 M1=DLAT 2380 N1=(DLAT-M1)*60.*SGN 2390 S1=((DLAT-M1)*60.*SGN-N1)*60. 2400 SGN=DLOG/DABS(DLOG) 2410 M2=DLOG 2420 N2=(DLOG-M2)*60.*SGN 2430 S2=((DLOG-M2)*60.*SGN-N2)*60. 2440 C NASR C PRINT THE FINAL ESTIMATE OF BOTH CARTESIAN AND CURVILINEAR COORDINAT-NASR C ES OF THE TRACKING STATION 'ICHOOZ'. NASR C NASR 5055 PRINT5060 NASR2560 5060 FORMAT(///5X,'***FINAL ADJUSTED COORDINATES OF THE RECEIVER POSITINASR2570 1ON BOTH GEODETIC (CURVILINEAR & CARTESIAN)'/20X,'AND GEOCENTRIC (CNASR2580 2ARTESIAN ONLY)***') NASR2590 PRINT5082 NASR2600 5082 FORMAT(/5X,'ST.NO.',7X,'LAT',18X,'LONG',14X,'HT',13X,'X',15X,'Y',1NASR2610 15X,'Z') NASR2620 PRINT 5079,IDENT(1),M1,N1,S1,M2,N2,S2,H,XRNAD,YRNAD,ZRNAD NASR 5079 FORMAT('0',4X,I4,3X,2(2I4,F8.3,5X),F8.3,5X,3(F12.2,3X)) NASR2640 PRINT5059,X(IN),Y(IN),Z(IN) NASR2650 5059 FORMAT(/23X,'**CORRESPONDING GEOCENTRIC COORDINATES** =',2X,3(F12.NASR2660 12,3X)) NASR2670 PRINT55552 NASR2680 55552 FORMAT(//50X,'==================='//) NASR2690 PRINT 45654 NASR 45654 FORMAT(08X,'NOTE : THE UNITS OF THE QUANTITIES GIVEN IN THE ABOVE NASR 1TABLE ARE AS FOLLOWS :-'/,15X,'FOR LAT : DEG ,MIN & SEC -(+VE NORNASR 2TH)'/,15X,'FOR LONG : DEG ,MIN & SEC -(+VE EAST)'/,15X,'FOR HT ,X NASR 3,Y & Z : METERS') NASR 5056 CONTINUE NASR2700 C NASR C PRINT THE VARIANCE-COVARIANCE MATRIX AND THE STANDARD DEVIATIONS OF NASR C THE RECEIVER'S (X,Y & Z) GEOCENTRIC COORDINATES. NASR C NASR PRINT 79821 NASR 79821 FORMAT('1'//) NASR PRINT 559 NASR PRINT 79822,MDFM NASR 79822 FORMAT(5X,'TOTAL ACCUMULATED NO. OF DEGREES OF FREEDOM =',I5//) NASR PRINT 66619,SIGM NASR PRINT79823 NASR 79823 FORMAT(5X,'VARIANCE-COVARIANCE MATRIX OF THE ADJUSTED (X,Y,Z) RECENASR 1IVER COORD.'/) NASR DO 79824 I=1,3 NASR 79824 PRINT150,(ANORM(I,J),J=1,3) NASR PRINT 6100,SIGX,SIGY,SIGZ NASR C NASR C PRINT-OUT THE THREE SEMI-AXES OF THE STANDARD CONFIDENCE REGION (OR NASR C THE STANDARD ERROR ELLIPSOID) OF THE FINAL ESTIMATED COORDINATES OF NASR C THE ICHOOZ-STATION ,AFTER THE FINAL PASS -IN METERS. NASR C NASR SEM1=DSQRT(EX(1,1)) NASR SEM2=DSQRT(EX(2,2)) NASR SEM3=DSQRT(EX(3,3)) NASR PRINT 55799 NASR 55799 FORMAT(//,5X,'SEMI-AXES OF THE STANDARD ERROR ELLIPSOID'/5X,'PRINTNASR 1ED IN METERS -AFTER THE LAST ACCEPTED PASS'/) NASR PRINT 44688,SEM1,SEM2,SEM3 NASR 44688 FORMAT(13X,F15.3/13X,F15.3/13X,F15.3) NASR IF(INISH.EQ.1) GO TO 576 2710 C NASR C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*NASR C NASR C COMPUTE AND PRINT THE VARIANCE-COVARIANCE MATRIX ,AND THE STANDARD NASR C DEVIATIONS OF THE CURVILINEAR COORDINATES (LAT,LONG,HT) OF THE TRAC-NASR C KING STATION ,FROM THE CORRESPONDING ONES OF (X,Y,Z) COORD. NASR C NASR C COMPUTE TRANSPOSE OF THE ROTATION MATRIX FOR TRANSFORMATION OF COORD,NASR2730 C FROM THE GEOCENTRIC SYSTEM TO THE TOPOCENTRIC SYSTEM NASR2740 C 2750 CR(2)=PI-DLOG*RADG 2760 CR(1)=PI/2.-DLAT*RADG 2770 LR(1)=2 2780 LR(2)=3 2790 CALL DROTAT(CR,LR,2,ROT) 2800 JK=3 2810 IF(NTEST(IN).EQ.12)JK=JK+2 2820 IF(NTEST(IN).EQ.13)JK=JK+1 2830 IF(NTEST(IN).EQ.14)JK=JK+1 2840 DO 597 I=1,6 2850 DO 597 J=1,6 2860 DAR(I,J)=0. 2870 VAR(I,J)=0. 2880 597 AROT(I,J)=0. 2890 JT=MACOL-NGRPTS+IN 2900 JL=JK+1 2910 DO 540 I=1,JK 2920 VAR(I,JL)=ANORM(JROW+I,JT) 2930 VAR(JL,I)=VAR(I,JL) 2940 DO 541 J=1,JK 2950 541 VAR(I,J)=ANORM(JROW+I,JROW+J) 2960 COUNT(I+1)=0. 2970 AROT(I+1,I+1)=1. 2980 540 CONTINUE 2990 VAR(JL,JL)=ANORM(JT,JT) 3000 DO 542 I=1,3 3010 DO 542 J=1,3 3020 542 AROT(I,J)=ROT(I,J) 3030 C NASR3040 C ROTATING THE VARIANCES AND COVARIANCES FROM (XYZ) TO LAT,LONG,HT) NASR3050 C NASR3060 C COMPUTE VAR=AROT'*VAR*AROT ,USING THE SUBROUTINE 'MULTP' TO PERFORM NASR C THE NECESSARY MULTIPLICATIONS. NASR C 3090 CALL MULTP(JL,JL,JL,6,6,VAR,AROT,DAR) 3100 CALL MULTP(JL,JL,JL,6,6,DAR,AROT,VAR) 3110 IDFM=MDFM 3120 C NASR C PRINT THE VAR.-COV. MATRIX OF THE ESTIMAED (LAT,LONG,HT) RECEIVER NASR C COORDINATES. NASR C NASR 576 PRINT 1561 NASR 1561 FORMAT(//,5X,'VARIANCE-COVARIANCE MATRIX OF THE ESTIMATED (LAT,LONNASR 1G,HT) RECEIVER COORD.'/) NASR JL=3 NASR DO 543 I=1,JL 3160 543 WRITE(IOUT,150) (VAR(I,J),J=1,JL) 3180 C 150 FORMAT('0',3(3X,3E13.5)) 3190 C NASR3200 C COMPUTE THE STANDARD DEVIATION FOR (LAT.,LONG.,HT.) IN METER UNITS NASR3210 C NASR3220 SIGLAT=DSQRT(VAR(1,1)) 3230 SIGLOG=DSQRT(VAR(2,2)) 3240 SIGHT=DSQRT(VAR(3,3)) 3250 C NASR C PRINT THE STANDARD DEVIATIONS OF THE LAT ,LONG & HT IN METERS. NASR C NASR PRINT 6173,SIGLAT,SIGLOG,SIGHT 3260 6173 FORMAT(//5X,'STANDARD DEVIATIONS OF (LAT ,LONG & HT.) COORD. - (UNNASR 1ITS ARE IN METERS)'//8X,'FOR LAT.',F12.3/8X,'FOR LONG.',F11.3/8X,'NASR 2FOR HEIFHT',F10.3) NASR C 3350 C PUNCH RECEIVER CO-ORDINATES (X,Y,Z) FOR EVERY INTEIG PASS IF REQUIREDNASR C (WHEN THE OPTION 'IPUNCH' = 1 OR 2 ) NASR C 3370 IF(IPUNCH.EQ.0)GO TO 8877 3380 DO 8878 IWW=1,KWW 3390 8878 PUNCH8879,(STORE(IWW,IWP),IWP=1,3) 3400 8879 FORMAT(3X,3(F12.2,3X)) 3410 8877 CONTINUE 3420 C 3430 IF(INISH.EQ.1)GO TO 544 3440 IF(SPEED(IN).EQ.0) GO TO 544 3450 C 3460 C***********************************************************************NASR C NASR C NOTE THAT THE FOLLOWING SEGMENT OF THIS PROGRAM IS DISREGARDED FOR NASR C OUR CASE OF PHASED-SOLUTION OF A SINGLE STATIONARY TRACKING STATION NASR C NASR C FOR SPEED NE 0 CONSTRAIN HEIGHT AT GIVEVEN H=H0 3480 C 3490 IF(MDFM.EQ.0)SIGM=1. 3500 C CONDITION UPON ADJUSTED PARAMETERS 3510 C PRINT 203,(NDESC(I),I=1,18),NPASS 3520 C 203 FORMAT('1',/,20X,18A4,10X,'PASS NO. ',I5,/) 3530 PRINT 587 3540 587 FORMAT('0',///,' CONSTRAINED HEIGHT SOLUTION ',///) 3550 CNST=VAR(3,3) 3560 C CONTRIBUTION TO VPV & MDFM 3570 IDFM=MDFM+1 3580 VPV=VPV+(H-H0)**2*SIGM**2/CNST 3590 C PRINT 161,VPV,VPV 3600 C 161 FORMAT('0',' VPV=',F15.3,' VPV1=',F15.3,' DVPV=',F15.3,//) 3610 IF(NROW.EQ.1)CNST=1. 3620 IF(NROW.EQ.1)SIGM=1. 3630 DSIGM=VPV/IDFM 3640 DO 546 I=1,JL 3650 COUNT(I)=-VAR(I,3)*(H-H0)/CNST 3660 DO 546 J=1,JL 3670 546 DAR(I,J)=(VAR(I,J)-VAR(I,3)*VAR(3,J)/CNST)*DSIGM/SIGM**2 3680 BLAT=BLAT+COUNT(1)*(1.-ESQ*DSIN(BLAT))/(RAN*(1.-ESQ)) 3690 BLOG=BLOG-COUNT(2)/(RAN*DCOS(BLAT)) 3700 H=H+COUNT(3) 3710 JT=4 3720 IF(NTEST(IN).LT.12) GO TO 548 3730 IF(NTEST(IN).EQ.14) GO TO 577 3740 SPEED(IN)=SPEED(IN)+COUNT(JT) 3750 JT=JT+1 3760 IF(NTEST(IN).EQ.13) GO TO 548 3770 577 DIR(IN)=DIR(IN)+COUNT(JT) 3780 JT=JT+1 3790 C FREQUENCY OFFSET UNKNOWN 3800 548 FOFF=FOFF+COUNT(JT)/600.0000 3810 WRITE(IOUT,532)IN,FOFF 3820 532 FORMAT('0',' STATION NO=',I4,' FREQ. OFFSET=',F15.4) 3830 INISH=1 3840 DO 589 I=1,JL 3850 DO 589 J=1,JL 3860 589 VAR(I,J)=DAR(I,J) 3870 DSIGM=DSQRT(DSIGM) 3880 PRINT 562,DSIGM 3890 C 562 FORMAT('0',' ESTIMATED VARIANCE FACTOR =',F15.3) 3900 GO TO 445 3910 C 3920 C***********************************************************************NASR C NASR 544 PRINT 578 NASR 578 FORMAT('1') NASR 509 CONTINUE NASR C NASR C***********************************************************************NASR C NASR STOP NASR END NASR C***********************************************************************NASR C* *NASR C* FUNCTION SUBPROGRAM 'FREAD' : *NASR C* --------------------------- *NASR C* THE PURPOSE OF THIS FUNCTION IS TO READ THE *NASR C* INPUTTED NUMBERS ON CARDS IN A 'F R E E'-FORMAT ,WHICH NEEDS AT *NASR C* LEAST ONE BLANK BETWEEN ANY TWO CONSECUTIVE NUMBERS. HOWEVER *NASR C* FREAD CANNOT READ ANY INPUTTED CHARACTERS. *NASR C* *NASR C* USAGE : A=FREAD(IN,IOUT) *NASR C* *NASR C* WHERE : A -IS THE VARIABLE NAME THAT REQUIRED TO READ-IN AND *NASR C* STORE ITS VALUE. *NASR C* IN -IS THE INPUT PARAMETERS FOR READING AND PRINTING DATA *NASR C* IN=1 : MEANS READ AND STOR THE INPUTTED NUMBERS ONLY *NASR C* IN=2 : MEANS READ-IN ,STORE ,AND PRINT-OUT THE INPUT *NASR C* DATA AS IT IS ,WITHOUT ANY CHANGES. *NASR C* AND IOUT-IS THE OUTPUT PARAMETER ,IT CAN BE ANY INTEGER VARIA- *NASR C* BLE NAME (I.E. STARTING WITH EITHER I,J,K,L,M OR N) *NASR C* *NASR C* NOTE :- FOR MORE DETAILS SEE THE REFERENCE BELOW. *NASR C* *NASR C***********************************************************************NASR C*********************************************************************** 8340 C* * 8350 C* FUNCTION FOR READING FROM CARDS FORMAT FREE * 8360 C* * 8370 C* R.L.PARKER (SEE BEDFORD INSTITUTE COMPUTER * 8380 C* NOTE 1969-8-C FOR DETAILS) * 8390 C* * 8400 C* MODIFICATIONS BY D.E.WELLS UNB JAN 1971 * 8410 C* * 8420 C*********************************************************************** 8430 C* * 8440 C**********PARAMETER DESCRIPTIONS********** * 8450 C* * 8460 C* INPUT IN =0 NO PRINT * 8470 C* =1 PRINT MESSAGE AND BAD NUMBERS ONLY * 8480 C* =2 PRINT ALL * 8490 C* =3 REREAD LAST NUMBER AND PRINT * 8500 C* * 8510 C* OUTPUT IOUT=0 IF SUCCESSFUL READ * 8520 C* = FIRST BAD CHARACTER IF ERROR HIT * 8530 C* * 8540 C* COLUMN COUNTER ICOL AVAILABLE AS ARGUMENT * 8550 C* USUALLY LEFT UNTOUCHED BY USER * 8560 C* HOWEVER FOR FLEXIBILITY IT CAN BE SET WHEN FUNCTION CALLED * 8570 C* ICOL=81 WILL SKIP TO END OF CURRENT CARD * 8580 C* ICOL.LT.81 WILL SKIP BACK OR FORWARD TO THAT COLUMN * 8590 C* * 8600 C**********POSSIBLE MODIFICATIONS********** * 8610 C* * 8620 C* IF 'IN' NOT REQUIRED - * 8630 C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * 8640 C* CALLING STATEMENTS * 8650 C* INITIALIZE TO DESIRED VALUE ( 0 1 2 OR 3) USING EITHER * 8660 C* DATA STATEMENT ( E.G. 'DATA IN/1/ ' ) OR * 8670 C* REPLACEMENT STATEMENT (E.G. 'IN=2 ' ) * 8680 C* * 8690 C* IF ' IOUT' NOT REQUIRED - * 8700 C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * 8710 C* CALLING STATEMENTS * 8720 C* * 8730 C* IF 'ICOL' NOT REQUIRED - * 8740 C* DROP FROM PARAMETER LISTS IN FUNCTION STATEMENT AND IN ALL * 8750 C* CALLING STATEMENTS * 8760 C* INITIALIZE TO 81 USING DATA STATEMENT ( DATA ICOL/81/ ) * 8770 C* * 8780 C* IF DOUBLE PRECISION DATA EXPECTED (MORE THAN SEVEN DIGITS) - * 8790 C* ADD TO BOTH FUNCTION AND TO CALLING PROGRAM THE STATEMENT * 8800 C* ' DOUBLE PRECISION FREAD ' * 8810 C* * C*********************************************************************** C 8820 FUNCTION FREAD(IN,IOUT) 8830 DOUBLE PRECISION X,XX,B 8840 DOUBLE PRECISION FREAD 8850 DIMENSION IA(80),IDIGIT(17) 8860 DATA IDIGIT/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, 8870 1 1H ,1H, ,1H.,1HE,1HD,1H-,1H+/ 8880 DATA IE/3/ 8890 DATA ICOL/81/ 8900 C*****TEST FOR REREAD MODE (IN=3) 8910 IF(IN.GE.3) GOTO(400,410,860),IE 8920 IOUT=0 8930 C*****TEST FOR EMPTY CARD BUFFER. READ NEW CARD IF EMPTY 8940 IF(ICOL.GT.80) GOTO 201 8950 C** IE IS ONE IF NO E YET HIT 8960 C IC1 ADDRESSES FIRST CHARACTER IN FREAD CALL 8970 C B IS POWER OF TEN BY WHICH FRACTIONAL PART MUST BE DIVIDED 8980 100 IE=1 8990 IC1=ICOL 9000 B=1.0 9010 C*****RETURN HERE TO PICK UP EXPONENT 9020 C IDEC IS ONE IF NO DECIMAL OR E YET HIT 9030 C X IS NUMBER EVALUATED (FREAD) 9040 C SIGN IS SIGN OF FREAD 9050 C NC IS LAST ELEMENT ON IDIGIT LIST 9060 110 IDEC=1 9070 X=0.0 9080 SIGN=+1.0 9090 NC=17 9100 C** IGNORE LEADING SPACES AND COMMAS 9110 DO 150 IC=ICOL,80 9120 IF(IA(IC).NE.IDIGIT(11).AND.IA(IC).NE.IDIGIT(12)) GOTO 202 9130 150 CONTINUE 9140 C*****END OF CARD. ERROR IF EXPONENTIAL PART EXPECTED. 9150 ICOL=ICOL-1 9160 IF(IE.GT.1) GOTO 800 9170 C** NO ERROR. READ ANOTHER CARD 9180 201 READ 2000,(IA(J),J=1,80) 9190 2000 FORMAT(80A1) 9200 ICOL=1 9210 GO TO 100 9220 C **************************************************************** 9230 C NO MORE LEADING SPACES AND COMMAS 9240 C START TESTING CHARACTER BY CHARACTER UNTIL TERMINATING SPACE 9250 C OR COMMA HIT 9260 202 DO 290 ICOL=IC,80 9270 JA=IA(ICOL) 9280 C COMPARE EACH CHARACTER TO IDIGIT LIST 9290 DO 205 I=1,NC 9300 IF(JA.EQ.IDIGIT(I)) GOTO 206 9310 205 CONTINUE 9320 C CHARACTER NOT ON LIST. ERROR 9330 GO TO 800 9340 C CHARACTER ON LIST. TEST WHETHER NUMERAL 9350 206 J=I-9 9360 IF(J.LE.0) GOTO 25 9370 GOTO(25,300,300,210,310,310,220,290),J 9380 C CHARACTER IS DECIMAL. ERROR IF PREVIOUS DECIMAL OR E HIT 9390 210 IDEC=IDEC+IE 9400 IF(IDEC.GT.2) GOTO 800 9410 GOTO 290 9420 C CHARACTER IS MINUS. RESET SIGN AND SHORTEN IDIGIT LIST 9430 220 SIGN=-1.0 9440 GO TO 290 9450 C CHARACTER IS NUMERAL 9460 25 DIGIT=I-1 9470 GOTO(251,252,800),IDEC 9480 C INTEGRAL PART OF NUMBER 9490 251 X=10.0D0*X+DIGIT 9500 GO TO 290 9510 C FRACTIONAL PART OF NUMBER 9520 252 B=B/10.0D0 9530 X=X+B*DIGIT 9540 C CHARACTER IS PLUS. SHORTEN IDIGIT LIST 9550 290 NC=15 9560 C ***************************************************************** 9570 C NO TERMINATING SPACE OR COMMA BEFORE END OF CARD HIT 9580 ICOL=81 9590 C CHARACTER IS TERMINATING SPACE OR COMMA. TEST FOR E OR F FORMAT 9600 300 GO TO(400,410),IE 9610 C CHARACTER IS E. SET XX=MANTISSA 9620 310 IE=IE+1 9630 XX=X*SIGN 9640 C ERROR IF NO MANTISSA 9650 IF(ICOL.EQ.IC) GOTO 800 9660 ICOL=ICOL+3-IE 9670 GO TO (110,110,800),IE 9680 C TERMINATION OF F FORMAT 9690 400 FREAD=X*SIGN 9700 GO TO 500 9710 C TERMINATION OF E FORMAT 9720 410 IX=X*SIGN 9730 FREAD=XX*10.0D0**IX 9740 C** PRINT IF IN=2 OR 3 9750 500 IC2=ICOL-1 9760 IF (IN .GT.1) PRINT 5000,(IA(J),J=IC1,IC2) 9770 5000 FORMAT (10X,80A1) 9780 RETURN 9790 C ***************************************************************** 9800 C ERROR HANDLING 9810 800 IOUT=IA(ICOL) 9820 C SEARCH FOR NEXT TERMINATING BLANK OR COMMA 9830 DO 810 I=ICOL,80 9840 IF(IA(I).EQ.IDIGIT(11).OR.IA(I).EQ.IDIGIT(12)) GO TO 815 9850 810 CONTINUE 9860 I=81 9870 815 ICOL=I 9880 IF(IN.EQ.0) GOTO 850 9890 IC2=ICOL-1 9900 PRINT 820,IOUT,(IA(J),J=IC1,IC2) 9910 820 FORMAT(25H MALFORMED NUMBER DUE TO ,A1,4H IN ,80A1) 9920 850 GO TO(400,410,400),IE 9930 860 PRINT 870 9940 870 FORMAT(25H ILLEGAL REREAD ATTEMPTED ) 9950 FREAD=0. 9960 RETURN 9970 END 9980 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'DGEOCA' : *NASR C* ------------------- *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE *NASR C* THE (X,Y,Z) GEODETIC CARTESIAN COORDINATES OF A SINGLE STATION AS *NASR C* REFERRED TO THE GEOMETRICAL CENTRE OF THE GIVEN REFERENCE ELLIPSOID*NASR C* AS THE ORIGIN OF THE COORDINATE SYSTEM ; FROM ITS GIVEN GEODETIC *NASR C* CURVILINEAR COORDINATES (LAT,LONG,HT) ,WHICH ARE REFERRING TO THE *NASR C* SAME REFERENCE ELLIPSOID. *NASR C* *NASR C* USAGE : CALL DGEOCA(X,Y,Z,A,ESQ,ALAT,ALONG,H) *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- *NASR C* *NASR C* 1-A : SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID -IN METERS. *NASR C* 2-ESQ : SQUARE OF THE FIRST ECCENTRICITY OF THE SAME ELLIPSOID. *NASR C* 3-ALAT: GEODETIC LATITUDE OF THE STATION -IN RADIANS. *NASR C* 4-ALONG:GEODETIC LONGITUDE OF THE STATION -IN RADIANS. *NASR C* 5-H : GEOMETRICAL HEIGHT OF THE STATION ABOVE THE ELLIPSOID IN *NASR C* METERS. *NASR C* *NASR C* THE OUTPUT ,WHICH WILL BE RETURNED BACK TO THE CALLING ROUTINE WILL*NASR C* BE THE (X,Y & Z) GEODETIC CARTESIAN COORDINATES OF THE REQUIRED *NASR C* STATION IN METERS. *NASR C* *NASR C* NOTE :- THIS SUBROUTINE IS VALID FOR ONLY ONE STATION AT A TIME. *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE DGEOCA(X,Y,Z,A,ESQ,ALAT,ALONG,H) 4040 IMPLICIT REAL *8(A-H,O-Z) 4100 C NASR C COMPUTE THE RADIUS OF CURVATURE (RAD) IN THE PRIME VERTICAL PLANE NASR C OF THE STATION. NASR C NASR DENSQ=1.-ESQ*DSIN(ALAT)**2 4110 DEN=DSQRT(DENSQ) 4120 RAD=A/DEN 4130 C NASR C COMPUTE THE (X,Y & Z) COORD. FROM (LAT,LONG & HT) COORD. NASR C NASR X=(RAD+H)*DCOS(ALAT)*DCOS(ALONG) 4140 Y=(RAD+H)*DCOS(ALAT)*DSIN(ALONG) 4150 Z=(RAD*(1.-ESQ)+H)*DSIN(ALAT) 4160 RETURN 4170 END 4180 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'CARGEO' : *NASR C* ------------------- *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE *NASR C* THE GEODETIC CURVILINEAL COORDINATES (LAT,LONG & HT) FOR A GIVEN *NASR C* STATION ,AS REFERRED TO A CHOSEN REFERENCE ELLIPSOID ; FROM THE *NASR C* INPUTTED (X,Y & Z) GEODETIC CARTESIAN COORDINATES OF THE INVOLVED *NASR C* STATION ,WHICH ARE REFERRING TO THE CENTRE OF THE GIVEN REFERENCE *NASR C* ELLIPSOID AS THE ORIGIN OF THE COORDINATE SYSTEM. *NASR C* *NASR C* USAGE : CALL CARGEO(X,Y,Z,A,ESQ,ALAT,ALOG,H,RAN) *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- *NASR C* *NASR C* 1-THE (X,Y,Z) GEODETIC CARTESIAN COORDINATES OF THE STATION - ALL *NASR C* ARE IN METERS. *NASR C* 2-A : THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID - IN METERS. *NASR C* 3-ESQ : THE SQUARE OF THE FIRST ECCENTRICITY OF THE ELLIPSOID. *NASR C* *NASR C* THE OUTPUT WILL BE ; THE GEODETIC LATITUDE (ALAT) AND THE GEODETIC *NASR C* LONGITUDE (ALOG) ,BOTH ARE IN RADIANS.- PLUS THE GEOMETRICAL HEIGHT*NASR C* (H) OF THE STATION ABOVE THE REFERENCE ELLIPSOID ,ALONG WITH THE *NASR C* RADIUS OF CURVATURE (RAN) IN THE PRIME VERTICAL PLANE OF THIS STA- *NASR C* TION ,BOTH ARE IN METERS. *NASR C* *NASR C* NOTE :- THIS SUBROUTINE IS VALID FOR ONE STATION ONLY AT A TIME *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE CARGEO(X,Y,Z,A,ESQ,ALAT,ALOG,H,RAN) NASR8120 IMPLICIT REAL*8(A-H,O-Z) 8130 PI=3.1415926535 8140 D=Z/(DSQRT(X**2+Y**2)) C NASR C COMPUTE AN INITIAL APPROXIMATE VALUE FOR THE LATITUDE (ALAT). NASR C NASR ALAT=DATAN(D) 8160 C NASR C STARTING THE ITERATIVE SOLUTION FOR THE LATITUDE (ALAT) ,AND THE NASR C RADIUS OF CURVATURE IN THE PRIME VERTICAL PLANE (RAN). NASR C NASR 100 BLAT=ALAT 8170 B=DSIN(ALAT) 8180 RAN=A/DSQRT(1.-ESQ*B**2) 8190 ALAT=D*(1.+ESQ*RAN*B/Z) 8200 ALAT=DATAN(ALAT) 8210 C NASR C TESTING THE DIFFERENCE BETWEEN EACH TWO SUCCESSIVE VALUES FOR THE NASR C LATITUDE ; TO STOP THE ITERATION PROCESS ,WHEN THE PRECISION LIMIT NASR C IS REACHED. NASR C NASR IF(DABS(ALAT-BLAT).LE.1.D-11) GO TO 102 8220 GO TO 100 8230 C NASR C COMPUTE THE LONGITUDE (ALOG) DIRECTLY ,IT DOES NOT NEED ITERATIONS. NASR C NASR 102 ALOG=DATAN(Y/X) 8240 IF(ALOG)103,104,104 8250 103 IF(Y.LT.0.) ALOG=2.*PI+ALOG 8260 IF(Y.GT.0.)ALOG=PI+ALOG 8270 GO TO 105 8280 104 IF(Y.LE.0.)ALOG=PI+ALOG 8290 C NASR C COMPUTE THE FINAL VALUE OF THE HEIGHT (H) ,USING THE FINAL VALUES NASR C OBTAINED FOR (ALAT ,AND RAN). NASR C NASR 105 H=DSQRT(X**2+Y**2)/DCOS(ALAT)-RAN 8300 IF(ALOG.GT.PI) ALOG=ALOG-2.*PI NASR8310 RETURN 8320 END 8330 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'DROTAT' : *NASR C* ------------------- *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE *NASR C* THE (3,3) PRODUCT ROTATION MATRIX (ROT) ,WHICH RESULTS FROM A SERI-*NASR C* ES OF UP TO NINE CONSECUTIVE SINGLE AXIS ROTATIONS ,REQUIRED FOR *NASR C* THE COORDINATE TRANSFORMATIONS BETWEEN THE DIFFERENT SYSTEMS OF *NASR C* COORDINATES *NASR C* *NASR C* USAGE : CALL DROTAT(ANGLE,NAXIS,NUM,ROT) *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- *NASR C* *NASR C* 1-ANGLE : IS A VECTOR OF UP TO 9 CONSECUTIVE ANGLES OF SINGLE-AXIS *NASR C* ROTATIONS - IN RADIANS. *NASR C* 2-NAXIS : IS THE VECTOR OF THE AXIS-OF-ROTATION NUMBERS (EITHER 1, *NASR C* 2 OR 3 FOR EACH GIVEN ANGLE) ,I.E. SAME SIZE AS ANGLE *NASR C* 3-NUM : IS THE NUMBER OF REQUIRED ROTATIONS ,WHICH REPRESENTS *NASR C* THE ACTUAL SIZE OF BOTH (ANGLE & NAXIS) VECTORS. *NASR C* *NASR C* THE OUTPUT WILL BE THE (3,3) PRODUCT ROTATION MATRIX (ROT) ,COMPUT-*NASR C* ED IN DOUBLE PRECISION *NASR C* *NASR C* NOTE :- THIS SUBROUTINE IS VALID ONLY UP TO 9 CONSECUTIVE SINGLE *NASR C* ROTATIONS. *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE DROTAT(ANGLE,NAXIS,NUM,ROT) 4190 IMPLICIT REAL*8(A-H,O-Z) 4270 DIMENSION ROT(3,3),R1(3,3),R2(3,3),ANGLE(NUM),NAXIS(NUM) 4280 DO 100 K=1,NUM 4300 C C FIND OTHER TWO AXES N2,N3 4310 C N2=NAXIS(K)+1 4320 IF(N2.GT.3) N2=N2-3 4330 N3=N2+1 4340 IF(N3.GT.3) N3=N3-3 4350 C C SET DIAGONAL ELEMENTS 4360 C R1(NAXIS(K),NAXIS(K))=1. 4370 R1(N2,N2)=DCOS(ANGLE(K)) 4380 R1(N3,N3)=R1(N2,N2) 4390 C C SET NONZERO OFF-DIAGONAL ELEMENTS 4400 C R1(N2,N3)=DSIN(ANGLE(K)) 4410 R1(N3,N2)=-R1(N2,N3) 4420 C C SET ZERO OFF-DIAGONAL ELEMENTS 4430 C R1(NAXIS(K),N2)=0. 4440 R1(NAXIS(K),N3)=0. 4450 R1(N2,NAXIS(K))=0. 4460 R1(N3,NAXIS(K))=0. 4470 C C IF FIRST ROTATION, SET ROT=R1 4480 C IF(K.GT.1) GOTO 20 4490 DO 10 I=1,3 4500 DO 10 J=1,3 4510 10 ROT(I,J)=R1(I,J) 4520 GOTO 100 4530 C C IF NOT FIRST ROTATION SET R2=R1*ROT AND ROT=R2 4540 C 20 DO 30 J=1,3 4550 DO 30 I=1,3 4560 R2(I,J)=0. 4570 DO 30 M=1,3 4580 30 R2(I,J)=R2(I,J)+R1(I,M)*ROT(M,J) 4590 DO 40 I=1,3 4600 DO 40 J=1,3 4610 40 ROT(I,J)=R2(I,J) 4620 100 CONTINUE 4630 RETURN 4640 END 4650 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'ADJUST' : *NASR C* -------------------- *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO ESTABLI-*NASR C* SH THE RELATIVE VARIANCE-COVARIANC MATRIX OF THE USED DOPPLERS , *NASR C* EITHER IDENTITY MATRIX OR WITH OFF-DIAGONALS ,AS DESIRED ; THEN TO *NASR C* PERFORM THE PARAMETRIC LEAST SQUARES ADJUSTMENT OF THE DOPPLER SAT-*NASR C* ELLITE OBSERVATIONS ,UTILIZING THE TECHNIQUE OF WEIGHTING THE UNKN-*NASR C* OWN PARAMETERS ,FOR A SINGLE MOTIONLESS TRACKING STATION (ORBIT KN-*NASR C* OWN MODE) *NASR C* *NASR C* USAGE : CALL ADJUST *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- (USING THE SPECIFIED COMMON BLOCK) *NASR C* *NASR C* 1-THE DESIGN MATRIX-A ,AND THE MISCLOSURE VECTOR-W *NASR C* 2-THE WEIGHT MATRIX-PX OF THE ESTIMATED PARAMETERS FROM THE PREVIO-*NASR C* US PHASE ADJUSTMENT. *NASR C* 3-THE VECTOR-KDPSI OF THE SEQUENCE NUMBERS ASSIGNED TO THE OBSERVED*NASR C* DOPPLERS AT THE TRACKING STATION 'ICHOOZ' *NASR C* 4-THE DIMENSIONS OF A-MATRIX ,MROW AND NACOL *NASR C* 5-THE NUMBER OF UNKNOWN GROUND STATIONS (NGRPTS) ,WHICH IS ONE HERE*NASR C* 6-THE OPTIONAL PARAMETER (IDNT) REQUIRED TO FORM THE RELATIVE VAR- *NASR C* COV. MATRIX OF THE OBSERVED DOPPLERS. *NASR C* 7-THE ACCUMULATED QUADRATIC FORM OF THE WEIGHTED RESIDUALS ,FROM *NASR C* THE PREVIOUS PHASE ADJUSTMENT.-(VPV & VPV1) *NASR C* *NASR C* THE OUTPUT WILL BE : *NASR C* ------------------ *NASR C* *NASR C* 1-THE RELATIVE VAR.-COV. MATRIX ,AND THE WEIGHT MATRIX OF THE USED *NASR C* DOPPLERS IN THIS PHASE ONLY. BOTH ARE PRINTED-OUT BY THE SUBROU-*NASR C* TINE 'ADJUST' ITSELF ,UNDER THE SAME NAME (PDOP)-MATRIX. *NASR C* 2-THE DETERMINANT (D1) OF THE NORMAL EQUATIONS MATRIX. *NASR C* 3-THE WEIGHT-COEFFICIENT MATRIX (ANORM) OF THE ADJUSTED PARAMETERS *NASR C* (WHICH IS THE INVERSE OF THE NORMAL EQUATION MATRIX) *NASR C* 4-THE ABSOLUTE TERMS OF THE NORMAL EQUATIONS : U-VECTOR. *NASR C* 5-THE SOLUTION VECTOR OF THE INVOLVED PARAMETERS. - X-VECTOR. *NASR C* 6-THE RESIDUAL VECTOR (CORRECTIONS TO THE USED DOPPLERS) IN THIS *NASR C* PHASE ONLY. - V-VECTOR. *NASR C* 7-THE ACCUMULATED QUADRATIC FORM OF THE WEIGHTED RESIDUALS (VPV) *NASR C* FROM ALL PHASES ,COMPUTED BY TWO DIFFERENT METHODS. -(VPV & VPV1)*NASR C* *NASR C* THE OTHER SUBROUTINES BEING CALLED BY 'ADJUST' ARE : *NASR C* -------------------------------------------------- *NASR C* *NASR C* 1-MULTP : FOR MULTIPLYING THE TRANSPOSE OF A GIVEN MATRIX BY *NASR C* ANOTHER GIVEN MATRIX. *NASR C* 2-ADMTR : FOR ADDITION OF TWO MATRICIES. *NASR C* 3-DINVPV : FOR INVERTING A GIVEN MATRIX *NASR C* *NASR C* NOTE :- THIS SUBROUTINE IS VALID ONLY FOR A SINGLE MOTIONLESS *NASR C* STATION ,I.E. UP TO 4 UNKNOWNS AND 8 OBSERVED DOPPLERS *NASR C* CAN BE ACCEPTED. *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE ADJUST 4780 IMPLICIT REAL*8(A-H,O-Z) 4800 DIMENSION PX(4,4),A(8,4),W(8),U(4),ANORM(4,4),X(4),PA(8,8),V(8), 4820 1 PW(8),XWT(4),KDPSI(7,14),KNON(1,14),VTP(8),PDOP(8,8) 4830 C NASR C SPECIFY A LABELLED COMMON BLOCK /DOPADJ/ ,FOR THE INPUT DATA TO ,AND NASR C THE OUTPUT DATA FROM THE SUBROUTINE 'ADJUST' ,TO SHARE THIS AREA NASR C WITH THE MAIN PROGRAM 'DOPSATS' ; THIS IS TO SIMPLIFY THE CALLING NASR C NASR COMMON /DOPADJ/ D1,ANORM,PX,A,W,X,U,V,VPV,VPV1,MACOL,MROW, NASR 1 NGRPTS,KDPSI,IDNT NASR C NASR C***********************************************************************NASR C 4910 C I N T I A L I Z A T I O N 4920 C ......................... 4930 C 4940 D1=0.D0 4950 DO 303 I=1,14 4960 303 KNON(1,I)=0 4970 DO 304 IDT=1,8 4980 DO 304 JDT=1,8 4990 304 PDOP(IDT,JDT)=0. 5000 C 5010 C F O R M T H E W E I G H T M A T R I X O F T H E D O P P L E R S 5020 C ..................................................................... 5030 C 5040 IF(IDNT.NE.1) GO TO 97 NASR C NASR C FORM THE WEIGHT MATRIX (PDOP) AS AN IDENTITY MATRIX ,WHEN DESIRED. NASR C NASR DO 5701 IDT=1,MROW 5100 5701 PDOP(IDT,IDT)=1. 5110 GO TO 70 NASR C NASR 97 CONTINUE NASR C NASR C FORM THE RELATIVE VARIANCE-COVARIANCE MATRIX OF THE USED DOPPLERS NASR C IN THIS PHASE ,NAMELY : PDOP-MATRIX. NASR C NASR C ESTABLISH PDOP-MATRIX ,WITH ITS OFF-DIAGONAL ELEMENTS ,TO INCORPO- NASR C RATE THE CORRELATION BETWEEN THE ADJACENT DOPPLERS. NASR C NASR C INSERT THE DIAGONAL ELEMENTS OF PDOP-MATRIX (ALL ARE 2'S) NASR C NASR DO 80 ICOR=1,MROW NASR 80 PDOP(ICOR,ICOR)=2. NASR C NASR C INSERT THE OFF-DIAGONAL ELEMENTS OF -1.0 ON BOTH SIDES OF THE MAIN NASR C DIAGONAL OF PDOP-MATRIX. NASR C NASR DO 5702 JDT=2,MROW 5160 5702 PDOP(JDT-1,JDT)=-1.0 NASR DO 5703 IDT=2,MROW 5180 5703 PDOP(IDT,IDT-1)=-1.0 NASR C 5200 C RE-SORT KDPSI KEEPING ONLY THE NON-ZERO DOPPLERS IN KNON 5210 C 5220 JNON=0 5230 DO 5704 I=1,NGRPTS 5240 DO 5704 J=1,8 5250 IF(KDPSI(I,J).EQ.0)GO TO 5704 5260 JNON=JNON+1 5270 KNON(I,JNON)=KDPSI(I,J) 5280 5704 CONTINUE 5290 C 5300 C REMOVE THE OFF-DIAGONALS OF -1.0 ,WHERE ZERO DOPPLER OCCUR ,I.E. NASR C CANCEL THE THE EXISTING CORRELATION BETWEEN THE NON-ADJACENT DOPPLERSNASR C 5320 JNON=MROW-1 5330 DO 5705 I=1,NGRPTS 5340 DO 5705 J=1,JNON 5350 IF(KNON(I,J)+1.EQ.KNON(I,J+1))GO TO 5705 5360 PDOP(J,J+1)=0. 5370 PDOP(J+1,J)=0. 5380 5705 CONTINUE 5390 C 5400 C PRINT THE RELATIVE VARIANCE-COVARIANCE MATRIX OF THE USED DOPPLERS NASR C IN THIS PHASE ONLY NASR C 5420 PRINT50 5430 50 FORMAT(//9X,'RELATIVE VAR.-COV. MATRIX OF THE USED DOPPLERS'/) NASR DO 60 IDT=1,MROW 5460 60 PRINT5023,(PDOP(IDT,JDT),JDT=1,MROW) 5470 C 5480 C INVERT THE RELATIVE VARIANCE-COVARIANCE MATRIX (OR THE WEIGHT COEFF- NASR C ICIENT MATRIX) OF THE USED DOPPLERS ,TO OBTAIN THEIR CORRESPONDING NASR C WEIGHT MATRIX (THE SAME NAME 'PDOP' WILL BE USED). NASR C 5510 CALL DINVPV(PDOP,MROW,DP,8) 5520 C 5530 C PRINT OUT THE WEIGHT MATRIX OF THE OBSERVATIONS 5540 C 5550 70 WRITE(6,154) 5560 154 FORMAT(//17X,'WEIGHT MATRIX OF THE USED DOPPLERS'/) NASR DO 5022 IDT=1,MROW 5580 5022 PRINT5023,(PDOP(IDT,JDT),JDT=1,MROW) 5590 5023 FORMAT(' ',10X,8F8.2) 5600 C 5610 C F O R M T H E N O R M A L E Q U A T I O N S ; A N O R M M A T R I 5620 C ........................................................................5630 C 5640 CALL MULTP(MROW,MROW,MACOL,08,08,PDOP,A,PA) 5650 CALL MULTP(MROW,MACOL,MACOL,08,04,A,PA,ANORM) 5660 C ANORM =A'*PDOP*A 5670 CALL ADMTR(MACOL,MACOL,04,04,PX,ANORM) 5680 C ANORM =PX+A'*PDOP*A 5690 C 5700 C F O R M T H E A B S O L U T E T E R M ; U - V E C T O R 5710 C ........................................................... 5720 C 5730 DO 203 I=1,MROW 5740 PW(I)=0. 5750 DO 203 J=1,MROW 5760 203 PW(I)=PW(I)+PDOP(I,J)*W(J) 5770 DO 103 I=1,MACOL 5780 U(I)=0. 5790 DO 103 J=1,MROW 5800 103 U(I)=U(I)+A(J,I)*PW(J) 5810 C U= A'*PDOP*W 5820 C 5870 C I N V E R T N O R M A L E Q U A T I O N S 5880 C ........................................... 5890 C 5900 112 CALL DINVPV(ANORM,MACOL,D1,04) NASR C NASR C CHECK ; IF THE NORMAL EQUATIONS MATRIX IS SINGULAR ? (I.E. ITS DETER-NASR C MINANT IS ZERO) ; -PRINT A MESSAGE ,STOP EXECUTION ,AND RETURN. NASR C NASR IF(D1.EQ.0.) GO TO 95 NASR C 5950 C C O M P U T E T H E S O L U T I O N V E C T O R ; X - V E C T O R 5960 C ..................................................................... 5970 C 5980 DO 116 I=1,MACOL 5990 X(I)=0. 6000 DO 116 J=1,MACOL 6010 116 X(I)=X(I)-ANORM(I,J)*U(J) 6020 C X=-ANORM-1 U 6030 C 6070 C C O M P U T E T H E R E S I D U A L S ; V - V E C T O R 6080 C ......................................................... 6090 C 6100 150 DO 120 I=1,MROW 6110 V(I)=0. 6120 DO 121 J=1,MACOL 6130 121 V(I)=V(I)+A(I,J)*X(J) 6140 C V = AX 6150 120 V(I)=V(I)+W(I) 6160 C V = AX + W 6170 C 6210 C T H E W E I G H T E D S Q U A R E S U M O F T H E R E S I D U A L 6220 C ........................................................................6230 C NASR C COMPUTE (V'PV) BY THE DIRECT METHOD ; V'PV = V'*PDOP*V + X'*PX*X. NASR C NOTE : THE SECOND TERM (X'PX X) ,IS DUE TO WEIGHTING THE PARAMETERS. NASR C NASR DO 400 J=1,MROW 6250 VTP(J)=0. 6260 DO 400 I=1,MROW 6270 400 VTP(J)=VTP(J)+V(I)*PDOP(I,J) 6280 DO 401 I=1,MROW 6290 401 VPV=VPV+VTP(I)*V(I) 6300 C 6310 XPX=0. 6320 DO 128 I=1,MACOL 6330 XWT(I)=0 NASR6340 DO 135 J=1,MACOL 6350 135 XWT(I)=XWT(I)+X(J)*PX(I,J) NASR6360 128 XPX=XPX+XWT(I)*X(I) NASR6370 VPV=VPV+XPX 6380 C NASR6390 C COMPUTE (V'PV) BY ANOTHER INDIRECT METHOD ,USING THE LAGRANGE MULT- NASR C TIPLIERS 'K' ; V'PV = -K'*W = V'*P*W. NASR C NASR6410 DO 126 I=1,MROW 6420 126 VPV1=VPV1+VTP(I)*W(I) 6430 GO TO 85 NASR 95 CONTINUE NASR PRINT 99 NASR 99 FORMAT(///,5X,'NORMAL EQUATIONS MATRIX IS SINGULAR - NO AVAILABLE NASR 1SOLUTION FOR THIS PASS'/) NASR 85 RETURN NASR END 6450 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'ADMTR' : *NASR C* ------------------ *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO ADD TWO *NASR C* SUBMATRICES A1(NR,NC) ,AND B1(NR,NC) ,WHICH WILL BE EXTRACTED FROM *NASR C* THE UPPER LEFT PARTS OF THE TWO INPUTTED GENERAL MATRICES A(ND1, *NASR C* NC1) ,AND B(ND2,NC2) ; THE RESULTED SUBMATRIX AFTER PERFORMING THE *NASR C* ADDITION ,WILL BE STORED IN THE SAME PLACE OF THE UPPER LEFT PART *NASR C* (NR,NC) ,OF THE ORIGINAL INPUTTED B-MATRIX. *NASR C* *NASR C* USAGE : CALL ADMTR(NR,NC,ND1,ND2,A,B) *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- *NASR C* *NASR C* 1-NR & NC : THE ROW AND COLUMN DIMENSIONS ,OF BOTH SUBMATRICES *NASR C* REQUIRED TO BE EXTRACTED FROM A & B MATRICES AND *NASR C* THEN ADDED,. *NASR C* 2-ND1 & ND2 : THE ROW DIMENSIONS OF THE GIVEN MATRICES A & B *NASR C* 3- A & B : THE TWO ORIGINAL GENERAL GIVEN MATRICES. *NASR C* *NASR C* THE OUTPUT WILL BE THE MATRIX (B) IN ITS GIVEN FULL SIZE ,AFTER *NASR C* REPLACING ITS UPPER LEFT PART (NR,NC) BY THE RESULTING SUBMATRIX *NASR C* FROM THE REQUIRED ADDITION. *NASR C* *NASR C* NOTE : BOTH SUBMATRICES REQUIRED TO BE ADDED SHOULD BE CONFOR- *NASR C* MABLE FOR ADDITION ,I.E. HAVING THE SAME SIZE (NR,NC). *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE ADMTR(NR,NC,ND1,ND2,A,B) 7520 IMPLICIT REAL*8 (A-H,O-Z) 7530 DIMENSION A(1),B(1) 7590 KA=-ND1 7600 KB=-ND2 7610 DO 100 J=1,NC 7620 KA=KA+ND1 7630 KB=KB+ND2 7640 DO 100 I=1,NR 7650 IA=KA+I 7660 IB=KB+I 7670 100 B(IB)=B(IB)+A(IA) 7680 RETURN 7690 END 7700 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'MULTP' : *NASR C* ------------------ *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO MULTIPLY*NASR C* THE TRANSPOSE OF THE SUBMATRIX A1(NRA,NCA) ,BY ANOTHER SUBMATRIX *NASR C* B1(NRA,NCB) ,WHICH WILL BE EXTRACTED FROM THE UPPER LEFT PARTS OF *NASR C* THE TWO INPUTTED GENERAL MATRICES A(ND1,NC1) ,AND B(ND1,NC2) , *NASR C* RESPECTIVELY. THE RESULTING PRODUCT SUBMATRIX RESU(NCA,NCB) WILL *NASR C* BE STORED IN THE SAME PLACE OF THE UPPER LEFT PART OF A DESIRED *NASR C* GENERAL MATRIX RES(ND2,NC3) . IT SHOULD BE NOTED THAT THE ORIGINAL*NASR C* INPUTTED A & B MATRICES SHOULD HAVE THE SAME ROW DIMENSION (ND1) *NASR C* *NASR C* USAGE : CALL MULTP(NRA,NCA,NCB,ND1,ND2,A,B,RES) *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- *NASR C* *NASR C* 1-NRA : THE ROW DIMENSION OF BOTH SUBMATRICES REQUIRED TO BE MULT- *NASR C* IPLIED AFTER TRANSPOSING THE FIRST ONE. *NASR C* 2-NCA : THE COLUMN DIMENSION OF THE FIRST SUBMATRIX ,BEFORE BEING *NASR C* TRANSPOSED. *NASR C* 3-NCB : THE COLUMN DIMENSION OF THE SECOND SUBMATRIX. *NASR C* 4-ND1 : THE ROW DIMENSION OF BOTH ORIGINALLY GIVEN GENERAL *NASR C* MATRICES A & B . *NASR C* 5-ND2 : THE ROW DIMENSION OF THE GENERAL MATRIX (RES) ,DESIRED *NASR C* TO INCLUDE THE PRODUCT SUBMATRIX IN ITS UPPER LEFT PART *NASR C* (NCA , NCB). *NASR C* 6- A : THE FIRST GENERAL MATRIX. *NASR C* 7- B : THE SECOND GENERAL MATRIX. *NASR C* *NASR C* THE OUTPUT WILL BE THE DESIRED MATRIX (RES) ,AFTER REPLACING ITS *NASR C* UPPER LEFT PART (NCA , NCB) BY THE REQUIRED PRODUCT SUBMATRIX. *NASR C* *NASR C* NOTE : THE CHOSEN SUBMATRICES TO BE EXTRACTED FROM THE ORIGINAL *NASR C* A & B MATRICES ,SHOULD HAVE THE SAME ROW DIMENSION (NRA) *NASR C* TO BE CONFORMABLE FOR MULTIPLICATION ,AFTER TRANSPOSING *NASR C* THE FIRST ONE. *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE MULTP(NRA,NCA,NCB,ND1,ND2,A,B,RES) 7710 IMPLICIT REAL*8 (A-H,O-Z) 7720 DIMENSION A(1),B(1),RES(1) 7790 KB=-ND1 7800 KR=-ND2 7810 DO 100 J=1,NCB 7820 KA=-ND1 7830 KB=KB+ND1 7840 KR=KR+ND2 7850 DO 100 I=1,NCA 7860 KA=KA+ND1 7870 LR=KR+I 7880 RES(LR)=0. 7890 DO 100 II=1,NRA 7900 LA=KA+II 7910 LB=KB+II 7920 100 RES(LR)=RES(LR)+A(LA)*B(LB) 7930 RETURN 7940 END 7950 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'DINVPV' : *NASR C* ------------------- *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO RELACE *NASR C* THE UPPER LEFT SUBMATRIX A1(N,N) , FROM A GENERAL MATRIX A(M,M) *NASR C* ; BY ITS INVERSE A1(-1)(N,N) ,IN ITS PLACE IN THE ORIGINAL GIVEN *NASR C* A-MATRIX. USING ''GAUSS-JORDAN'' METHOD OF THE MAXIMUM PIVOTAL *NASR C* ELEMENT. (REFERENCE : THE U.N.B. COMPUTING CENTRE LIBRARY. ) *NASR C* *NASR C* USAGE : CALL DINVPV(A,N,DETER,M) *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- *NASR C* *NASR C* 1-A : THE GENERAL MATRIX OF ORDER (M). *NASR C* 2-N : THE ORDER OF THE UPPER LEFT SUBMATRIX ,REQUIRED TO BE *NASR C* EXTRACTED FROM THE A-MATRIX ,AND THEN INVERTED. *NASR C* 3-M : THE ORDER OF THE GENERAL GIVEN A-MATRIX. *NASR C* *NASR C* NOTE THAT IN CASE OF INVERTING THE WHOLE ORIGINAL MATRIX (A) *NASR C* THE ORDERS M & N WILL BE EQUAL. *NASR C* *NASR C* THE OUTPUT WILL BE THE MATRIX A(M,M) ,AFTER REPLACING ITS UPPER *NASR C* LEFT PART (N,N) BY THE CORRESPONDING REQUIRED INVERSE. THE DETER- *NASR C* MINANT (DETER) OF THE SUBMATRIX A1(N,N) ,BEFORE INVERSION WILL *NASR C* BE OLSO COMPUTED AND GIVEN AS AN OUTPUT. *NASR C* *NASR C* NOTE :- THIS SUBROUTINE IS VALID FOR INVERTING A MATRIX OF *NASR C* ORDER ONLY UP TO ( 100 * 100 ). *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE DINVPV(A,N,DETER,M) 00006460 IMPLICIT REAL*8 (A-H,O-Z) 00006480 DIMENSION A(M , M),Y(100),ICJ(100),IRR(100),JCC(100) 00006490 IF(M.GT.100)GO TO 26 00006500 2 DETER=1.0 00006510 DO 16 K=1,N 00006520 BIGA=0.0 00006530 C 00006540 C LOOP TO FIND MAXIMUM PIVOTAL ELEMENT AVAILABLE 00006550 C DO 7 I=1,N 00006560 DO 7 J=1,N 00006570 IF(K-1)58,58,13 00006580 13 K5=K-1 00006590 C 00006600 C CHECK THAT SAME PIVOTAL ROW NOT USED TWICE. 00006610 C DO 54 I1=1,K5 00006620 IF (I-IRR(I1))54,7,54 00006630 54 CONTINUE 00006640 C 00006650 C CHECK THAT SAME PIVOTAL COLUMN NOT USED TWICE. 00006660 C 55 DO 57 J1=1,K5 00006670 IF (J-JCC(J1))57,7,57 00006680 57 CONTINUE 00006690 C 00006700 C LARGEST PIVOT ELEMENT STORED AS BIGA 00006710 C 58 ZZ=DABS(A(I,J)) 00006720 IF(ZZ-BIGA)7,7,6 00006730 6 BIGA=ZZ 00006740 IRR(K)=I 00006750 JCC(K)=J 00006760 7 CONTINUE 00006770 C 00006780 C DECIDES MATRIX IS NEAR-SINGULAR IF LARGEST AVAILABLE PIVOT IS 00006790 C LESS THAN 0.1E-50 00006800 C IF(BIGA-0.1E-50)8,9,9 00006810 8 PRINT 51 00006820 51 FORMAT(' DINVPV ERROR.. ARRAY IS NEAR SINGULAR') 6830 DETER=0.0 00006840 RETURN 00006850 9 IR=IRR(K) 00006860 JC=JCC(K) 00006870 BIGA=A(IR,JC) 00006880 DETER=DETER*BIGA 00006890 C 00006900 C NORMALIZES ITH PIVOTAL ROW. 00006910 C DO 10 J=1,N 00006920 10 A(IR,J)=A(IR,J)/BIGA 00006930 A(IR,JC)=1.0/BIGA 00006940 C 00006950 C ROW AND COLUMN OPERATIONS FOR CLEARING. 00006960 C DO 16 I=1,N 00006970 AJCK=A(I,JC) 00006980 IF (I-IR)12,16,12 00006990 12 A(I,JC)=-AJCK/BIGA 00007000 DO 15 J=1,N 00007010 IF (J-JC)14,15,14 00007020 14 A(I,J)=A(I,J)-(AJCK*A(IR,J)) 00007030 15 CONTINUE 00007040 16 CONTINUE 00007050 C 00007060 C COUNT ROW AND COLUMN INTERCHANGES NECESSARY TO REORDER 00007070 C DO 17 I=1,N 00007080 IR=IRR(I) 00007090 JC=JCC(I) 00007100 17 ICJ(IR)=JC 00007110 ICOUNT=0 00007120 N5=N-1 00007130 DO 19 I=1,N5 00007140 I5=I+1 00007150 DO 19 J=I5,N 00007160 IF(ICJ(J)-ICJ(I))18,19,19 00007170 18 TEMP=ICJ(J) 00007180 ICJ(J)=ICJ(I) 00007190 ICJ(I)=TEMP 00007200 ICOUNT=ICOUNT+1 00007210 19 CONTINUE 00007220 C 00007230 C DETERMINE SIGN OF DETERMINENT. 00007240 C IF(ICOUNT/2*2-ICOUNT)20,21,21 00007250 20 DETER=-DETER 00007260 C 00007270 C COLUMNS STORED IN Y AS ARRAY IS REARRANGED IN PLACE 00007280 C 21 DO 23 J=1,N 00007290 DO 22 I=1,N 00007300 IR=IRR(I) 00007310 JC=JCC(I) 00007320 22 Y(JC)=A(IR,J) 00007330 DO 23 I=1,N 00007340 23 A(I,J)=Y(I) 00007350 C 00007360 C ROWS STORED IN Y AS ARRAY IS REARRANGED IN PLACE. 00007370 C DO 25 I=1,N 00007380 DO 24 J=1,N 00007390 IR=IRR(J) 00007400 JC=JCC(J) 00007410 24 Y(IR)=A(I,JC) 00007420 DO 25 J=1,N 00007430 25 A(I,J)=Y(J) 00007440 GO TO 11 00007450 26 WRITE(6,27) 00007460 27 FORMAT(5X,'MATRIX SHOULD NOT BE LARGER THAN 100*100') 00007470 11 RETURN 00007480 END 00007490 C***********************************************************************NASR C* *NASR C* SUBROUTINE 'MEIGEN' : *NASR C* ------------------- *NASR C* (REFERENCE : APPLIED NUMERICAL METHODS - *NASR C* BY CARNAHAN , PAGES 250 - 260 ) *NASR C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE EIGEN-VECTORS AND *NASR C* THE EIGEN-VALUES OF A (N,N) REAL SYMMETRIC MATRIX (A) ; USING *NASR C* THE 'JACOBI-METHOD' ,WHOSE ALGORITM ITERATIVELY REDUCES THE ORIGIN-*NASR C* AL INPUTTED A-MATRIX TO A NEARLY DIAGONAL MATRIX (A) ,REPLACING THE*NASR C* ORIGINAL ONE ,BY SUCCESSIVELY ANNIHILATING OFF-DIAGONAL ELEMENTS *NASR C* THROUGH A SERIES OF ROTATIONS. THE FINAL PRODUCT OF THESE ROTAT- *NASR C* IONS WILL BE THE ORTHOGONAL MATRIX (R) ,WHOSE N COLUMNS ARE THE N *NASR C* EIGEN-VECTORS OF THE ORIGINAL INPUTTED MATRIX-A. *NASR C* THE N DIAGONAL ELEMENTS OF THE RESULTING DIAGONAL MATRIX (A) WILL *NASR C* BE THE EIGEN-VALUES (N VALUES) OF THE ORIGINAL MATRIX -A. *NASR C* *NASR C* HOWEVER THIS SUBROUTINE IS DESIGNED TO HANDLE THE GENERAL CASE , *NASR C* NAMELY ; IT ACCEPTS THE GENERAL MATRIX A(NA,CA) ,REQUIRED TO DIAG-*NASR C* ONALIZE ITS UPPER LEFT SUBMATRIX A1(N,N) ,TO OBTAIN THE N EIGEN- *NASR C* VALUES ; THEN THE RESULTING CORRESPONDING ORTHOGONAL SUBMATRIX *NASR C* R1(N,N) ,CAN BE STORED IN THE UPPER LEFT PART OF THE INPUTTED *NASR C* DESIRED GENERAL MATRIX R(NR,CR) . *NASR C* *NASR C* USAGE : CALL MEIGEN(A,R,N,NA,NR) *NASR C* *NASR C* THE REQUIRED INPUT DATA ARE : (THROUGH THE CALLING ROUTINE) *NASR C* --------------------------- *NASR C* *NASR C* 1-A : THE ORIGINAL GENERAL MATRIX ,WHICH CONTAINS THE SUBMATRIX *NASR C* REQUIRED TO BE DIAGONALIZED. *NASR C* 2-R : THE GENERAL MATRIX DESIRED TO STORE THE RESULTING ORTHOGONAL*NASR C* SUBMATRIX IN ITS UPPER LEFT PART. *NASR C* 3-N : THE ORDER OF THE SUBMATRIX REQUIRED TO BE DIAGONALIZED. *NASR C* 4-NA : THE ROW DIMENSION OF THE ORIGINAL GENERAL MATRIX (A). *NASR C* 5-NR : THE ROW DIMENSION OF THE GENERAL DESIRED MATRIX (R). *NASR C* *NASR C* NOTE : WHEN DESIRED TO DIAGONALIZE THE WHOLE MATRIX (A) ; *NASR C* N = NA = NR *NASR C* *NASR C* THE OUTPUT WILL BE THE A-MATRIX ,AFTER REPLACING ITS UPPER LEFT *NASR C* PART (N,N) BY THE CORRESPONDING REQUIRED DIAGONAL SUBMATRIX. AND *NASR C* ALSO THE R-MATRIX ,AFTER REPLACING ITS UPPER LEFT PART (N,N) ,BY *NASR C* THE RESULTING ORTHOGONAL SUBMATRIX ,WHICH GIVES THE N EIGEN-VECTORS*NASR C* *NASR C* NOTE :- THIS SUBROUTINE IS VALID FOR A MATRIX OF ORDER UP TO 100 *NASR C* *NASR C***********************************************************************NASR C SUBROUTINE MEIGEN(A,R,N,NA,NR) 0120 IMPLICIT REAL*8(A-H,O-Z) 0130 DIMENSION A(NA,N),R(NR,N),Y(100) 0140 C C SPECIFY ITMAX. MAXIMUM NUMBER OF ITERATIONS 0150 C EPS1. ROTATION ANGLE TEST(ROTATE BY PI/4 IF 0160 C EQUAL DIAGONAL ELEMENTS) 0170 C EPS2. IGNORE OFF-DIAGONAL ELEMENTS .LT.EPS2 0180 C EPS3. BREAKOUT TEST (SUM OF SQUARES OF DIAGONAL 0190 C ELEMENTS UNCHANGED BY LAST ITERATION) 0200 C 0210 ITMAX=50 0220 EPS1=0.1D-12 0230 EPS2=0.1D-12 0240 EPS3=0.1D-10 0250 NM1=N-1 0260 C C SET R=IDENTITY AND UPPER TRIANGULATE A 0270 C DO 16 I=1,N 0280 DO 16 J=1,N 0290 R(I,J)=0. 0300 IF(I.EQ.J) R(I,J)=1. 0310 IF(I.GT.J) A(I,J)=0. 0320 16 CONTINUE 0330 C C COMPUTE SIGMA1 AND S 0340 C SIGMA1=0. 0350 OFFDSQ=0. 0360 DO 17 I=1,N 0370 SIGMA1=SIGMA1+A(I,I)*A(I,I) 0380 IP1=I+1 0390 IF(I.GE.N) GOTO 18 0400 DO 17 J=IP1,N 0410 17 OFFDSQ=OFFDSQ+A(I,J)*A(I,J) 0420 18 S=2.0*OFFDSQ+SIGMA1 0430 C 0440 C BEGIN JACOBI ITERATION 0450 C DO 29 ITER=1,ITMAX 0460 DO 27 I=1,NM1 0470 IP1=I+1 0480 DO 27 J=IP1,N 0490 Q=DABS(A(I,I)-A(J,J)) 0500 C 0510 C COMPUTE SINE AND COSINE OF ROTATION ANGLE 0520 C IF(Q.LE.EPS1) GOTO 19 0530 IF(DABS(A(I,J)).LE.EPS2) GOTO 27 0540 P=2.0*A(I,J)*Q/(A(I,I)-A(J,J)) 0550 SPQ=DSQRT(P*P+Q*Q) 0560 CSA=DSQRT((1.0+Q/SPQ)/2.0) 0570 SNA=P/(2.0*CSA*SPQ) 0580 GOTO 20 0590 19 CSA=1.0/DSQRT(2.0D0) 0600 SNA=CSA 0610 20 CONTINUE 0620 C 0630 C UPDATE COLUMNS I AND J OF R - EQUIVALENT TO 0640 C MULTIPLICATION BY THE ANNIHILATION MATRIX 0650 C DO 21 K=1,N 0660 HOLDKI=R(K,I) 0670 R(K,I)=HOLDKI*CSA+R(K,J)*SNA 0680 21 R(K,J)=HOLDKI*SNA-R(K,J)*CSA 0690 C 0700 C COMPUTE NEW ELEMENTS OF A IN ROWS I AND J 0710 C DO 24 K=1,N 0720 IF(K.GT.J) GOTO 23 0730 Y(K)=A(I,K) 0740 A(I,K)=CSA*Y(K)+SNA*A(K,J) 0750 IF(K.NE.J) GOTO 22 0760 A(J,K)=SNA*Y(K)-CSA*A(J,K) 0770 22 GOTO 24 0780 23 HOLDIK=A(I,K) 0790 A(I,K)=CSA*HOLDIK+SNA*A(J,K) 0800 A(J,K)=SNA*HOLDIK-CSA*A(J,K) 0810 24 CONTINUE 0820 C 0830 C COMPUTE NEW ELEMENTS OF A IN COLUMNS I AND J 0840 C Y(J)=SNA*Y(I)-CSA*Y(J) 0850 C 0860 C WHEN K IS LARGER THAN I 0870 C DO 26 K=1,J 0880 IF(K.LE.I) GOTO 25 0890 A(K,J)=SNA*Y(K)-CSA*A(K,J) 0900 GOTO 26 0910 25 HOLDKI=A(K,I) 0920 A(K,I)=CSA*HOLDKI+SNA*A(K,J) 0930 A(K,J)=SNA*HOLDKI-CSA*A(K,J) 0940 26 CONTINUE 0950 27 A(I,J)=0. 0960 C 0970 C FIND SIGMA2 FOR TRANSFORMED A AND TEST FOR CONVERGENCE 0980 C SIGMA2=0. 0990 DO 28 I=1,N 1000 28 SIGMA2=SIGMA2+A(I,I)*A(I,I) 1010 IF(1.0-SIGMA1/SIGMA2.LT.EPS3) RETURN 1020 29 SIGMA1=SIGMA2 1030 WRITE(6,203) ITER,S,SIGMA1,SIGMA2 1040 203 FORMAT('0','EIGEN NOT CONVERGED'/1X,'ITER =',I5/1X,'S =',F10.5/1X,NASR 1'SIGMA1=',F10.5,/1X,'SIGMA2=',F10.5) NASR RETURN 1070 END