C DATA SET ONESTN AT LEVEL TMP AS OF 05/05/81 C ******** 00001 C *ONESTN* 00002 C ******** 00003 C VERSION JULY 1975 00004 C COMPUTE SINGLE STATION MULTIPASS PHASED SOLUTION FOR AVERAGE 00005 C TERRESTRIAL STATION COORDINATES AND THE FOLLOWING QUASI-OBSERVABLES00006 C RECEIVER FREQUENCY OFFSET 00007 C RECEIVER FREQUENCY DRIFT (NOT YET IMPLEMENTED) 00008 C RECEIVER DELAY ERROR (NOT YET IMPLEMENTED) 00009 C THREE ORBIT TRANSLATIONS 00010 C THREE ORBIT VELOCITIES (NOT YET IMPLEMENTED) 00011 C TROPOSPHERIC REFRACTION SCALING (NOT YET IMPLEMENTED) 00012 C**MACHINE DEPENDENT STATEMENT 00013 IMPLICIT REAL*8(A-H,O-Z) 00014 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00015 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS 00016 COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), 00017 $ DUMMY(260,9) 00018 COMMON /IBR/ STN,IBREAK 00019 COMMON /INIT/ APA(4,4),XGI(3),XGU(3),VPV,IDF 00020 COMMON /CDPX/ DPX(3,3) 00021 COMMON /CPHASE/ PX(3,3),PY(200),PZ(4),AX(200,4),AZ(200,4), 00022 $ Y(200),Z(4),APA1(4,4),DELX(3),VPV1,NX,NY,NZ,NXW,IRC2,IRPX 00023 DIMENSION X0(3),R0(3),PLH(3),SPLH(3,3) 00024 DIMENSION XG(3),SXYZ(3,3) 00025 EQUIVALENCE (FHEAD(2),XG(1)),(FHEAD(5),SXYZ(1,1)) 00026 $ ,(FHEAD(36),A),(FHEAD(37),B),(FHEAD(38),X0(1)),(FHEAD(41),R0(1)),00027 $ (FHEAD(44),S0) 00028 C INITIALIZE RUN 00029 CALL CARDF 00030 NZ = 0 00031 NZ = 3 00032 NZ = 4 00033 C INITIALIZE PHASED SOLUTION 00034 10 CONTINUE 00035 NPASS = 0 00036 VPV = 0. 00037 IDF = 0 00038 DO 15 I = 1,4 00039 DO 15 J = 1,4 00040 15 APA(I,J) = 0D0 00041 DO 20 I = 1,3 00042 XGI(I) = XG(I) 00043 DO 20 J = 1,3 00044 20 APA(I,J) = SXYZ(I,J) 00045 CALL CHOLD(APA,4,3,DETAPA,JCHK) 00046 IF(JCHK .EQ. 0) GO TO 25 00047 WRITE(6,1010) 00048 1010 FORMAT(5X,30HAPA MATRIX IS SINGULAR ) 00049 STOP 00050 25 CONTINUE 00051 WRITE(IPR,1008) XGI,STN,SXYZ 00052 WRITE(IPR,1007) 00053 WRITE(IPR2,1006) STN 00054 C INITIALIZE PASS 00055 30 CONTINUE 00056 NULL = 0 00057 CALL PINIT 00058 CALL RPASS 00059 IF(IEDS .NE. 0) GO TO 40 00060 CALL FORM1 00061 IF(IBREAK .EQ. 0 .OR. NPASS .LT. IBREAK) GO TO 30 00062 C OUTPUT COMPLETED PHASED SOLUTION 00063 40 CONTINUE 00064 IF(NPASS .EQ. 0) RETURN 00065 WRITE(IPR,1003) XGU,STN,DPX 00066 WRITE(IPN,1004) XGU,STN,DPX 00067 WRITE(IPR2,1005) 00068 CALL RECGEO(A,B,X0,R0,S0,XGU,DPX,PLH,SPLH) 00069 WRITE(IPR,1009) PLH,SPLH 00070 IF(IEDS .EQ. 0) GO TO 10 00071 RETURN 00072 C FORMAT STATEMENTS 00073 1003 FORMAT(5X,20HFINAL COORDINATES ,3F15.3,10X,A4/3(/25X,3F15.4)) 00074 1004 FORMAT(10X,3F15.3,10X,A4,3(/3F15.4)) 00075 1005 FORMAT(/) 00076 1006 FORMAT(1H1,10X,22HPHASED SOLUTION FOR ,A4//5X,20(5H----| )) 00077 1007 FORMAT(///50H PASS KEPT NDOP VF NDOPT VFT SAT EL QUAD LOCKON ,00078 * 31H DX DY DZ SIGX SIGY SIGZ,38X,12H95% CONF REG /) 00079 1008 FORMAT(5X,20HINITIAL COORDINATES ,3F15.3,10X,A4/3(/25X,3F15.4)) 00080 1009 FORMAT(5X,20HFINAL COORDINATES ,2F15.8,F15.3/3(/25X,3F15.4)) 00081 END 00082 SUBROUTINE CARDF 00083 C INITIALIZE DATA SET REFERENCE NUMBERS, CONSTANTS AND POINTERS, 00084 C READ INPUT CARDS ONE TO THREE, AND READ ARRAY FHEAD 00085 C**MACHINE DEPENDENT STATEMENT 00086 IMPLICIT REAL*8(A-H,O-Z) 00087 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00088 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS 00089 COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), 00090 $ DUMMY(260,9) 00091 COMMON /IBR/ STN,IBREAK 00092 DIMENSION TITLE(20) 00093 DIMENSION IUP(12),UP(12) 00094 DIMENSION SD(3),SXYZ(3,3),SEPHSD(3),VEPHSD(3) 00095 EQUIVALENCE (FHEAD(5),SXYZ(1,1)),(FHEAD(61),DOPSIG),(FHEAD(65), 00096 $ FGOSIG),(FHEAD(66),DFGSIG),(FHEAD(59),DELSIG),(FHEAD(46), 00097 $ SEPHSD(1)),(FHEAD(49),VEPHSD(1)),(FHEAD(85),TROPSD) 00098 SQRT(PI) = DSQRT(PI) 00099 C DATA SET REFERENCE NUMBERS 00100 C MET=WEATHER INPUT, IN=MJV INPUT, IOUT=OUTPUT RESULTS, IOUT2=SCRATCH,00101 C ICR=CARD READER, IPR=PRINTER1, IPN=PUNCH, IPR2=PRINTER2, IPR3=PRINT300102 MET = 1 00103 IN = 2 00104 IOUT = 3 00106 IOUT2= 4 00107 ICR = 5 00108 IPR = 6 00109 IPN=3 **2 IPR2 = 8 00111 IPR3 = 9 00112 C MACHINE DEPENDENT PARAMETERS. VEL=VACUUM VELOCITY OF PROP (M/SEC) 00113 C WE=ROTATION RATE OF EARTH (RAD/MIN). FO=GROUND REFERENCE FREQ (HZ) 00114 PI = 3.141592653589793D0 00115 VEL = 299792.5D3 00116 WE = 4.3752695D-3 00117 FO = 400D6 00118 C DERIVED CONSTANTS 00119 RADG = PI / 180D0 00120 WAVENO = FO / VEL 00121 C BIT PERIODS FOR BR(SATELLITE) AND CBR(RECEIVER) TIMING 00122 BITS = 120D0 / 6103D0 00123 BITR = 7865000D0 / FO 00124 C POINTERS 00125 IKEPT = 0 00126 IEDS = 0 00127 C RESET FILE HEADER ARRAY 00128 DO 10 I = 1,200 00129 10 FHEAD(I) = 0. 00130 CARD ONE CONTAINS LISTING OPTION,FOUR CHARACTER STATION NAME, AND 00131 C NUMBER OF PASSES FOR EACH PHASED SOLUTION 00132 C ILIST = LEVEL OF OUTPUT LISTING 00133 C 0 - ONE LINE SUMMARY OF EACH PASS ONLY 00134 C 1 - SUMMARY PLUS ERROR MESSAGES 00135 C 2 - COMPLETE LISTING OF COMPUTED DATA 00136 C 3 - COMPLETE INPUT AND OUTPUT LISTING 00137 C IBREAK = NUMBER OF PASSES TO BE USED IN EACH PHASED SOLUTION 00138 C = 0 IF ALL PASSES ARE USED IN ONE PHASED SOLUTION 00139 READ(ICR,1004) ILIST,STN,IBREAK 00140 IF(ILIST .GE. 1) WRITE(IPR,1000) 00141 CARD TWO DESCRIBES THIS COMPUTER RUN IN UP TO 80 CHARACTERS (20A4) 00142 READ(ICR,1001) TITLE 00143 WRITE(IPR,2001) TITLE 00144 IF(IBREAK .NE. 0) WRITE(IPR,2002) IBREAK 00145 IF(IBREAK .EQ. 0) WRITE(IPR,2003) 00146 CARDS THREE AND FOUR CONTAIN UP TO 12 UPDATES TO FHEAD PARAMETERS 00147 C IUP(I) = INDEX OF FHEAD CONTAINING VALUE TO BE UPDATED 00148 C UP(I) = NEW VALUE 00149 READ(ICR,1003) (IUP(I),UP(I),I=1,12) 00150 C READ ARRAY FHEAD 00151 READ(IN,1002) FHEAD 00152 IF(ILIST .GE. 3) WRITE(IPR,1002) FHEAD 00153 C UPDATE FHEAD 00154 NUP = 0 00155 DO 20 I = 1,12 00156 IF(IUP(I) .EQ. 0) GO TO 30 00157 NUP = NUP + 1 00158 WRITE(IPR,1005) IUP(I),FHEAD(IUP(I)),UP(I) 00159 FHEAD(IUP(I)) = UP(I) 00160 20 CONTINUE 00161 30 CONTINUE 00162 C LIST STANDARD DEVIATIONS 00163 DO 40 I = 1,3 00164 40 SD(I) = SQRT(SXYZ(I,I)) 00165 WRITE(IPR,1006) SD,DOPSIG,FGOSIG,DFGSIG,DELSIG,SEPHSD,VEPHSD, 00166 $ TROPSD 00167 RETURN 00168 C FORMAT STATEMENTS 00169 1000 FORMAT(1X,30HSUBROUTINE CARDF ENTERED ) 00170 1001 FORMAT(20A4) 00171 1002 FORMAT(1X,F10.0,3F13.3,3F10.0/1X,6F10.0,F13.8/1X,F13.8,F13.3,5F10. **2 $ /1X,4F10.0,3F13.10/1X,6F13.10/1X,F13.3,2F10.1/1X,3F9.3,4F9.4/ 00173 $ 1X,F3.0,6F5.2,F3.0,F5.2,F3.0,F6.0,F3.0,F5.0,F6.0,F6.0,F3.0,F6.2/00174 $ 1X,F8.1,F6.1,F5.0,F8.0,F6.0,F5.1,F5.0, F5.1,F7.1,F5.2/ 00175 $ 1X,F5.1,F4.0,F4.0,F5.3,F4.1,F6.0,F5.0,F4.0,F5.2,F3.0,F5.0,F6.0, 00176 $ F4.0,F5.2/1X,F4.0,12F6.0/1X,13F6.0/1X,F6.0,F4.0,11F6.0/1X,13F6.000177 $ /1X,2F6.0,15F4.0/1X,10F5.0,5F5.1/1X,15F5.1/1X,15F5.0/1X,F5.0) 00178 1003 FORMAT(6(I3,F10.0)) 00179 1004 FORMAT(I1,A4,I3) 00180 1005 FORMAT(5X,13HUPDATE FHEAD(I5,7H ) FROM ,F20.10,4H TO ,F20.10) 00181 1006 FORMAT(1X,30HSTANDARD DEVIATIONS ASSIGNED //5X,11HCOORDINATES,9X,00182 $ 3F15.3//5X,20HDOPPLER OBSERVATIONS ,F15.3//5X,10HQUASI-OBS / 00183 $ 10X,15HFREQ OFFSET ,F15.3/10X,15HFREQ DRIFT ,F15.3/10X, 00184 $ 15HRCVR DELAY ,F15.3/10X,15HORBIT TRANS ,3F15.3/10X, 00185 $ 15HORBIT VELOC ,3F15.3/10X,15HTROPO SCALING ,F15.3/) 00186 2001 FORMAT(1H1//50X,30HPROGRAM **ONESTN** OUTPUT //,20X,20A4) 00187 2002 FORMAT(5X,I3,40H PASSES PER INDEPENDENT PHASED SOLUTION ) 00188 2003 FORMAT(5X,40HALL PASSES USED IN ONE PHASED SOLUTION ) 00189 END 00190 SUBROUTINE PINIT 00191 C INITIALIZE PASS ARRAYS 00192 IMPLICIT REAL*8(A-H,O-Z) 00193 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00194 COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), 00195 $ DUMMY(260,9) 00196 IF(ILIST .GT. 0) WRITE(IPR,1000) 00197 DO 10 I = 1,100 00198 10 PHEAD(I) = 0. 00199 DO 20 I = 1,20 00200 ORB1(I) = 0. 00201 20 ORB3(I) = 0. 00202 DO 30 I = 1,70 00203 30 ORB2(I) = 0. 00204 DO 50 I = 1,260 00205 DO 50 J = 1,9 00206 50 DUMMY(I,J) = 0D0 00207 RETURN 00208 1000 FORMAT(1X,30HSUBROUTINE PINIT ENTERED ) 00209 END 00210 SUBROUTINE RPASS 00211 C READ PASS HEADER AND DOP1 AND DOP2 VECTORS 00212 C**MACHINE DEPENDENT STATEMENT 00213 IMPLICIT REAL*8(A-H,O-Z) 00214 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00215 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS 00216 COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), 00217 $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), 00218 $ DOP7(260),DOP8(260),DOP9(260) 00219 EQUIVALENCE (PHEAD(26),DOPNUM) 00220 READ(IN,1001,END=100,ERR=110) PHEAD 00221 IF(ILIST .GE. 1) WRITE(IPR,1001) PHEAD 00222 READ(IN ,1002,END=100,ERR=110) ORB1 00223 IF(ILIST .GE. 1) WRITE(IPR,1002) ORB1 00224 READ(IN ,1004,END=100,ERR=110) ORB3 00225 IF(ILIST .GE. 1) WRITE(IPR,1004) ORB3 00226 NDOP = DOPNUM 00227 NDOP1 = NDOP + 1 00228 DO 10 J = 1,NDOP1 00229 READ(IN,1008,END=100,ERR=110) I,DOP1(I),DOP2(I),DOP3(I),DOP4(I),00230 $ DOP5(I),DOP6(I),DOP7(I),DOP8(I),DOP9(I) 00231 IF(ILIST .GE. 1) WRITE(IPR,1008) I,DOP1(I),DOP2(I),DOP3(I), 00232 $ DOP4(I),DOP5(I),DOP6(I),DOP7(I),DOP8(I),DOP9(I) 00233 10 CONTINUE 00234 RETURN 00235 C END OF DATA 00236 100 CONTINUE 00237 WRITE(IPR,1006) 00238 IEDS = 1 00239 RETURN 00240 C INPUT FORMAT ERROR 00241 110 WRITE(IPR,1009) 00242 STOP 00243 C FORMAT STATEMENTS 00244 1001 FORMAT(1X,F5.0,F4.0,F5.0,F6.0,F5.0,F6.0,F12.6,F5.1,2F5.0,F3.0, 00245 $ F6.1,F7.1,F5.1/1X,2F5.1,6F4.0,6F5.0,F3.0/1X,3F13.3/1X,9F8.0/ 00246 $ 1X,6F13.10/1X,3F13.3,F4.0,F3.0,F12.6,F7.3,F4.0/ 00247 $ 1X,F3.0,F5.0,F5.1,6F7.1/1X,36F2.0) 00248 1002 FORMAT(1X,5F15.7) 00249 1003 FORMAT(1X,11F7.0) 00250 1004 FORMAT(1X,8F9.4) 00251 1006 FORMAT(1X,30HEND OF DATA IN RPASS ) 00252 1007 FORMAT(1X,3(F10.6,F11.3,F5.0)) 00253 1008 FORMAT(1X,I3,2F11.2,F10.6,F12.3,F7.3,F8.3,F5.1,F5.0,F7.3) 00254 1009 FORMAT(1X,30HINPUT FORMAT ERROR IN RPASS ) 00255 END 00256 C 00257 SUBROUTINE FORM1 00258 C COMPUTE DESIGN MATRICES 00259 C CALL PHASE1 00260 C UPDATE PHASED SOLUTION FOR ACCEPTED PASSES 00261 C PRINT PASS SUMMARY 00262 C PLOT PASS 00263 C**MACHINE DEPENDENT STATEMENT 00264 IMPLICIT REAL*8(A-H,O-Z) 00265 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00266 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS 00267 COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), 00268 $ DOP1(260),DOP2(260),DOP3(260),DOP4(260),DOP5(260),DOP6(260), 00269 $ DOP7(260),DOP8(260),DOP9(260) 00270 COMMON /CDPX/ DPX(3,3) 00271 COMMON /INIT/ APA(4,4),XGI(3),XGU(3),VPV,IDF 00272 COMMON /CPHASE/ PX(3,3),PY(200),PZ(4),AX(200,4),AZ(200,4), 00273 $ Y(200),Z(4),APA1(4,4),DELX(3),VPV1,NX,NY,NZ,NXW,IRC2,IRPX 00274 DIMENSION DSDX1(3),DSDX2(3),DSDC1(3),DSDC2(3) 00275 DIMENSION VECTOR(3,3),EIG(3) 00276 DIMENSION SD(3) 00277 DIMENSION IPLOT(100) 00278 DIMENSION IX(3),IZ(3) 00279 DATA IX(1),IX(2),IX(3),IZ(1),IZ(2),IZ(3)/1HX,1HY,1HZ,70,40,10/ 00280 DATA IBLANK,IS,IDOT,IMIN,ID/1H ,1HS,1H.,1H-,1HD/ 00281 DIMENSION SEPHSD(3) 00282 DIMENSION KWAD(4) 00283 DATA KWAD(1),KWAD(2),KWAD(3),KWAD(4)/4HW NS,4HW SN,4HE NS,4HE SN/ 00284 EQUIVALENCE (FHEAD(46),SEPHSD(1)),(FHEAD(58),DELAY),(FHEAD(60), 00285 $ TCODE),(FHEAD(61),DOPSIG),(FHEAD(65),FGOSIG),(FHEAD(70),YMISCL),00286 $ (FHEAD(73),DOPMIN),(PHEAD(26),DOPNUM),(ORB1(13),FG),(ORB1(14), 00287 $ OFFREQ),(PHEAD(2),SAT),(PHEAD(3),DLOCK),(PHEAD(4),TLOCK), 00288 $ (PHEAD(8),ELEV),(PHEAD(11),QUAD) 00289 $ ,(PHEAD(61),GUIERF) 00290 C**MACHINE DEPENDENT STATEMENT 00291 ABS(PI) = DABS(PI) 00292 COS(PI) = DCOS(PI) 00293 SIN(PI) = DSIN(PI) 00294 SQRT(PI) = DSQRT(PI) 00295 AMOD(PI,VEL) = DMOD(PI,VEL) 00296 ATAN2(PI,VEL) = DATAN2(PI,VEL) 00297 C INITIALIZATION 00298 IF(ILIST .GE. 1) WRITE(IPR,1000) 00299 IRC2 = 4 00300 IRPX = 3 00301 NX = 3 00302 NXW = NX + 1 00303 X60 = 60D0 00304 ISAT = SAT 00305 IELEV = ELEV 00306 IQUAD = QUAD 00307 IQUAD = KWAD(IQUAD) 00308 LODAY = DLOCK 00309 LOHR = TLOCK / 60D0 00310 LOMIN = AMOD(TLOCK,X60) 00311 MINDOP = DOPMIN 00312 NDOP = DOPNUM 00313 FCA = (OFFREQ + FG + GUIERF) * FO / 1D10 00314 WAVE = (1D0 + (FG + GUIERF) / 1D10) * FO / VEL 00315 IF(ILIST .GE. 3) WRITE(IPR,1002) FCA,WAVE 00316 IF(ILIST .GE. 2) WRITE(IPR,1001) 00317 C COMPUTE DESIGN MATRICES 00318 NOBS = 0 00319 T = 0. 00320 IF(TCODE .EQ. 2.) T = DELAY/1D6 00321 CALL DERIV(T,DSDX1,DSDC1,S1) 00322 DO 80 I = 1,NDOP 00323 T = T + DOP3(I) 00324 CALL DERIV(T,DSDX2,DSDC2,S2) 00325 IF(DOP8(I) .GT. 0.) GO TO 10 00326 GO TO 60 00327 10 XMISCL = DOP4(I) - DOP3(I) * FCA - WAVE * (S2 - S1) 00328 IF(ABS(XMISCL) .LT. YMISCL) GO TO 20 00329 DOP8(I) = -10. 00330 GO TO 60 00331 20 NOBS = NOBS + 1 00332 C COORDINATE DESIGN MATRIX 00333 DO 30 K = 1,NX 00334 30 AX(NOBS,K) = -WAVE * (DSDX2(K) - DSDX1(K)) 00335 AX(NOBS,NXW) = - XMISCL 00336 C QUASI-OBSERVABLE DESIGN MATRIX 00337 DO 40 K = 1,3 00338 40 AZ(NOBS,K) = WAVE * (DSDC2(K) - DSDC1(K)) 00339 AZ(NOBS,4) = (DOP3(I) + (S2 - S1) / VEL) * FO / 1D10 00340 C OBSERVABLE WEIGHT MATRIX 00341 PY(NOBS) = 1D0 / DOPSIG**2 00342 IF(ILIST .GE. 3) WRITE(IPR,1013) NOBS,(AX(NOBS,K),K=1,NXW), 00343 $ (AZ(NOBS,K),K=1,NZ),PY(NOBS) 00344 60 S1 = S2 00345 DO 70 K = 1,3 00346 DSDX1(K) = DSDX2(K) 00347 70 DSDC1(K) = DSDC2(K) 00348 80 CONTINUE 00349 C QUASI-OBSERVABLE WEIGHT MATRIX 00350 DO 90 K = 1,3 00351 90 PZ(K) = 1D0 / SEPHSD(K)**2 00352 PZ(4) = 1D0 / FGOSIG**2 00353 IF(ILIST .GE. 3) WRITE(IPR,1015) PZ 00354 C REJECT PASS IF LT MINDOP DOPPLERS 00355 IF(NOBS .LT. MINDOP) GO TO 220 00356 C PERFORM PHASED ADJUSTMENT 00357 NY = NOBS 00358 DO 100 I = 1,NXW 00359 DO 100 J = 1,NXW 00360 100 APA1(I,J) = APA(I,J) 00361 IF(ILIST .GE. 2) WRITE(IPR,1014) APA 00362 CALL PHASE1(NULL) 00363 IF(NULL .NE. 0) RETURN 00364 C COMPUTE DEGREES OF FREEDOM FOR THIS PASS 00365 JDF = NY 00366 IF(NPASS .EQ. 0) JDF = JDF - 3 00367 DVPV = VPV1 - VPV 00368 IF(ILIST .GE. 1) WRITE(IPR,2004) VPV,DVPV,VPV1 00369 DSIGMA = DVPV/JDF 00370 C REJECT PASS IF IT FAILS CHI-SQUARE TEST ON RESIDUAL NORM 00371 IF(DSIGMA .GT. CHISQ(0.975D0,JDF)/JDF) GO TO 200 00372 C ACCEPT PASS AND UPDATE NUMBER OF PASSES, RESIDUAL NORM, 00373 C DEGREES OF FREEDOM, SOLUTION VECTOR AND ITS WEIGHT MATRIX 00374 NPASS = NPASS + 1 00375 IDF = IDF + JDF 00376 VPV = VPV1 00377 SIGMA = VPV / IDF 00378 DO 136 I = 1,NXW 00379 DO 136 J = 1,NXW 00380 136 APA(I,J) = APA1(I,J) 00381 DO 137 I = 1,NX 00382 DO 137 J = 1,NX 00383 PX(I,J) = SIGMA * PX(I,J) 00384 DPX(I,J) = PX(I,J) 00385 137 CONTINUE 00386 DO 138 I = 1,NX 00387 XGU(I) = XGI(I) + DELX(I) 00388 138 SD(I) = SQRT(PX(I,I)) 00389 C PRINT ONE LINE SUMMARY OF PASS 00390 IF(NZ .EQ. 0) WRITE(IPR,1006) NPASS,NY,DSIGMA,IDF,SIGMA,ISAT, 00391 $ IELEV,IQUAD,LODAY,LOHR,LOMIN,DELX,SD 00392 IF(NZ .GT. 0) WRITE(IPR,1006) NPASS,NY,DSIGMA,IDF,SIGMA,ISAT, 00393 $ IELEV,IQUAD,LODAY,LOHR,LOMIN,DELX,SD,(Z(I),I=1,NZ) 00394 IF(ILIST .GE. 1) WRITE(IPR,1007) XGU,IDF,VPV,SIGMA,PX 00395 C COMPUTE SEMI-AXES OF 95 PERCENT CONFIDENCE REGION 00396 DO 145 I = 1,NX 00397 145 EIG(I) = 0. 00398 CALL EIGEN(PX,VECTOR,NX,KCHK) 00399 IF(KCHK .NE. 0) GO TO 230 00400 DO 146 I = 1,NX 00401 146 EIG(I) = SQRT(PX(I,I)) * 2.79D0 00402 WRITE(IPR,1008) EIG 00403 C PLOT RESULTS ON DATASET IPR2 00404 DO 150 J = 1,100 00405 150 IPLOT(J) = IBLANK 00406 IAX = IDOT 00407 IF(MOD(NPASS,10) .EQ. 0) IAX = IMIN 00408 DO 160 J = 1,3 00409 J1 = DELX (J) + IZ(J) 00410 J2 = SD(J) + IZ(J) 00411 J3 = EIG(J) + IZ(J) 00412 IF(J1 .GT. 100) J1 = 100 00413 IF(J2 .GT. 100) J2 = 100 00414 IF(J3 .GT. 100) J3 = 100 00415 IF(J1 .LE. 0) J1 = 1 00416 IF(J2 .LE. 0) J2 = 1 00417 IF(J3 .LE. 0) J3 = 1 00418 IPLOT(J2) = IS 00419 IPLOT(J1) = IX(J) 00420 IPLOT(J3) = IDOT 00421 IPLOT(IZ(J)) = IAX 00422 160 CONTINUE 00423 WRITE(IPR2,1004) IPLOT 00424 RETURN 00425 C ERROR RETURNS 00426 200 WRITE(IPR,1009) DSIGMA,ISAT,IELEV,IQUAD,LODAY,LOHR,LOMIN 00427 NULL = 14 00428 RETURN 00429 210 WRITE(IPR,1010) 00430 NULL = 16 00431 RETURN 00432 220 WRITE(IPR,1011) MINDOP,ISAT,IELEV,IQUAD,LODAY,LOHR,LOMIN 00433 NULL = 15 00434 RETURN 00435 230 WRITE(IPR,1012) 00436 NULL = 17 00437 RETURN 00438 C FORMAT STATEMENTS 00439 1000 FORMAT(1X,30HSUBROUTINE FORM1 ENTERED ) 00440 1001 FORMAT(5X,7HMSG ROW,5X,1HS,5X,4HNOBS,5X,3HDOP,1HF,5X,2HAX,5X,2HAZ)00441 1002 FORMAT(5X,8HFCA WAVE,2F20.10) 00442 1003 FORMAT(1H+,22X,I3,F8.0,7F11.5,F11.3) 00443 1004 FORMAT(5X,100A1,5X) 00444 1006 FORMAT(5X,I4,I5,F5.1,I6,F5.1,2I3,1X,A4,I4,2I3,12F7.3) 00445 1007 FORMAT(5X,15HSOLUTION VECTOR,3F15.5/5X,12HDF VPV SIGMA,I10,2F15.5/00446 * 5X,10HCOV MAT ,3E15.2/15X,3E15.2/15X,3E15.2) 00447 1008 FORMAT(1H+,118X,3F4.0) 00448 1009 FORMAT(5X,5X,15HREJECT-CHISQ ON,F5.0,2I3,1X,A4,I4,2I3) 00449 1010 FORMAT(10X,25HSINGULAR MATRIX IN CHOLD ) 00450 1011 FORMAT(5X,5X,9HREJECT-LT,I6,5H DOPS,2I3,1X,A4,I4,2I3) 00451 1012 FORMAT(5X,15HPASS ABORTED ) 00452 1013 FORMAT(5X,I5,9E12.3) 00453 1014 FORMAT(5X,3HAPA/4(5X,4E10.3/)) 00454 1015 FORMAT(5X,10E12.3) 00455 2004 FORMAT(5X,3HVPV/5X,F10.3,3H + ,F10.3,3H = ,F10.3) 00456 END 00457 SUBROUTINE DERIV(SEC,DSDX,DSDC,S) 00458 C**MACHINE DEPENDENT STATEMENT 00459 IMPLICIT REAL*8(A-H,O-Z) 00460 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00461 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS 00462 COMMON /ARRAYS/ FHEAD(200),PHEAD(100),ORB1(20),ORB2(70),ORB3(20), 00463 $ DUMMY(260,9) 00464 DIMENSION R(5),ROT(3,3),NAX(5),NAX1(3),X(3),XS(3),CDE(4),CDA(4), 00465 $ CET(4) 00466 DIMENSION XSO(3),XSG(3),DXSODC(3,3),DXSGDC(3,3),DELX(3),DSDX(3), 00467 $ DSDC(3) 00468 DIMENSION XG(3),TOPROT(3,3),GPROT(3,3) 00469 EQUIVALENCE (ORB1(1),TP),(ORB1(2),XN),(ORB1(3),WNOT),(ORB1(4), 00470 $ WDOT),(ORB1(5),EXC),(ORB1(6),AO),(ORB1(7),RA),(ORB1(8),RADOT), 00471 $ (ORB1(9),XI),(ORB1(10),GAST),(ORB1(11),DAYTP),(FHEAD(2),XG(1)), 00472 $ (FHEAD(26),TOPROT(1,1)),(PHEAD(3),DLOCK),(PHEAD(4),TLOCK), 00473 $ (PHEAD(15),XPOLE),(PHEAD(16),YPOLE),(PHEAD(17),STAT),(PHEAD(42),00474 $ GPROT(1,1)),(ORB3(1),CDE(1)),(ORB3(5),CDA(1)),(ORB3(9),CET(1)) 00475 DATA NAX(1),NAX(2),NAX(3),NAX(4),NAX(5)/3,1,3,2,1/ 00476 C**MACHINE DEPENDENT STATEMENT 00477 SIN(PI) = DSIN(PI) 00478 COS(PI) = DCOS(PI) 00479 SQRT(PI) = DSQRT(PI) 00480 IF(ILIST .GE. 3) WRITE(IPR,1000) 00481 C COMPUTE EPHEMERAL PARAMETERS AT T 00482 T = SEC / 60D0 00483 CT = COS(2. * T * XN * RADG) 00484 ST = SIN(2. * T * XN * RADG) 00485 DE = (CDE(1) + CDE(2) * CT + CDE(3) * ST + CDE(4) * T) * RADG/1D4 00486 DA = (CDA(1) + CDA(2) * CT + CDA(3) * ST + CDA(4) * T) * 1D1 00487 ET = (CET(1) + CET(2) * CT + CET(3) * ST + CET(4) * T) * 1D1 00488 C COMPUTE SATELLITE COORDINATES IN ORBITAL COORDINATE SYSTEM 00489 UT = T + TLOCK 00490 DT = UT - TP + (DLOCK - DAYTP) * 1440. 00491 XE = XN * RADG * DT + EXC * SIN(XN * RADG * DT) + DE 00492 XA = AO + DA 00493 C COMPUTE ROTATION MATRIX FROM ORBITAL TO AVERAGE TERRESTRIAL SYSTEM 00494 R(1) = - (WNOT - WDOT * DT) * RADG 00495 R(2) = - XI * RADG 00496 R(3) = -(RA + RADOT * DT) * RADG + GAST * RADG + WE * DT 00497 R(4) = - XPOLE / 6356D3 00498 R(5) = -YPOLE / 6356D3 00499 CALL ROTREF(5,NAX,R,ROT) 00500 COSXE = COS(XE) 00501 SINXE = SIN(XE) 00502 XSO(1) = XA * (COSXE - EXC) 00503 XSO(2) = XA * SINXE 00504 XSO(3) = ET 00505 C COMPUTE SATELLITE COORDINATES IN GEOCENTRIC SYSTEM 00506 DO 30 K = 1,3 00507 XSG(K) = 0. 00508 DO 30 L = 1,3 00509 30 XSG(K) = XSG(K) + ROT(K,L) * XSO(L) 00510 C COMPUTE DERIVATIVES OF SATELLITE ORBIT COORDINATES WITH 00511 C RESPECT TO ORBIT BIAS COEFFICIENTS 00512 DO 40 K = 1,3 00513 DO 40 L = 1,3 00514 40 DXSODC(K,L) = 0. 00515 DXSODC(1,1) = -XA * SINXE * RADG / 1D4 00516 DXSODC(1,2) = XA * COSXE * RADG / 1D4 00517 DXSODC(2,1) = (COSXE - EXC) * 1D1 00518 DXSODC(2,2) = SINXE * 1D1 00519 DXSODC(3,3) = 10D0 00520 C COMPUTE DERIVATIVES OF SATELLITE GEOCENTRIC COORDINATES 00521 C WITH RESPECT TO ORBIT BIAS COEFFICIENTS 00522 DO 50 K = 1,3 00523 DO 50 L = 1,3 00524 DXSGDC(K,L) = 0. 00525 DO 50 M = 1,3 00526 50 DXSGDC(K,L) = DXSGDC(K,L) + ROT(K,M) * DXSODC(M,L) 00527 C COMPUTE SLANT RANGE 00528 S = 0. 00529 DO 60 K = 1,3 00530 DELX(K) = XSG(K) - XG(K) 00531 60 S = S + DELX(K) * DELX(K) 00532 S = SQRT(S) 00533 C COMPUTE DERIVATIVES OF SLANT RANGE WITH RESPECT TO GEOCENTRIC COORDS 00534 DO 70 K = 1,3 00535 70 DSDX(K) = DELX(K) / S 00536 C COMPUTE DERIVATIVES OF SLANT RANGE WITH RESPECT TO ORBIT BIAS COEFFS 00537 DO 80 K = 1,3 00538 DSDC(K) = 0. 00539 DO 80 L = 1,3 00540 80 DSDC(K) = DSDC(K) + DXSGDC(K,L) * DSDX(L) 00541 IF(ILIST .GE. 3) WRITE(IPR,1002) SEC,DSDX,DSDC,S 00542 RETURN 00543 1000 FORMAT(1X,30HSUBROUTINE DERIV ENTERED ) 00544 1002 FORMAT(5X,F10.6,6E12.3,F13.3) 00545 END 00546 SUBROUTINE PHASE1(ICHK) 00547 C PHASED SOLUTION FOR ONE STATION ( 3 COORDINATES) 00548 C ALGORITHM GIVEN IN UNBSE REPORT 29 PAGES 98/99 00549 C INPUT 00550 C APA1= 4X4 ARRAY CONTAINING N(I-1) IN UPPER LEFT 3X3 00551 C U(I-1) IN RIGHT AND BOTTOM 1X3 00552 C RHO(W)(I-1) IN BOTTOM RIGHT 1X1 00553 C AX = NYX4 ARRAY CONTAINING A(I) IN LEFT 3 COLUMNS 00554 C W(I) IN RIGHT COLUMN 00555 C AZ = NYXNZ ARRAY CONTAINING D(I) 00556 C PY = NY VECTOR CONTAINING DIAGONAL ELEMENTS OF PL(I) MATRIX 00557 C OUTPUT 00558 C APA1 = UPDATED APA (REPLACES APA ONLY IF PASS ACCEPTED) 00559 C PX = 3X3 ARRAY CONTAINING N(I) INVERSE 00560 C DELX = 3 VECTOR CONTAINING - PX * U(I) 00561 C PZ = NZ VECTOR CONTAINING DIAGONAL ELEMENTS OF PZ(I) MATRIX 00562 C VPV1 = SCALAR CONTAINING RHO(W)(I) + U(I) TRANS * DELX 00563 C Y = NY VECTOR CONTAINING L(I) 00564 C Z = NZ VECTOR CONTAINING Z(I) 00565 C**MACHINE DEPENDENT STATEMENT 00566 IMPLICIT REAL*8(A-H,O-Z) 00567 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00568 COMMON /CPHASE/ PX(3,3),PY(200),PZ(4),AX(200,4),AZ(200,4), 00569 $ Y(200),Z(4),APA1(4,4),DELX(3),VPV1,NX,NY,NZ,NXW,IRC2,IRPX 00570 DIMENSION C1(4,200),C2(4,4),C3(4,4),C4(4,4) 00571 IF(ILIST .GE. 1) WRITE(IPR,1000) 00572 ICHK = 0 00573 C COMPUTE APA1 = APA + AX(T) * PY * AX 00574 DO 122 I = 1,NXW 00575 DO 122 J = 1,NXW 00576 DO 122 K = 1,NY 00577 122 APA1(I,J) = APA1(I,J) + AX(K,I) * PY(K) * AX(K,J) 00578 IF(ILIST .GE. 2) WRITE(IPR,2001) APA1 00579 IF(NZ .LE. 0) GO TO 128 00580 C COMPUTE APA1 = APA1 - AX(T)*PY*AZ*(PZ+AZ(T)*PY*AZ)(-1)*AZ(T)*PY*AX 00581 C SET C1 = AZ(T) * PY 00582 DO 123 I = 1,NZ 00583 DO 123 J = 1,NY 00584 123 C1(I,J) = AZ(J,I) * PY(J) 00585 C SET C2 = PZ + AZ(T)*PY*AZ 00586 DO 124 I = 1,NZ 00587 DO 124 J = 1,NZ 00588 C2(I,J) = 0. 00589 IF(I .EQ. J) C2(I,J) = PZ(I) 00590 DO 124 K = 1,NY 00591 124 C2(I,J) = C2(I,J) + C1(I,K) * AZ(K,J) 00592 IF(ILIST .GE. 3) WRITE(IPR,2007) C2 00593 C INVERT C2 00594 CALL CHOLD(C2,IRC2,NZ,DET,JCHK) 00595 IF(JCHK .NE. 0) GO TO 210 00596 IF(ILIST .GE. 3) WRITE(IPR,2007) C2 00597 C SET C3 = AZ(T)*PY*AX 00598 DO 125 I = 1,NZ 00599 DO 125 J = 1,NXW 00600 C3(I,J) = 0. 00601 DO 125 K = 1,NY 00602 125 C3(I,J) = C3(I,J) + C1(I,K) * AX(K,J) 00603 IF(ILIST .GE. 3) WRITE(IPR,2008) C3 00604 C SET C4 = (PZ+AZ(T)*PY*AZ)(-1)*AZ(T)*PY*AX 00605 DO 126 I = 1,NZ 00606 DO 126 J = 1,NXW 00607 C4(I,J) = 0. 00608 DO 126 K = 1,NZ 00609 126 C4(I,J) = C4(I,J) + C2(I,K) * C3(K,J) 00610 IF(ILIST .GE. 3) WRITE(IPR,2009) C4 00611 DO 127 I = 1,NXW 00612 DO 127 J = 1,NXW 00613 DO 127 K = 1,NZ 00614 127 APA1(I,J) = APA1(I,J) - C3(K,I) * C4(K,J) 00615 IF(ILIST .GE. 2) WRITE(IPR,2010) APA1 00616 128 CONTINUE 00617 C COMPUTE PX 00618 DO 129 I = 1,NX 00619 DO 129 J = 1,NX 00620 129 PX(I,J) = APA1(I,J) 00621 IF(ILIST .GE. 2) WRITE(IPR,2002) PX 00622 C INVERT PX 00623 CALL CHOLD(PX,IRPX,NX,DET,JCHK) 00624 IF(JCHK .NE. 0) GO TO 210 00625 IF(ILIST .GE. 2) WRITE(IPR,2002) PX 00626 C COMPUTE DELX 00627 DO 130 I = 1,NX 00628 DELX(I) = 0. 00629 DO 130 J = 1,NX 00630 130 DELX(I) = DELX(I) - PX(I,J) * APA1(J,NXW) 00631 IF(ILIST .GE. 2) WRITE(IPR,2003) DELX 00632 C COMPUTE VPV1 00633 VPV1 = APA1(NXW,NXW) 00634 DO 131 I = 1,NX 00635 131 VPV1 = VPV1 + APA1(NXW,I) * DELX(I) 00636 C COMPUTE Y 00637 DO 132 I = 1,NY 00638 Y(I) = AX(I,NXW) 00639 DO 132 J = 1,NX 00640 132 Y(I) = Y(I) + AX(I,J) * DELX(J) 00641 IF(NZ .EQ. 0) GO TO 135 00642 DO 133 I = 1,NZ 00643 Z(I) = -C4(I,NXW) 00644 DO 133 J = 1,NX 00645 133 Z(I) = Z(I) - C4(I,J) * DELX(J) 00646 IF(ILIST .GE. 1) WRITE(IPR ,2011) (Z(I),I=1,NZ) 00647 DO 134 I = 1,NY 00648 DO 134 J = 1,NZ 00649 134 Y(I) = Y(I) + AZ(I,J) * Z(J) 00650 135 CONTINUE 00651 IF(ILIST .GE. 1) WRITE(IPR ,2005) (Y(I),I=1,NY) 00652 RETURN 00653 C ERROR RETURN 00654 210 WRITE(IPR ,1010) 00655 ICHK = 1 00656 RETURN 00657 C FORMAT STATEMENTS 00658 1000 FORMAT(1X,30HSUBROUTINE PHASE1 ENTERED ) 00659 1010 FORMAT(10X, 30HSINGULAR WT MATRIX IN CHOLD ) 00660 2001 FORMAT(5X,4HAPA1/4(5X,4E10.3/)) 00661 2002 FORMAT(5X,2HPX/3(5X,3E10.3/)) 00662 2003 FORMAT(5X,4HDELX/5X,3E10.3) 00663 2005 FORMAT(5X,15F7.3) 00664 2007 FORMAT(5X,2HC2/4(5X,4E10.3/)) 00665 2008 FORMAT(5X,2HC3/4(5X,4E10.3/)) 00666 2009 FORMAT(5X,2HC4/4(5X,4E10.3/)) 00667 2010 FORMAT(5X,4HAPA1/4(5X,4E10.3/)) 00668 2011 FORMAT(5X,1HZ,4E10.2/) 00669 END 00670 C 00671 SUBROUTINE EIGEN(A,R,N,KCHK) 00672 C COMPUTE EIGENVALUES AND EIGENVECTORS OF REAL SYMMETRIC MATRIX 00673 C INPUT ARGUMENT 00674 C A(N,N) = INPUT MATRIX (MAX 100X100) 00675 C OUTPUT ARGUMENTS 00676 C A(N,N) = MATRIX WHOSE DIAGONAL ELEMENTS ARE EIGENVALUES OF INPUT 00677 C R(N,N) = ORTHOGONAL MATRIX WHOSE COLUMNS ARE EIGENVECTORS OF INPUT 00678 C KCHK = 0 SUCCESSFUL RETURN 00679 C 1 NO CONVERGENCE 00680 C 2 A DIMENSIONED TOO LARGE 00681 C JACOBI METHOD IS USED - INPUT MATRIX ITERATIVELY REDUCED TO DIAGONAL 00682 C MATRIX BY SUCCESSIVELY ANNIHILATING OFF-DIAGONAL ELEMENTS THROUGH 00683 C SERIES OF ROTATIONS. R IS FINAL PRODUCT OF THESE ROTATIONS. 00684 C REFERENCE - CARNAHAN, LUTHER, WILKES 'APPLIED NUMERICAL METHODS' P250 00685 DIMENSION A(N,N),R(N,N),Y(100) 00686 C**MACHINE DEPENDENT STATEMENTS 00687 DOUBLE PRECISION A,R,Y,EPS1,EPS2,EPS3,SIGMA1,SIGMA2,OFFDSQ, 00688 * S,Q,SPQ,CSA,SNA,HOLDKI,ABS,SQRT,DABS,DSQRT,P,HOLDIK 00689 ABS(S) = DABS(S) 00690 SQRT(S) = DSQRT(S) 00691 C ITMAX = MAXIMUM NUMBER OF ITERATIONS 00692 C EPS1 = ROTATION ANGLE TEST (ROTATE BY PI/4 IF EQUAL DIAGONAL ELEMENTS)00693 C EPS2 = ROUNDOFF TEST (IGNORE OFFDIAGONAL ELEMENTS LT EPS2) 00694 C EPS3 = BREAKOUT TEST (SUM OF SQUARES OF DIAGONAL ELEMENTS UNCHANGED 00695 C BY LAST ITERATION) 00696 IY = 100 00697 ITMAX=50 00698 EPS1=0.1D-12 00699 EPS2=0.1D-12 00700 EPS3=0.1D-10 00701 NM1=N-1 00702 KCHK = 0 00703 IF(N .GT. IY) GO TO 30 00704 C SET R = IDENTITY AND UPPER TRIANGULATE A 00705 DO 16 I=1,N 00706 DO 16 J=1,N 00707 R(I,J)=0. 00708 IF(I.EQ.J) R(I,J)=1. 00709 IF(I.GT.J) A(I,J)=0. 00710 16 CONTINUE 00711 C COMPUTE SIGMA1 AND S 00712 SIGMA1=0. 00713 OFFDSQ=0. 00714 DO 17 I=1,N 00715 SIGMA1=SIGMA1+A(I,I)*A(I,I) 00716 IP1=I+1 00717 IF(I.GE.N) GOTO 18 00718 DO 17 J=IP1,N 00719 17 OFFDSQ = OFFDSQ + A(I,J) * A(I,J) 00720 18 S = 2D0 * OFFDSQ + SIGMA1 00721 C BEGIN JACOBI ITERATION 00722 DO 29 ITER=1,ITMAX 00723 DO 27 I=1,NM1 00724 IP1=I+1 00725 DO 27 J=IP1,N 00726 Q= ABS(A(I,I)-A(J,J)) 00727 C COMPUTE SINE AND COSINE OF ROTATION ANGLE 00728 IF(Q.LE.EPS1) GOTO 19 00729 IF( ABS(A(I,J)).LE.EPS2) GOTO 27 00730 P=2.0*A(I,J)*Q/(A(I,I)-A(J,J)) 00731 SPQ= SQRT(P*P+Q*Q) 00732 CSA= SQRT((1.0+Q/SPQ)/2.0) 00733 SNA=P/(2.0*CSA*SPQ) 00734 GOTO 20 00735 19 CSA = 1D0 / SQRT(2D0) 00736 SNA=CSA 00737 20 CONTINUE 00738 C UPDATE COLUMNS I AND J OF R (EQUIV TO MULT BY ANNIHILATION MATRIX) 00739 DO 21 K=1,N 00740 HOLDKI=R(K,I) 00741 R(K,I)=HOLDKI*CSA+R(K,J)*SNA 00742 21 R(K,J) = HOLDKI * SNA - R(K,J) * CSA 00743 C COMPUTE NEW ELEMENTS OF A IN ROWS I AND J 00744 DO 24 K=1,N 00745 IF(K.GT.J) GOTO 23 00746 Y(K)=A(I,K) 00747 A(I,K)=CSA*Y(K)+SNA*A(K,J) 00748 IF(K.NE.J) GOTO 22 00749 A(J,K)=SNA*Y(K)-CSA*A(J,K) 00750 22 GO TO 24 00751 23 HOLDIK = A(I,K) 00752 A(I,K)=CSA*HOLDIK+SNA*A(J,K) 00753 A(J,K)=SNA*HOLDIK-CSA*A(J,K) 00754 24 CONTINUE 00755 C COMPUTE NEW ELEMENTS OF A IN COLUMNS I AND J 00756 Y(J)=SNA*Y(I)-CSA*Y(J) 00757 C WHEN K IS LARGER THAN I 00758 DO 26 K=1,J 00759 IF(K.LE.I) GOTO 25 00760 A(K,J)=SNA*Y(K)-CSA*A(K,J) 00761 GOTO 26 00762 25 HOLDKI = A(K,I) 00763 A(K,I)=CSA*HOLDKI+SNA*A(K,J) 00764 A(K,J)=SNA*HOLDKI-CSA*A(K,J) 00765 26 CONTINUE 00766 27 A(I,J) = 0. 00767 C FIND SIGMA2 FOR TRANSFORMED A AND TEST FOR CONVERGENCE 00768 SIGMA2=0. 00769 DO 28 I=1,N 00770 28 SIGMA2 = SIGMA2 + A(I,I) * A(I,I) 00771 IF(1.0-SIGMA1/SIGMA2.LT.EPS3) RETURN 00772 29 SIGMA1 = SIGMA2 00773 C ERROR RETURNS 00774 WRITE(6,1001) ITER,S,SIGMA1,SIGMA2 00775 KCHK = 1 00776 RETURN 00777 30 WRITE(6,1002) 00778 KCHK = 2 00779 RETURN 00780 C FORMAT STATEMENTS 00781 1001 FORMAT(20H0EIGEN NOT CONVERGED /1X,6HITER = ,I5/1X,3HS = ,F10.5/ 00782 * 1X,7HSIGMA1= ,F10.5/1X,7HSIGMA2= ,F10.5) 00783 1002 FORMAT(41H EIGEN INPUT MATRIX DIMENSIONED TOO LARGE ) 00784 END 00785 SUBROUTINE CHOLD(A,IRDA,NA,DETA,JCHK) 00786 C MATRIX INVERSION USING CHOLESKI DECOMPOSITION 00787 C INPUT ARGUMENTS 00788 C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX 00789 C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX 00790 C NA = SIZE OF INPUT MATRIX 00791 C OUTPUT ARGUMENTS 00792 C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED) 00793 C DETA = DETERMINANT OF INPUT MATRIX 00794 C JCHK = 0 SUCCESSFUL RETURN 00795 C 1 SINGULAR MATRIX 00796 C 2 MATRIX OF ZERO DIMENSION 00797 DIMENSION A(IRDA,NA) 00798 C**MACHINE DEPENDENT STATEMENTS 00799 DOUBLE PRECISION A,DETA,SUM,SQRT,ABS,SING,DSQRT,DABS 00800 SQRT(DETA) = DSQRT(DETA) 00801 ABS(DETA) = DABS(DETA) 00802 SING = 1D-20 00803 JCHK = 0 00804 C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 00805 IF(NA .LT. 1) GO TO 18 00806 DETA = A(1,1) 00807 A(1,1) = SQRT(A(1,1)) 00808 IF(NA .EQ. 1) GO TO 6 00809 DO 1 I = 2,NA 00810 1 A(I,1) = A(I,1) / A(1,1) 00811 DO 5 J = 2,NA 00812 SUM = 0. 00813 J1 = J - 1 00814 DO 2 K = 1,J1 00815 2 SUM = SUM + A(J,K) ** 2 00816 DETA = DETA * (A(J,J) - SUM) 00817 A(J,J) = SQRT(A(J,J) - SUM) 00818 IF(J .EQ. NA) GO TO 5 00819 J2 = J + 1 00820 DO 4 I = J2,NA 00821 SUM = 0. 00822 DO 3 K = 1,J1 00823 3 SUM = SUM + A(I,K) * A(J,K) 00824 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 00825 5 CONTINUE 00826 6 IF(ABS(DETA) .LT. SING) GO TO 16 00827 C INVERSION OF LOWER TRIANGULAR MATRIX 00828 DO 7 I = 1,NA 00829 7 A(I,I) = 1. / A(I,I) 00830 IF(NA .EQ. 1) GO TO 10 00831 N1 = NA - 1 00832 DO 9 J = 1,N1 00833 J2 = J + 1 00834 DO 9 I = J2,NA 00835 SUM = 0. 00836 I1 = I - 1 00837 DO 8 K = J,I1 00838 8 SUM = SUM + A(I,K) * A(K,J) 00839 9 A(I,J) = - A(I,I) * SUM 00840 C CONSTRUCTION OF INVERSE OF INPUT MATRIX 00841 10 DO 15 J = 1,NA 00842 IF(J .EQ. 1) GO TO 12 00843 J1 = J - 1 00844 DO 11 I = 1,J1 00845 11 A(I,J) = A(J,I) 00846 12 DO 14 I = J,NA 00847 SUM = 0. 00848 DO 13 K = I,NA 00849 13 SUM = SUM + A(K,I) * A(K,J) 00850 14 A(I,J) = SUM 00851 15 CONTINUE 00852 RETURN 00853 C ERROR RETURNS 00854 16 WRITE(6,17) DETA 00855 JCHK = 1 00856 RETURN 00857 18 WRITE(6,19) 00858 JCHK = 2 00859 RETURN 00860 C FORMAT STATEMENTS 00861 17 FORMAT(33H SINGULAR MATRIX IN CHOLD. DET = ,E20.5) 00862 19 FORMAT(35H MATRIX OF DIMENSION ZERO IN CHOLD ) 00863 END 00864 C 00865 FUNCTION CHISQ(A,N) 00866 C PERCENTILES OF THE CHI-SQUARE DISTRIBUTION X(N) 00867 C INPUT 00868 C A = PROBABILITY INTEGRAL FROM ZERO TO CHISQ 00869 C N = NUMBER OF DEGREES OF FREEDOM 00870 C OUTPUT 00871 C CHISQ = ABSCISSA VALUE OF X(N) CORRESPONDING TO PROBABILITY A 00872 C ABRAMOWITZ AND STEGUN EQN 26.4.17 IS USED (WILSON-HILFERTY APPROX) 00873 C ACCURACY BETTER THAN 0.04 FOR N .GT. 1 00874 C**MACHINE DEPENDENT STATEMENTS 00875 DOUBLE PRECISION A,CHISQ,XNORM,SQRT,DSQRT 00876 SQRT(A) = DSQRT(A) 00877 CHISQ =N*(1-2./(9.*N)+XNORM(A,0D0,1D0)* SQRT(2D0/(9.*N)))**3 00878 RETURN 00879 END 00880 C 00881 FUNCTION XNORM(ALF,XMEAN,SIG) 00882 C PERCENTILES OF THE NORMAL DISTRIBUTION N(XMEAN,SIG) 00883 C INPUT 00884 C ALF = PROBABILITY INTEGRAL FROM NEGATIVE INFINITY TO XNORM 00885 C XMEAN = POPULATION MEAN 00886 C SIG = POPULATION STANDARD DEVIATION 00887 C OUTPUT 00888 C XNORM = ABSCISSA VALUE OF N(XMEAN,SIG) CORRESPONDING TO 00889 C PROBABILITY ALF 00890 C ABRAMOWITZ AND STEGUN EQUATION 26.2.23 IS USED 00891 C ACCURACY BETTER THAN THAN 0.00045 00892 DIMENSION C(6) 00893 C**MACHINE DEPENDENT STATEMENTS 00894 DOUBLE PRECISION XNORM,ALF,XMEAN,SIG,P,T,T2,T3,XP,SQRT,ALOG,DSQRT,00895 * DLOG 00896 SQRT(P) = DSQRT(P) 00897 ALOG(P) = DLOG(P) 00898 DATA C(1),C(2),C(3),C(4),C(5),C(6)/2.515517,.802853,.010328, 00899 * 1.432788,.189269,.001308/ 00900 IF(ALF .GE. 1. .OR. ALF .LE. 0.) GO TO 20 00901 SIGN = -1. 00902 P = ALF 00903 IF(ALF .LT. .5) GO TO 10 00904 SIGN = 1. 00905 P = 1. - ALF 00906 10 T = SQRT(ALOG(1./P**2)) 00907 T2 = T * T 00908 T3 = T * T2 00909 XP = T - (C(1)+C(2)*T+C(3)*T2) / (1.+C(4)*T+C(5)*T2+C(6)*T3) 00910 XP = XP * SIGN 00911 XNORM = XMEAN + XP * SIG 00912 RETURN 00913 C ERROR RETURN 00914 20 WRITE(6,1001) ALF 00915 XNORM = 0. 00916 RETURN 00917 C FORMAT STATEMENTS 00918 1001 FORMAT(25H XNORM INPUT PROBABILITY= ,E20.10) 00919 END 00920 C 00921 C 00922 SUBROUTINE ROTREF(NUM,NAXIS,ANGLE,ROT) 00923 C COMPUTE PRODUCT MATRIX OF A SEQUENCE OF ROTATIONS AND REFLECTIONS 00924 C INPUT ARGUMENTS 00925 C NUM = NUMBER OF ROTATIONS AND REFLECTIONS IN SEQUENCE 00926 C NAXIS(NUM) = SEQUENCE OF ROTATION AND REFLECTION AXES 00927 C FOR ROTATIONS USE 1,2, OR 3 00928 C FOR REFLECTIONS USE -1, -2, OR -3 00929 C ANGLE(NUM) = SEQUENCE OF ROTATION ANGLES IN RADIANS 00930 C FOR REFLECTIONS THIS ANGLE IS IGNORED (SET TO ZERO) 00931 C OUTPUT ARGUMENT 00932 C ROT(3,3) = PRODUCT MATRIX 00933 DIMENSION ROT(3,3),R1(3,3),R2(3),ANGLE(NUM),NAXIS(NUM) 00934 C**MACHINE DEPENDENT STATEMENTS 00935 DOUBLE PRECISION ROT,R1,R2,ANGLE,EPS,COS,SIN,ABS,DCOS,DSIN,DABS 00936 COS(EPS) = DCOS(EPS) 00937 SIN(EPS) = DSIN(EPS) 00938 ABS(EPS) = DABS(EPS) 00939 EPS = 1D-15 00940 C SET ROT = IDENTITY MATRIX 00941 DO 1 I = 1,3 00942 DO 1 J = 1,3 00943 ROT(I,J) = 0. 00944 IF(I .EQ. J) ROT(I,J) = 1. 00945 1 CONTINUE 00946 C CHECK ELEMENTS OF NAXIS AND SET REFLECTION ELEMENTS OF ANGLE = 0 00947 DO 2 N = 1,NUM 00948 IF(NAXIS(N) .EQ. 0 .OR. NAXIS(N) .LT. -3 .OR. NAXIS(N) .GT. 3) 00949 * GO TO 5 00950 IF(NAXIS(N) .LT. 0) ANGLE(N) = 0. 00951 2 CONTINUE 00952 C PROCESS SEQUENCE OF ROTATIONS AND REFLECTIONS ONE AT A TIME 00953 DO 4 N = 1,NUM 00954 C DEFINE THREE AXES FOR CURRENT ROTATION OR REFLECTION 00955 N1 = IABS(NAXIS(N)) 00956 N2 = MOD(N1,3) + 1 00957 N3 = MOD(N2,3) + 1 00958 C DEFINE DIAGONAL ELEMENTS 00959 R1(N1,N1) = 1. 00960 IF(NAXIS(N) .LT. 0.) R1(N1,N1) = -1. 00961 R1(N2,N2) = COS(ANGLE(N)) 00962 R1(N3,N3) = R1(N2,N2) 00963 C DEFINE NON-ZERO OFF-DIAGONAL ELEMENTS 00964 R1(N2,N3) = SIN(ANGLE(N)) 00965 R1(N3,N2) = - R1(N2,N3) 00966 C DEFINE ZERO OFF-DIAGONAL ELEMENTS 00967 R1(N1,N2) = 0. 00968 R1(N1,N3) = 0. 00969 R1(N2,N1) = 0. 00970 R1(N3,N1) = 0. 00971 C FORM PRODUCT ( SET ROT = R1 * ROT) 00972 DO 4 J = 1,3 00973 DO 3 I = 1,3 00974 R2(I) = 0. 00975 DO 3 K = 1,3 00976 3 R2(I) = R2(I) + R1(I,K) * ROT(K,J) 00977 DO 4 I = 1,3 00978 ROT(I,J) = R2(I) 00979 IF(ABS(ROT(I,J)) .LT. EPS) ROT(I,J) = 0. 00980 4 CONTINUE 00981 RETURN 00982 C ERROR RETURN 00983 5 WRITE(6,6) N,NAXIS(N) 00984 RETURN 00985 C FORMAT STATEMENTS 00986 6 FORMAT(21H ILLEGAL VALUE NAXIS( ,I3,4H) = ,I5,10H IN ROTREF ) 00987 END 00988 SUBROUTINE RECGEO(A,B,X0,R0,S0,XG,SXYZ,PLH,SPLH) 00989 C AVERAGE TERRESTRIAL TO GEODETIC CONVERSION 00990 C INPUT AVERAGE TERRESTRIAL CARTESIAN COORDINATES IN M 00991 C SEVEN PARAMETER TRANSFORMATION PERFORMED 00992 C OUTPUT GEODETIC COORDINATES LAT/LONG/HT IN DEGREES AND M 00993 C REFERRED TO ELLIPSOID WITH AXES A,B 00994 C ITERATIVE COMPUTATION CONTINUED UNTIL CONVERGENCE TO 1E-10 RADIANS IN 00995 C PHI AND A*1E-10 METRES IN HT, OR FOR MAXIMUM OF 20 ITERATIONS 00996 C**MACHINE DEPENDENT STATEMENT 00997 IMPLICIT REAL*8(A-H,O-Z) 00998 COMMON /DSRN/ MET,IN,IOUT,IOUT2,ICR,IPR,IPN,IPR2,IPR3,ILIST 00999 COMMON /CNST/ PI,VEL,WE,FO,RADG,WAVENO,BITS,BITR,NPASS,NULL,IEDS 01000 DIMENSION X(3),XX(3),XG(3),X0(3),R0(3),PLH(3),SPLH(3,3),SXYZ(3,3),01001 $ R(3),ROT(3,3),NAX1(3),NAX2(3) 01002 DATA NAX1(1),NAX1(2),NAX1(3)/3,2,1/ 01003 DATA NAX2(1),NAX2(2),NAX2(3)/3,2,-2/ 01004 ABS(PI) = DABS(PI) 01005 SQRT(PI) = DSQRT(PI) 01006 SIN(PI) = DSIN(PI) 01007 COS(PI) = DCOS(PI) 01008 ATAN(PI) = DATAN(PI) 01009 ATAN2(X1,Y1) = DATAN2(X1,Y1) 01010 IF(ILIST .GE. 1) WRITE(IPR,1000) 01011 C INITIALIZE 01012 EPS2 = 1E-10 01013 EPS1 = A * EPS2 01014 ITER = 20 01015 BOA = (B/A)**2 01016 IF(B .LT. 6D6) BOA = (1. - 1. / B) **2 01017 E2=1.-BOA 01018 C THREE TRANSLATIONS 01019 DO 10 I = 1,3 01020 10 X(I) = XG(I) - X0(I) 01021 C SCALE CHANGE 01022 DO 20 I = 1,3 01023 20 X(I) = X(I) / (1D0 + S0 * 1D-6) 01024 C THREE ROTATIONS 01025 DO 30 I = 1,3 01026 30 R(I) = R0(I) * RADG / 3600. 01027 CALL ROTREF(3,NAX1,R,ROT) 01028 DO 40 I = 1,3 01029 XX(I) = 0D0 01030 DO 40 J = 1,3 01031 40 XX(I) = XX(I) + ROT(I,J) * X(J) 01032 X1 = XX(1) 01033 Y1 = XX(2) 01034 Z1 = XX(3) 01035 C COMPUTE LONGITUDE NON-ITERATIVELY 01036 P = SQRT(X1**2 + Y1**2) 01037 XLAMB = ATAN2(Y1,X1) 01038 C INITIAL APPROXIMATIONS FOR ITERATION 01039 XN=A 01040 HT = SQRT(X1**2 + Y1**2 + Z1**2) - A 01041 PHI = ATAN((Z1/P)/(1.-E2*XN/(XN+HT))) 01042 C ITERATION LOOP 01043 DO 50 I = 1,ITER 01044 HT1=HT 01045 PHI1=PHI 01046 XN = A / SQRT(COS(PHI)**2 + BOA * SIN(PHI)**2) 01047 HT = P / COS(PHI) - XN 01048 PHI = ATAN((Z1/P)/(1.-E2*XN/(XN+HT))) 01049 IF(ABS(HT1-HT) .LT. EPS1 .AND. ABS(PHI1-PHI) .LT. EPS2) GO TO 60 01050 50 CONTINUE 01051 GO TO 70 01052 60 PLH(1) = PHI / RADG 01053 PLH(2) = XLAMB / RADG 01054 PLH(3) = HT 01055 C TRANSFORM COVARIANCE MATRIX 01056 R(1) = XLAMB - PI 01057 R(2) = PHI - PI / 2. 01058 R(3) = 0. 01059 CALL ROTREF(3,NAX2,R,ROT) 01060 CALL PROP(SPLH,ROT,SXYZ,3,3) 01061 RETURN 01062 C ERROR RETURN 01063 70 WRITE(IPR,1001) 01064 STOP 01065 C FORMAT STATEMENTS 01066 1000 FORMAT(1X,30HSUBROUTINE RECGEO ENTERED ) 01067 1001 FORMAT(1X,30HTOO MANY ITERATIONS IN RECGEO ) 01068 END 01069 SUBROUTINE PROP(QX,T,QY,NX,NY) 01070 C PROPAGATION OF COVARIANCE MATRICES 01071 C IF X = T Y 01072 C THEN QX = T QY T(T) 01073 C**MACHINE DEPENDENT STATEMENT 01074 IMPLICIT REAL*8(A-H,O-Z) 01075 DIMENSION QX(NX,NX),T(NX,NY),QY(NY,NY) 01076 DO 10 I = 1,NX 01077 DO 10 J = 1,NX 01078 QX(I,J) = 0D0 01079 DO 10 K = 1,NY 01080 DO 10 L = 1,NY 01081 10 QX(I,J) = QX(I,J) + T(I,K) * T(J,L) * QY(K,L) 01082 RETURN 01083 END 01084