FTN MODELS C ADJUSTMENT OF MODELS SEPTEMBER 1980 G.H.S. 7001 C 7002 C ******* ADJUSTMENT ******* MAIN ROUTINE ******* 7003 C 7004 C *************************************************************** 7005 C * SET DIMENSIONS AND 8 CONSTANTS AS FOLLOWS, OR LARGER * 7006 C * MMO NUMBER OF MODELS IN BLOCK * 7007 C * MIM MAXIMUM NUMBER OF MODEL POINTS FOR ONE TERRAIN POINT * 7008 C * OR ONE LAKE LEVEL * 7009 C * MW BANDWIDTH * 7010 C * MY NUMBER OF ADDED PARAMETERS, BUT NOT SMALLER THAN 1 * 7011 C * MZ NUMBER OF ADDED PARAMETERS + 1 * 7012 C * MP IF MIM .LE. MW: 7 * MW + MZ * 7013 C * MP IF MIM .GT. MW: 7 * MIM + MZ * 7014 C * MQ (MW * (MW + 1)) / 2 * 7015 C * MH MAXIMUM NUMBER OF SPECIAL HEIGHT CONTROL POINTS * 7016 C * LORI(8,MMO),V(7,MMO),LT(4,MIM),RES(3,MIM),C(3,MZ),LSPH(2,MH)* 7017 C * P(3,MP),Q(7,7,MQ),S(7,MZ,MW),Z(MY,MZ),ROT(3,4,MW) * 7018 C *************************************************************** 7019 C 7020 COMMON FLP(11), N1,N2,N3,ITER,ITR,NEWIT,LIST,JD31,JD32,JD33,IPR, 7021 1 NL,LAX(2,40),LTER(4),JL,ILK,K,KK,KDIF,INC,KKMAX,KAUXL,KLAKE, 7022 2 LLAC, N,NIM,NW,NZ,NY,NP,NQ,NH, IL,IM,IT 7023 EQUIVALENCE (P(1,1),V(1,1)) 7024 C *************************************************************** 7025 C DIMENSIONS FOR CRYSLER TEST BLOCK ADJUSTMENT OF MODELS 7026 C ********************************* 7027 DIMENSION LORI(8, 26),V(7, 26),LT(4, 6),RES(3, 6), 7028 1 C(3, 7),LSPH(2, 1) 7029 DOUBLE PRECISION P(3,49),Q(7,7,21),S(7,7, 6),Z( 6, 7),ROT(3,4, 6) 7030 DATA MMO,MIM,MW, MY,MZ, MP,MQ, MH / 26,6,6, 6,7, 49,21, 1 / 7031 C *************************************************************** 7032 C 7033 ITR = 0 7034 100 CALL MODL1 (LORI,LSPH,Q,S,Z,C, MMO,MIM,MW,MY,MZ,MP,MQ,MH) 7035 C 7036 C FORM THE NORMAL EQUATIONS, ELIMINATE THE TERRAIN POINT CORRECTIONS 7037 200 CALL MODL2 (LORI,LSPH, ROT,LT,RES,C,S,Z,P,Q, N,NIM,NW, 7038 1 NY,NZ,NP,NQ,NH) 7039 C 7040 C SOLVE THE REDUCED NORMAL EQUATIONS, CORRECT THE ORIENTATIONS 7041 CALL MODL3 (Q,S,Z, NW,NY,NZ,NQ) 7042 IF (K .GT. 0) GO TO 200 7043 CALL MODL4 (Q,S,Z,LORI,V, N,NW,NY,NZ,NQ) 7044 C 7045 C CORRECT THE TERRAIN POINT COORDINATES 7046 CALL MODL5 (LORI,ROT,LT,RES,Z,C, N,NIM,NW,NY,NZ) 7047 IF (NEWIT .GT. 0) GO TO 100 7048 STOP 77 7049 END 7050 FTN SUBM1 SUBROUTINE MODL1 (LORI,LSPH,Q,S,Z,C, MMO,MIM,MW,MY,MZ,MP,MQ,MH) 7201 C 7202 C READ CONTROL CARDS AND RECORD WITH APPROX. ORIENTATION PARAMETERS 7203 C 7204 COMMON WT(9),RT,CURV, N1,N2,N3,ITER,ITR,NEWIT,LIST,JD31,JD32,JD33, 7205 1 IPR,NL,LAX(2,40),LTER(4),JKI(6),KKMAX,KAUXL,KLAKE,LLAC, 7206 2 N,NIM,NW,NZ,NY,NP,NQ,NH 7207 DIMENSION A(19), LORI(8,MMO), LSPH(2,MH), C(3,MZ), WEQ(9) 7208 DOUBLE PRECISION Q(7,7,MQ), S(7,MZ,MW), Z(MY,MZ) 7209 EQUIVALENCE (W1EN,WEQ(1)),(W1H,WEQ(2)),(W2EN,WEQ(3)), 7210 1 (W2H,WEQ(4)),(WSPH,WEQ(5)),(WMODZ,WEQ(6)), 7211 2 (WPC12,WEQ(7)),(WPC3,WEQ(8)),(WLAKE,WEQ(9)) 7212 1 FORMAT (19A4, I4 / 8I4, 12X, 9F4.2) 7213 4 FORMAT (I10, 20X, I10) 7214 8 FORMAT (10I8) 7215 10 FORMAT (1H1 19A4 //I6,14H ITERATION(S),I3,19H ADDED PARAMETER(S)) 7216 11 FORMAT (30H0 EARTH CURVATURE CORRECTION) 7217 12 FORMAT (25H0 REFRACTION CORRECTION) 7218 13 FORMAT (11H0 RELEASE (1H T12, 10I8)) 7219 14 FORMAT ( 1H0 I8,' SPECIAL HEIGHTS') 7220 15 FORMAT (11H0 WEIGHTS 5F8.2, 4X, 4F8.2) 7221 16 FORMAT (45H0 INITIAL TERRAIN COORDINATES AND RESIDUALS) 7222 18 FORMAT (28H0 AUXILIARY CONTROL POINTS) 7223 19 FORMAT (21H0 LAKE LEVEL POINTS) 7224 21 FORMAT (13H0 MORE THAN I3,' RELEASED POINTS') 7225 22 FORMAT (13H0 MORE THAN I4,' SPECIAL HEIGHTS') 7226 23 FORMAT (19H0 NO DATA IN FILE I3) 7227 24 FORMAT (55H0 THE VARIABLE DIMENSIONS MMO MIM MW MY MZ MP MQ 7228 1 / 17H ARE PRESENTLY T28, 7I4 7229 2 / 20H BUT MUST BE .GE. T28, 7I4) 7230 IF (ITR .GT. 0) GO TO 1025 7231 C 7232 C DEFINE INPUT-OUTPUT DEVICES 7233 IRD = 1 7234 IPR = 3 7235 JD31 = 31 7236 JD32 = 32 7237 JD33 = 33 7238 C 7239 C READ THE CONTROL CARDS 7240 N1 = 3 7241 N2 = 7 7242 C NLMAX IS MAXIMUM NUMBER OF RELEASED POINTS, SEE ARRAY LAX 7243 NLMAX = 40 7244 READ (IRD,1) A, LIST, NEWIT, N3, KECC,KREF, NL,NH, KAUXL,KLAKE,WT 7245 WRITE (IPR,10) A, NEWIT, N3 7246 ERAD = 6378.E5 7247 CURV = 0. 7248 IF (KECC .LT. 0) GO TO 1001 7249 WRITE (IPR,11) 7250 CURV = 1. / ERAD 7251 1001 IF (KREF .LT. 0) GO TO 1010 7252 WRITE (IPR,12) 7253 CURV = CURV - .1 / ERAD 7254 C 7255 C READ RECORD WITH ORIENTATION PARAMETERS 7256 1010 READ (JD31,END=1053) ITER,N,NIM,NW,LLAC, ((LORI(J1,J2),J1=1,8), 7257 1 J2=1,N) 7258 IF (KAUXL.GE.0 .AND. LLAC.GE.2) WRITE (IPR,18) 7259 IF (KLAKE.GE.0 .AND. (LLAC.EQ.1 .OR. LLAC.EQ.3)) WRITE (IPR,19) 7260 NZ = N3 + 1 7261 NY = N3 7262 IF (NY .EQ. 0) NY = 1 7263 NP = N2 * NIM + NZ 7264 NQ = (NW * (NW + 1)) / 2 7265 IF (N.GT.MMO.OR.NIM.GT.MIM.OR.NW.GT.MW) GO TO 1054 7266 IF (NY.GT.MY.OR.NZ.GT.MZ.OR.NP.GT.MP.OR.NQ.GT.MQ) GO TO 1054 7267 KKMAX = N - NW + 1 7268 C READ LIST OF POINTS TO BE RELAXED 7269 IF (NL .LE. 0) GO TO 1011 7270 IF (NL .GT. NLMAX) GO TO 1051 7271 READ (IRD,8) ((LAX(J1,J2), J1=1,2), J2=1,NL) 7272 WRITE (IPR,13) ((LAX(J1,J2), J1=1,2), J2=1,NL) 7273 C READ SPECIAL HEIGHTS 7274 1011 IF (NH .LE. 0) GO TO 1020 7275 DO 1012 J2 = 1,NH 7276 IF (J2 .GT. MH) GO TO 1052 7277 READ (IRD,4, END=1013) (LSPH(J,J2), J=1,2) 7278 IF (LSPH(1,J2) .EQ. 0) GO TO 1013 7279 1012 CONTINUE 7280 GO TO 1014 7281 1013 NH = J2 - 1 7282 1014 WRITE (IPR,14) NH 7283 C 7284 C WEIGHTS ATTACHED TO GROUND CONTROL 7285 1020 W1EN = 25. 7286 W1H = 25. 7287 W2EN = 2. 7288 W2H = 2. 7289 WSPH = .025 7290 C WEIGHTS ATTACHED TO CORRECTION EQUATIONS 7291 WMODZ = 1. 7292 WPC12 = .25 7293 WPC3 = .25 7294 WLAKE = 1. 7295 J1 = 0 7296 DO 1022 J = 1,9 7297 IF (WT(J) .GT. 0.) GO TO 1021 7298 WT(J) = WEQ(J) 7299 GO TO 1022 7300 1021 J1 = 1 7301 1022 CONTINUE 7302 IF (J1 .EQ. 1) WRITE (IPR,15) WT 7303 IF (LIST.GT.0 .OR. NEWIT.LE.0) WRITE (IPR,16) 7304 C ROOTS OF WEIGHTS 7305 DO 1023 J = 6,9 7306 IF (WT(J) .NE. 1.) WT(J) = WT(J)**.5 7307 1023 CONTINUE 7308 GO TO 1030 7309 C 7310 1025 READ (JD31) 7311 C 7312 C ZERO THE NORMAL EQUATIONS 7313 1030 DO 1031 J = 1,NQ 7314 DO 1031 J2 = 1,N2 7315 DO 1031 J1 = 1,N2 7316 1031 Q(J1,J2,J) = 0. 7317 C NOTE: DO NOT REPLACE MZ AND MY BY NZ AND NY 7318 DO 1032 J = 1,NW 7319 DO 1032 J2 = 1,MZ 7320 DO 1032 J1 = 1,N2 7321 1032 S(J1,J2,J) = 0. 7322 IF (NZ .EQ. 1) GO TO 1040 7323 DO 1034 J2 = 1,NZ 7324 C(3,J2) = 0. 7325 DO 1034 J1 = 1,MY 7326 1034 Z(J1,J2) = 0. 7327 C ZERO K, KK, KDIF, INC 7328 1040 DO 1041 J = 3,6 7329 1041 JKI(J) = 0 7330 RETURN 7331 C 7332 1051 WRITE (IPR,21) NLMAX 7333 GO TO 1060 7334 1052 WRITE (IPR,22) MH 7335 GO TO 1060 7336 1053 WRITE (IPR,23) JD31 7337 GO TO 1060 7338 1054 WRITE (IPR,24) MMO,MIM,MW,MY,MZ,MP,MQ, N,NIM,NW,NY,NZ,NP,NQ 7339 1060 STOP 1 7340 END 7341 FTN SUBM2 SUBROUTINE MODL2 (LORI,LSPH, ROT,LT,RES,C,S,Z,P,Q, N,NIM,NW, 7401 1 NY,NZ,NP,NQ,NH) 7402 C 7403 C FORM NORMAL EQUATIONS AND ELIMINATE TERRAIN POINT CORRECTIONS 7404 C 7405 COMMON WT(9),ROOTW,CU, N1,N2,N3,ITER,ITR,NEWIT,LIST,JD31,JD32, 7406 1 JD33,IPR, NL,LAX(2,40),LTER(4),JL,ILK,K,KK,KDIF,INC,KKMAX, 7407 2 KAUXL,KLAKE, INT(9),IL,IM,IT 7408 DIMENSION LORI(8,N),LT(4,NIM),C(3,NZ),RES(3,NIM),LSPH(2,NH) 7409 DOUBLE PRECISION U(3),P(3,NP),Q(7,7,NQ),S(7,NZ,NW),Z(NY,NZ), 7410 1 ROT(3,4,NW), ABCD(4),AR,BR,CR,DR 7411 EQUIVALENCE (AR,ABCD(1)),(BR,ABCD(2)),(CR,ABCD(3)),(DR,ABCD(4)) 7412 28 FORMAT (1H I3,I7, 2I10,I8, (1H T40, 4(F9.0,2F7.0))) 7413 29 FORMAT (1H I3,I7, 2I10,I8, (1H T40, 4(16X, F7.0))) 7414 IF (K .GT. 0) GO TO 1110 7415 C 7416 C READ A RECORD, FORM THE NN ORIENTATION MATRICES 7417 1105 READ (JD31,END=1150) LTER, ILK, NN, ((LT(J,J1),J=1,4),J1=1,NN) 7418 IF (K .EQ. LT(1,1)) GO TO 1120 7419 K = LT(1,1) 7420 IF (K .GT. KKMAX) INC = K - KKMAX 7421 IF (KK .EQ. KKMAX) GO TO 1120 7422 IF (KK .GT. 0) GO TO 1107 7423 KK = K 7424 L1 = 1 7425 GO TO 1115 7426 C 7427 C TO ELIMINATE ONE OR MORE SETS OF N2 PARAMETERS, GO TO MODL3 7428 1107 KDIF = K - KK 7429 IF (K .GT. KKMAX) KDIF = KKMAX - KK 7430 IF (NEWIT) 1110, 1110, 1151 7431 C 7432 C UPON RETURN FROM MODL3, SHIFT THE AVAILABLE ORIENTATION MATRICES 7433 1110 KK = KK + KDIF 7434 M1 = KDIF + 1 7435 DO 1112 J2 = M1,NW 7436 M2 = J2 - KDIF 7437 DO 1112 J1 = 1,3 7438 DO 1112 J = 1,3 7439 1112 ROT(J,J1,M2) = ROT(J,J1,J2) 7440 L1 = NW - KDIF + 1 7441 C 7442 C FORM THE MISSING MATRICES 7443 1115 DO 1118 L = L1,NW 7444 J = L + KK - 1 7445 DO 1116 J1 = 1,4 7446 1116 ABCD(J1) = LORI(J1+4,J) * 1D-7 7447 ROT(2,3,L) = 2D0*(BR*CR - AR*DR) 7448 ROT(3,1,L) = 2D0*(CR*AR - BR*DR) 7449 ROT(1,2,L) = 2D0*(AR*BR - CR*DR) 7450 ROT(3,2,L) = 2D0*(BR*CR + AR*DR) 7451 ROT(1,3,L) = 2D0*(CR*AR + BR*DR) 7452 ROT(2,1,L) = 2D0*(AR*BR + CR*DR) 7453 DO 1117 J1 = 1,4 7454 1117 ABCD(J1) = ABCD(J1)**2 7455 ROT(1,1,L) = DR + AR - BR - CR 7456 ROT(2,2,L) = DR - AR + BR - CR 7457 ROT(3,3,L) = DR - AR - BR + CR 7458 ROT(3,4,L) = (ROT(1,1,L)**2 + ROT(2,1,L)**2)**.5 7459 DO 1118 J1 = 1,2 7460 1118 ROT(J1,4,L) = ROT(J1,1,L) / ROT(3,4,L) 7461 C 7462 C CHECK POINT AGAINST LIST OF SPECIAL HEIGHTS 7463 1120 IF (ITR .GT. 0) GO TO 1124 7464 C ERASE ANY PREVIOUS SPECIAL HEIGHTS, INSERT NEW ONES 7465 IF (ILK .GT. 10) ILK = ILK - 2 7466 IF (ILK .NE. 7) GO TO 1121 7467 ILK = 0 7468 IF (NN .GT. 1) ILK = 8 7469 1121 IF (NH .EQ. 0) GO TO 1124 7470 DO 1123 J = 1,NH 7471 IF (LTER(1) .NE. LSPH(1,J)) GO TO 1123 7472 LTER(4) = LSPH(2,J) 7473 IF (ILK .GE. 9) GO TO 1122 7474 ILK = 7 7475 GO TO 1124 7476 1122 ILK = ILK + 2 7477 GO TO 1124 7478 1123 CONTINUE 7479 C CHECK POINT AGAINST LIST OF RELAXED POINTS 7480 1124 IF (NL .EQ. 0) GO TO 1126 7481 DO 1125 JL = 1,NL 7482 IF (LTER(1) .EQ. LAX(1,JL)) GO TO 1127 7483 1125 CONTINUE 7484 1126 JL = 0 7485 GO TO 1130 7486 1127 IF (LAX(2,JL)) 1105, 1128, 1130 7487 1128 IF (NN .GT. 1) GO TO 1130 7488 LAX(2,JL) = 0 7489 1129 ILK = 0 7490 C 7491 C FORM CORR EQ AND CONTRIBUTIONS TO NORMAL EQ, 7492 1130 IF (ILK.EQ.0 .OR. (ILK.EQ.9 .AND. NN.EQ.1)) GO TO 1141 7493 IF (ILK.EQ.10 .AND. KLAKE.LT.0) GO TO 1141 7494 IF (ILK.GE.4 .AND.ILK.LE.6 .AND.NN.EQ.1 .AND.KAUXL.LT.0)GO TO 1141 7495 ND = NN * N2 + NZ 7496 IL = 1 7497 IM = 3 7498 ROOTW = WT(6) 7499 C DEAL WITH PARTIAL CONTROL POINTS, LAKE LEVELS POINTS, PROJ CENTRES 7500 IF (NN .GT. 1) GO TO 1131 7501 IF (ILK.EQ.1 .OR. ILK.EQ.4 .OR. ILK.GE.7) IL = 3 7502 IF (ILK.EQ.2 .OR. ILK.EQ.5) IM = 2 7503 GO TO 1133 7504 1131 IF (ILK.NE.10 .AND. ILK.NE.12) GO TO 1133 7505 ND = N2 + NZ 7506 DO 1132 J = 2,NN 7507 IF (LT(1,J) .NE. LT(1,J-1)) ND = ND + N2 7508 1132 CONTINUE 7509 IL = 3 7510 ROOTW = WT(9) 7511 1133 IT = 3 7512 IF (ILK.EQ.9 .OR. ILK.EQ.11) IT = 1 7513 DO 1137 J1 = 1,3 7514 U(J1) = 0. 7515 DO 1137 J = 1,ND 7516 1137 P(J1,J) = 0. 7517 CALL MOD2A (LORI,LT,RES,ROT,S,Z,U,P,Q,C,N,NW,NY,NZ,NQ,ND,NN) 7518 IF (ITR .GT. 0) GO TO 1140 7519 IF (LIST.LE.0 .AND. NEWIT.GT.0) GO TO 1140 7520 C LIST RESIDUALS WITH RESPECT TO INITIAL VALUES 7521 IF (ILK.LT.10 .OR. ILK.EQ.11) GO TO 1138 7522 WRITE (IPR,29) ILK, LTER, (RES(3,J1), J1=1,NN) 7523 GO TO 1139 7524 1138 WRITE (IPR,28) ILK, LTER, ((RES(J,J1), J=1,3), J1=1,NN) 7525 1139 IF (NEWIT .LE. 0) GO TO 1105 7526 C 7527 C ELIMINATE THE COORDINATE CORRECTIONS 7528 1140 CALL MOD2B (U,P,Q,LT,S,Z, NW,NY,NZ,NQ,ND,NN) 7529 1141 WRITE (JD32) LTER, ILK, NN, ((LT(J,J1), J=1,4), J1=1,NN) 7530 GO TO 1105 7531 C 7532 C END OF DATA FILE JD31 7533 1150 IF (NEWIT .LE. 0) STOP 77 7534 K = 0 7535 KDIF = NW 7536 1151 RETURN 7537 END 7538 FTN SUBM2A SUBROUTINE MOD2A (LORI,LT,RES,ROT,S,Z,U,P,Q,C,N,NW,NY,NZ,NQ,ND,NN) 7601 C 7602 C CORRECTION EQS AND CONTRIBUTION TO THE NORMAL EQS 7603 C 7604 COMMON WT(11), N1,N2,N3, INT1(2),NEWIT,INT2(5), 7605 1 NL,LAX(2,40),LTER(4),JL,ILK,K,INT3(2),INC,INT4(12),IL,IM,IT 7606 DIMENSION LORI(8,N),LT(4,NN), B(3,4),C(3,NZ),RES(3,NN), T(3),A(3) 7607 DOUBLE PRECISION U(3),P(3,ND), Q(7,7,NQ),S(7,NZ,NW),Z(NY,NZ) 7608 DOUBLE PRECISION ROT(3,4,NW), TT(3), WW 7609 EQUIVALENCE (WT7,WT(7)),(WT8,WT(8)),(ROOTW,WT(10)) 7610 C 7611 DO 1202 J1 = IL,IM 7612 1202 B(J1,J1) = 0. 7613 C 7614 C FORM THE NEEDED CORR EQS FOR EACH MODEL POINT 7615 MMP = 1 7616 DO 1249 MM = 1,NN 7617 IF (JL .EQ. 0) GO TO 1211 7618 IF (MM .EQ. LAX(2,JL)) GO TO 1249 7619 1211 J2 = LT(1,MM) 7620 IF (MM .EQ. 1) GO TO 1212 7621 C THE NEXT STATEMENT CAN BE TRUE FOR LAKE LEVEL POINTS ONLY 7622 IF (J2 .EQ. LT(1,MM-1)) GO TO 1212 7623 MMP = MMP + 1 7624 1212 M = J2 - K + 1 + INC 7625 M1 = (M-1) * NW - (M * (M-3)) / 2 7626 CALL MTRAN (LT,ROT,C,TT, NN,NW,NZ, MM) 7627 DO 1215 J = IL,IM 7628 C(J,NZ) = TT(J) + (LORI(J+1,J2) - LTER(J+1)) 7629 1215 RES(J,MM) = C(J,NZ) 7630 IF (NEWIT .LE. 0) GO TO 1249 7631 DO 1216 J = 1,3 7632 T(J) = TT(J) * 1D-5 7633 A(J) = 1. 7634 1216 B(J,4) = -T(J) 7635 B(1,2) = -T(3) 7636 B(1,3) = T(2) 7637 B(2,1) = T(3) 7638 B(2,3) = -T(1) 7639 B(3,1) = -T(2) 7640 B(3,2) = T(1) 7641 IF (IM .EQ. 2) GO TO 1240 7642 C 7643 C APPLY WEIGHTS TO THE CORR EQS (FOR PROJ CENTRES ONLY, IT = 1) 7644 IF (IT .EQ. 1) ROOTW = WT7 7645 DO 1228 J1 = IT,3 7646 IF (IT.EQ.1 .AND. J1.EQ.3) ROOTW = WT8 7647 IF (ROOTW .EQ. 1.) GO TO 1228 7648 A(J1) = ROOTW 7649 DO 1226 J2 = 1,4 7650 1226 B(J1,J2) = B(J1,J2) * ROOTW 7651 DO 1227 J2 = 1,NZ 7652 1227 C(J1,J2) = C(J1,J2) * ROOTW 7653 1228 CONTINUE 7654 C 7655 C FORM CONTRIBUTION OF THE CORR EQ TO THE NORMAL EQ, 7656 C 7657 C ON-DIAGONAL BAND 7658 1240 J2 = N2 * (MMP - 1) 7659 DO 1241 J = IL,IM 7660 WW = A(J)**2 7661 U(J) = U(J) + WW 7662 JP = J2 + J 7663 P(J,JP) = P(J,JP) - WW 7664 1241 Q(J,J,M1) = Q(J,J,M1) + WW 7665 DO 1243 J = 1,4 7666 JQ = N1 + J 7667 JP = JQ + J2 7668 DO 1242 L = IL,IM 7669 WW = A(L) * B(L,J) 7670 P(L,JP) = P(L,JP) + WW 7671 1242 Q(L,JQ,M1) = Q(L,JQ,M1) - WW 7672 DO 1243 L = 1,J 7673 DO 1243 L1 = IL,IM 7674 1243 Q(L+3,JQ,M1) = Q(L+3,JQ,M1) + B(L1,L) * B(L1,J) 7675 C VERTICAL BAND 7676 JP = ND - NZ 7677 DO 1248 J = 1,NZ 7678 JP = JP + 1 7679 DO 1245 L = IL,IM 7680 WW = A(L) * C(L,J) 7681 P(L,JP) = P(L,JP) + WW 7682 1245 S(L,J,M) = S(L,J,M) - WW 7683 DO 1246 L = 1,4 7684 DO 1246 L1 = IL,IM 7685 1246 S(L+3,J,M) = S(L+3,J,M) + B(L1,L) * C(L1,J) 7686 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1248 7687 DO 1247 L = 1,J 7688 IF (L .EQ. NZ) GO TO 1248 7689 DO 1247 L1 = IL,IM 7690 1247 Z(L,J) = Z(L,J) + C(L1,L) * C(L1,J) 7691 1248 CONTINUE 7692 1249 CONTINUE 7693 RETURN 7694 END 7695 FTN SUBM2B SUBROUTINE MOD2B (U,P,Q,LT,S,Z, NW,NY,NZ,NQ,ND,NN) 7801 C 7802 C ELIMINATE THE COORDINATE CORRECTIONS FOR ONE TERRAIN POINT 7803 C 7804 COMMON WT(11), N1,N2,N3, IN1(2),NEWIT,IN2(5), 7805 1 NL,LAX(2,40),IN3(5),ILK,K,IN4(2),INC,KKMAX,KAUXL,IN5(10),IL,IM 7806 DIMENSION LT(4,NN) 7807 DOUBLE PRECISION U(3),P(3,ND),Q(7,7,NQ),S(7,NZ,NW),Z(NY,NZ) 7808 DOUBLE PRECISION W(7,3),V(3) 7809 EQUIVALENCE (W(1,1),V(1)), (W1EN,WT(1)),(W1H,WT(2)), 7810 1 (W2EN,WT(3)),(W2H,WT(4)),(WSPH,WT(5)) 7811 NA = ND - NZ 7812 C 7813 C ADD CONTRIBUTION OF CONDITION EQS TO THE NORMAL EQS 7814 IF (ILK .GE. 7) GO TO 1266 7815 GO TO (1261,1262,1261, 1263,1264,1263), ILK 7816 1261 U(3) = U(3) + W1H 7817 IF (ILK .EQ. 1) GO TO 1268 7818 1262 U(1) = U(1) + W1EN 7819 GO TO 1265 7820 1263 IF (KAUXL .LT. 0) GO TO 1268 7821 U(3) = U(3) + W2H 7822 IF (ILK .EQ. 4) GO TO 1268 7823 1264 IF (KAUXL .LT. 0) GO TO 1268 7824 U(1) = U(1) + W2EN 7825 1265 U(2) = U(1) 7826 GO TO 1268 7827 1266 IF (ILK.EQ.7 .OR. ILK.GE.11) U(3) = U(3) + WSPH 7828 C 7829 C INVERT ON-DIAGONAL SUBMATRIX 7830 1268 DO 1269 J = IL,IM 7831 1269 U(J) = 1D0 / U(J) 7832 C 7833 C REDUCE EACH FOLLOWING ROW IN TURN, EXCEPT THE LAST ONE 7834 L2P = 1 7835 DO 1279 L2 = 1,NN 7836 L4 = LT(1,L2) 7837 IF (L2 .EQ. 1) GO TO 1271 7838 C THE FOLLOWING STATAMENT CAN BE TRUE ONLY FOR LAKE LEVEL POINTS 7839 IF (L4 .EQ. LT(1,L2-1)) GO TO 1279 7840 L2P = L2P + 1 7841 1271 M = L4 - K + 1 + INC 7842 M1 = (M-1) * NW - (M * (M-3)) / 2 7843 JJ = (L2P - 1) * N2 7844 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 7845 DO 1272 J = IL,IM 7846 DO 1272 J1 = 1,N2 7847 J2 = JJ + J1 7848 1272 W(J1,J) = P(J,J2) * U(J) 7849 C THEN, REDUCE SUBMATRICES IN DIAGONAL BAND 7850 L3P = L2P 7851 DO 1275 L3 = L2,NN 7852 IF (L3 .EQ. L2) GO TO 1273 7853 C THE FOLLOWING STATAMENT CAN BE TRUE ONLY FOR LAKE LEVEL POINTS 7854 IF (LT(1,L3) .EQ. LT(1,L3-1)) GO TO 1275 7855 L3P = L3P + 1 7856 1273 M2 = M1 + LT(1,L3) - L4 7857 J2 = N2 7858 DO 1274 J1 = 1,N2 7859 IF (L3 .EQ. L2) J2 = J1 7860 L = (L3P-1) * N2 + J1 7861 DO 1274 J = 1,J2 7862 DO 1274 J3 = IL,IM 7863 1274 Q(J,J1,M2) = Q(J,J1,M2) - W(J,J3) * P(J3,L) 7864 1275 CONTINUE 7865 C THEN, REDUCE SUBMATRIX AND VECTOR IN VERTICAL BAND 7866 DO 1277 L = 1,NZ 7867 J2 = NA + L 7868 DO 1277 J = 1,N2 7869 DO 1277 J3 = IL,IM 7870 1277 S(J,L,M) = S(J,L,M) - W(J,J3) * P(J3,J2) 7871 1279 CONTINUE 7872 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1290 7873 C 7874 C REDUCE THE LAST ROW 7875 DO 1283 L = 1,N3 7876 L2 = NA + L 7877 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 7878 DO 1281 J = IL,IM 7879 1281 V(J) = P(J,L2) * U(J) 7880 C THEN, REDUCE SUBMATRIX OF ADDED PARAMETERS 7881 C AND REDUCE ADJOINING VECTOR OF SECOND PARTS 7882 DO 1282 J = L,NZ 7883 J2 = NA + J 7884 DO 1282 J3 = IL,IM 7885 1282 Z(L,J) = Z(L,J) - V(J3) * P(J3,J2) 7886 1283 CONTINUE 7887 C 7888 1290 RETURN 7889 END 7890 FTN SUBM3 SUBROUTINE MODL3 (Q,S,Z, NW,NY,NZ,NQ) 8001 C 8002 C ELIMINATE THE ORIENTATION PARAMETERS, SOLVE FOR ADDED PARAMETERS 8003 C 8004 COMMON FLP(11), N1,N2,N3, IN1(2),NEWIT,LIST, JD31,JD32,JD33,IPR, 8005 1 NL,LAX(2,40),IN2(6),K,KK,KDIF 8006 DOUBLE PRECISION Q(7,7,NQ),S(7,NZ,NW),Z(NY,NZ),U(7,7),W(7,7),T(7) 8007 EQUIVALENCE (W(1,1),T(1)) 8008 C 8009 C ELIMINATE ONE OR MORE SETS OF N2 PARAMETERS 8010 KOUNT = KK 8011 DO 1340 LL = 1,KDIF 8012 L1 = (LL-1) * NW - ((LL-1) * (LL-2)) / 2 + 1 8013 C INVERT PIVOTAL SUBMATRIX 8014 DO 1301 J2 = 1,N2 8015 DO 1301 J1 = 1,J2 8016 1301 U(J1,J2) = Q(J1,J2,L1) 8017 CALL MOD3A (U,N2,N2,KOUNT,IPR) 8018 IF (NW .EQ. 1) GO TO 1320 8019 C 8020 C REDUCE EACH FOLLOWING ROW IN TURN, EXCEPT THE LAST ONE 8021 DO 1317 L2 = 2,NW 8022 MM = LL + L2 - 1 8023 M1 = (MM-1) * NW - ((MM-1) * (MM-2)) / 2 + 1 8024 IF (MM .GT. NW) GO TO 1320 8025 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 8026 K1 = L1 + L2 - 1 8027 DO 1311 J2 = 1,N2 8028 DO 1311 J1 = 1,N2 8029 W(J1,J2) = 0. 8030 DO 1311 M = 1,N2 8031 1311 W(J1,J2) = W(J1,J2) + Q(M,J1,K1) * U(M,J2) 8032 C THEN, REDUCE SUBMATRICES IN ON-DIAG BAND 8033 L4 = M1 + NW - MM 8034 DO 1312 L3 = M1,L4 8035 L5 = K1 + L3 - M1 8036 J3 = N2 8037 DO 1312 J2 = 1,N2 8038 IF (L3 .EQ. M1) J3 = J2 8039 DO 1312 J1 = 1,J3 8040 DO 1312 M = 1,N2 8041 1312 Q(J1,J2,L3) = Q(J1,J2,L3) - W(J1,M) * Q(M,J2,L5) 8042 C THEN, REDUCE SUBMATRIX AND VECTOR IN VERTICAL BAND 8043 DO 1316 L = 1,NZ 8044 DO 1316 J = 1,N2 8045 DO 1316 M = 1,N2 8046 1316 S(J,L,MM) = S(J,L,MM) - W(J,M) * S(M,L,LL) 8047 C FINALLY, STORE MULTIPLIER IN PIVOTAL ROW 8048 DO 1317 J = 1,N2 8049 DO 1317 L = 1,N2 8050 1317 Q(L,J,K1) = W(J,L) 8051 C 8052 C REDUCE THE LAST ROW 8053 1320 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1330 8054 DO 1323 L = 1,N3 8055 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 8056 DO 1321 J = 1,N2 8057 T(J) = 0. 8058 DO 1321 M = 1,N2 8059 1321 T(J) = T(J) + S(M,L,LL) * U(M,J) 8060 C THEN, REDUCE SUBMATRIX OF ADDED PARAMETERS 8061 C AND REDUCE ADJOINING VECTOR OF SECOND PARTS 8062 DO 1322 J = L,NZ 8063 DO 1322 M = 1,N2 8064 1322 Z(L,J) = Z(L,J) - T(M) * S(M,J,LL) 8065 C FINALLY, STORE MULTIPLIER IN PIVOTAL ROW 8066 DO 1323 J = 1,N2 8067 1323 S(J,L,LL) = T(J) 8068 C 8069 C REDUCE VECTOR OF SECOND PARTS IN PIVOTAL ROW 8070 1330 DO 1331 L = 1,N2 8071 T(L) = 0. 8072 DO 1331 J = 1,N2 8073 1331 T(L) = T(L) + U(L,J) * S(J,NZ,LL) 8074 DO 1332 J = 1,N2 8075 1332 S(J,NZ,LL) = T(J) 8076 KOUNT = KOUNT +1 8077 1340 CONTINUE 8078 IF (K .EQ. 0) GO TO 1390 8079 C 8080 C STORE AND SHIFT REDUCED EQUATIONS 8081 KD1 = KDIF + 1 8082 NW1 = NW + 1 8083 DO 1369 L = 1,NW1 8084 IF (L .EQ. 1) GO TO 1368 8085 IF (L .EQ. NW1) GO TO 1364 8086 L3 = NW1 + NW1 - L 8087 L1 = ((L-1) * L3) / 2 + 1 8088 L2 = L1 + NW - L 8089 L4 = 1 8090 IF (L .GT. KD1) L4 = ((L-KD1) * (L3 + KDIF)) / 2 + 1 8091 DO 1362 L3 = L1,L2 8092 DO 1361 J2 = 1,N2 8093 DO 1361 J1 = 1,N2 8094 1361 Q(J1,J2,L4) = Q(J1,J2,L3) 8095 1362 L4 = L4 + 1 8096 L5 = L4 8097 IF (L .GT. KD1) L5 = L4 + KDIF - 1 8098 GO TO 1365 8099 1364 L4 = L5 + 1 8100 L5 = NQ 8101 1365 DO 1366 L3 = L4,L5 8102 DO 1366 J2 = 1,N2 8103 DO 1366 J1 = 1,N2 8104 1366 Q(J1,J2,L3) = 0. 8105 IF (L .GT. KDIF) GO TO 1369 8106 1368 WRITE (JD33) (((Q(J1,J2,J), J1=1,N2), J2=1,N2), J=2,NW), 8107 1 ((S(J1,J2,L), J1=1,N2), J2=1,NZ) 8108 1369 CONTINUE 8109 DO 1379 L4 = 1,NW 8110 IF (L4 .GT. NW-KDIF) GO TO 1372 8111 L3 = L4 + KDIF 8112 DO 1371 J2 = 1,NZ 8113 DO 1371 J1 = 1,N2 8114 1371 S(J1,J2,L4) = S(J1,J2,L3) 8115 GO TO 1379 8116 1372 DO 1373 J2 = 1,NZ 8117 DO 1373 J1 = 1,N2 8118 1373 S(J1,J2,L4) = 0. 8119 1379 CONTINUE 8120 RETURN 8121 C 8122 1390 END FILE JD32 8123 REWIND JD32 8124 REWIND JD31 8125 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1391 8126 C SOLVE FOR THE ADDED PARAMETERS 8127 CALL MOD3A (Z,NY,NZ,0,IPR) 8128 1391 RETURN 8129 END 8130 FTN SUBM3A SUBROUTINE MOD3A (S,N,N1,KOUNT,IPR) 8201 C 8202 C INVERT, OR INVERT AND SOLVE IN N BY N+1 ARRAY (N = 7 OR N = NY) 8203 C 8204 DOUBLE PRECISION S(N,N1), REC, FAC, TEST 8205 31 FORMAT (36H BREAKDOWN OF SOLUTION AT MODEL # I4) 8206 32 FORMAT (37H ILL-CONDITIONED ADDED PARAMETER # I3, 16H EQUATED 8207 1TO ZERO) 8208 C 8209 IF (N .GT. 1) GO TO 1410 8210 IF (S(1,1) .LE. 0.) GO TO 1401 8211 S(1,1) = 1D0 / S(1,1) 8212 S(1,2) = S(1,2) * S(1,1) 8213 GO TO 1452 8214 1401 S(1,1) = 1D0 8215 S(1,2) = 0. 8216 WRITE (IPR,32) N 8217 GO TO 1452 8218 1410 DO 1411 K = 2,N 8219 DO 1411 J = K,N 8220 1411 S(J,K-1) = 0. 8221 TEST = 0. 8222 IF (N1 .EQ. N) GO TO 1420 8223 DO 1412 J = 1,N 8224 1412 TEST = TEST + S(J,J) 8225 TEST = 1D-8 * TEST / N 8226 C 8227 1420 DO 1449 J = 1,N 8228 IF (S(J,J) .GT. TEST) GO TO 1430 8229 IF (N1 .GT. N) GO TO 1421 8230 WRITE (IPR,31) KOUNT 8231 STOP 2 8232 1421 WRITE (IPR,32) J 8233 DO 1423 K = 1,J 8234 1423 S(K,J) = 0. 8235 DO 1424 K = 1,N1 8236 1424 S(J,K) = 0. 8237 S(J,J) = 1D0 8238 GO TO 1449 8239 C 8240 1430 J1 = J + 1 8241 REC = 1D0 / S(J,J) 8242 S(J,J) = 1D0 8243 DO 1439 K = 1,N1 8244 IF (K-J) 1436, 1439, 1431 8245 C ELIMINATION 8246 1431 FAC = S(J,K) * REC 8247 C TO PREVENT UNDERFLOW, 8248 IF (FAC.GT.1D-12 .OR. FAC.LT.-1D-12) GO TO 1432 8249 S(J,K) = 0. 8250 GO TO 1439 8251 1432 IF (K .GT. N) GO TO 1435 8252 DO 1433 L = 1,J 8253 1433 S(K,L) = S(K,L) - FAC * S(J,L) 8254 DO 1434 L = K,N1 8255 1434 S(K,L) = S(K,L) - FAC * S(J,L) 8256 1435 S(J,K) = FAC 8257 GO TO 1439 8258 C BACK SUBSTITUTION 8259 1436 FAC = S(K,J) * REC 8260 DO 1437 L = 1,K 8261 1437 S(K,L) = S(K,L) - FAC * S(J,L) 8262 IF (J1 .GT. N1) GO TO 1439 8263 DO 1438 L = J1,N1 8264 1438 S(K,L) = S(K,L) - FAC * S(J,L) 8265 1439 CONTINUE 8266 DO 1441 K = 1,J 8267 IF (K .EQ. J) GO TO 1442 8268 1441 S(J,K) = S(J,K) * REC 8269 1442 S(J,J) = REC 8270 1449 CONTINUE 8271 C 8272 DO 1451 K = 2,N 8273 J1 = K - 1 8274 DO 1451 J = 1,J1 8275 1451 S(J,K) = S(K,J) 8276 1452 RETURN 8277 END 8278 FTN SUBM4 SUBROUTINE MODL4 (Q,S,Z,LORI,V, N,NW,NY,NZ,NQ) 8401 C 8402 C SOLVE FOR THE ORIENTATION CORRECTIONS, CORRECT THE ORIENTATIONS 8403 C 8404 COMMON FLP(11), N1,N2,N3,ITER,ITR,NEWIT,LIST,JD31,JD32,JD33,IPR, 8405 1 NL,LAX(2,40),INT(13),LLAC,INT1,NIM 8406 DIMENSION LORI(8,N), V(7,N) 8407 DOUBLE PRECISION Q(7,7,NQ),S(7,NZ,NW),Z(NY,NZ), W(12), SD,SDEV 8408 41 FORMAT (1H1) 8409 42 FORMAT (24H0 RESULTS OF ITERATION I2 /) 8410 43 FORMAT (48H ADDED PAR, ST. DEV, CORRELATION MATRIX) 8411 44 FORMAT (1H 2D12.3, (1H T29, 12F7.3)) 8412 45 FORMAT (26H ORIENTATION PARAMETERS // (1H 8I10)) 8413 IF (LIST .GT. 0) WRITE (IPR,41) 8414 ITER = ITER + 1 8415 WRITE (IPR,42) ITER 8416 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1510 8417 C 8418 C FIRST, COMPUTE VARIANCE AND CORRELATION OF ADDED PARAMETERS 8419 WRITE (IPR,43) 8420 IF (N3 .EQ. 1) GO TO 1505 8421 DO 1502 J2 = 2,NY 8422 J = J2 - 1 8423 DO 1502 J1 = 1,J 8424 1502 Z(J1,J2) = Z(J1,J2) / (Z(J1,J1) * Z(J2,J2)) ** .5D0 8425 C ASSUME: STANDARD DEVIATION OF INPUT COORD. IS 5 UNITS OF L.S.D. 8426 1505 SD = 5.D0 8427 DO 1507 J1 = 1,NY 8428 SDEV = Z(J1,J1)**.5D0 * SD 8429 Z(J1,J1) = 1D0 8430 1507 WRITE (IPR,44) Z(J1,NZ), SDEV, (Z(J,J1), J=1,J1) 8431 WRITE (IPR,44) 8432 C 8433 C BACK SUBSTITUTION FOR ORIENTATION CORRECTIONS 8434 1510 DO 1529 JJ = 1,N 8435 MV = N + 1 - JJ 8436 IF (JJ .GT. NW) GO TO 1512 8437 M = NW + 1 - JJ 8438 IF (JJ .EQ. 1) GO TO 1525 8439 M1 = (M-1) * NW - (M * (M-3)) / 2 8440 L2 = JJ 8441 C NOTE: L2.LE.NW, M.GE.1 8442 GO TO 1520 8443 C SHIFT VECTORS AND READ A ROW OF STORED SUBMATRICES 8444 1512 DO 1513 J1 = 2,NW 8445 J = NW - J1 + 2 8446 DO 1513 L = 1,N2 8447 1513 S(L,NZ,J) = S(L,NZ,J-1) 8448 BACKSPACE JD33 8449 READ (JD33) (((Q(J1,J2,J), J1=1,N2), J2=1,N2), J=2,NW), 8450 1 ((S(J1,J2,1), J1=1,N2), J2=1,NZ) 8451 BACKSPACE JD33 8452 C CORRECTIONS FOR ONE PHOTOGRAPH 8453 1520 DO 1522 L3 = 2,L2 8454 L4 = NW - L2 + L3 8455 L7 = M1 + L3 - 1 8456 DO 1521 L = 1,N2 8457 DO 1521 J = 1,N2 8458 1521 S(L,NZ,M) = S(L,NZ,M) - Q(L,J,L7) * S(J,NZ,L4) 8459 1522 CONTINUE 8460 1525 IF (N3 .EQ. 0) GO TO 1527 8461 DO 1526 L = 1,N2 8462 DO 1526 J = 1,NY 8463 1526 S(L,NZ,M) = S(L,NZ,M) - S(L,J,M) * Z(J,NZ) 8464 1527 DO 1528 L = 1,N2 8465 1528 V(L,MV) = S(L,NZ,M) 8466 1529 CONTINUE 8467 C 8468 C CORRECT ALL ORIENTATION PARAMETERS 8469 DO 1554 J2 = 1,N 8470 DO 1551 J = 1,3 8471 COR = V(J,J2) + .5D0 8472 IF (COR .LT. .5) COR = COR - 1. 8473 KOR = COR 8474 1551 LORI(J+1,J2) = LORI(J+1,J2) + KOR 8475 DO 1552 J = 5,8 8476 1552 W(J) = LORI(J,J2) * 1D-7 8477 DO 1553 J = 4,N2 8478 1553 W(J-3) = V(J,J2) * .5E-5 8479 W(4) = 1D0 + W(4) 8480 W( 9) = W(4)*W(5) - W(3)*W(6) + W(2)*W(7) + W(1)*W(8) 8481 W(10) = W(3)*W(5) + W(4)*W(6) - W(1)*W(7) + W(2)*W(8) 8482 W(11) = -W(2)*W(5) + W(1)*W(6) + W(4)*W(7) + W(3)*W(8) 8483 W(12) = -W(1)*W(5) - W(2)*W(6) - W(3)*W(7) + W(4)*W(8) 8484 DO 1554 J = 9,12 8485 W(J) = W(J) * 1D7 + .5D0 8486 IF (W(J) .LT. .5D0) W(J) = W(J) - 1D0 8487 1554 LORI(J-4,J2) = W(J) 8488 WRITE (IPR,45) ((LORI(J1,J2), J1=1,8), J2=1,N) 8489 WRITE (JD31) ITER,N,NIM,NW,LLAC, ((LORI(J1,J2), J1=1,8), J2=1,N) 8490 RETURN 8491 END 8492 FTN SUBM5 SUBROUTINE MODL5 (LORI,ROT,LT,RES,Z,C, N,NIM,NW,NY,NZ) 8601 C 8602 C CORRECT THE APPROXIMATE COORDINATES OF THE TERRAIN POINTS 8603 C 8604 COMMON FLP(11), N1,N2,N3,ITER,ITR,NEWIT,LIST,JD31,JD32,JD33,IPR, 8605 1 NL,LAX(2,40),LTER(4),JL,ILK,K,KK,KDIF,INC,KKMAX 8606 DIMENSION LORI(8,N),LT(4,NIM),C(3,NY),RES(3,NIM),COR(3) 8607 DOUBLE PRECISION ROT(3,4,NW), Z(NY,NZ), ABCD(4),AR,BR,CR,DR 8608 DOUBLE PRECISION TSUM(3), TT(3), TT1, TER(3) 8609 EQUIVALENCE (AR,ABCD(1)),(BR,ABCD(2)),(CR,ABCD(3)),(DR,ABCD(4)) 8610 EQUIVALENCE (LTER1,LTER(1)),(TSUM(1),TER(1)) 8611 50 FORMAT (37H0 TERRAIN COORDINATES AND RESIDUALS) 8612 51 FORMAT (1H I2,I7, 2F12.2,F9.2, (1H T44, 4(F8.0,2F7.0))) 8613 52 FORMAT (1H I2,I7, 2F12.2,F9.2, (1H T44, 4(15X, F7.0))) 8614 53 FORMAT (1H 10X, 2F12.0, F9.0) 8615 54 FORMAT (1H 34X, F9.0) 8616 60 FORMAT (16H0 FILE UPDATED /) 8617 IF (LIST.LE.0 .AND. NEWIT.GT.1) GO TO 1602 8618 WRITE (IPR,50) 8619 NN = -1 8620 GO TO 1660 8621 1602 K = 0 8622 KK = 0 8623 INC = 0 8624 C 8625 C READ EACH RECORD 8626 1610 READ (JD32, END=1680) LTER, ILK, NN, ((LT(J,J1), J=1,4), J1=1,NN) 8627 IF (K .EQ. LT(1,1)) GO TO 1620 8628 C 8629 C SET THE ORIENTATION MATRICES 8630 C 8631 IF (LIST.GT.0 .OR. NEWIT.EQ.1) WRITE (IPR,51) 8632 K = LT(1,1) 8633 IF (K .GT. KKMAX) INC = K - KKMAX 8634 IF (KK .EQ. KKMAX) GO TO 1620 8635 IF (KK .GT. 0) GO TO 1611 8636 KK = K 8637 L1 = 1 8638 GO TO 1615 8639 C SHIFT THE AVAILABLE ORIENTATION MATRICES 8640 1611 KDIF = K - KK 8641 IF (K .GT. KKMAX) KDIF = KKMAX - KK 8642 KK = KK + KDIF 8643 M1 = KDIF + 1 8644 DO 1612 J2 = M1,NW 8645 M2 = J2 - KDIF 8646 DO 1612 J1 = 1,3 8647 DO 1612 J = 1,3 8648 1612 ROT(J,J1,M2) = ROT(J,J1,J2) 8649 L1 = NW - KDIF + 1 8650 C FORM THE MISSING MATRICES 8651 1615 DO 1618 L = L1,NW 8652 J = L + KK - 1 8653 DO 1616 J1 = 1,4 8654 1616 ABCD(J1) = LORI(J1+4,J) * 1D-7 8655 ROT(2,3,L) = 2D0*(BR*CR - AR*DR) 8656 ROT(3,1,L) = 2D0*(CR*AR - BR*DR) 8657 ROT(1,2,L) = 2D0*(AR*BR - CR*DR) 8658 ROT(3,2,L) = 2D0*(BR*CR + AR*DR) 8659 ROT(1,3,L) = 2D0*(CR*AR + BR*DR) 8660 ROT(2,1,L) = 2D0*(AR*BR + CR*DR) 8661 DO 1617 J1 = 1,4 8662 1617 ABCD(J1) = ABCD(J1)**2 8663 ROT(1,1,L) = DR + AR - BR - CR 8664 ROT(2,2,L) = DR - AR + BR - CR 8665 ROT(3,3,L) = DR - AR - BR + CR 8666 ROT(3,4,L) = (ROT(1,1,L)**2 + ROT(2,1,L)**2)**.5 8667 DO 1618 J1 = 1,2 8668 1618 ROT(J1,4,L) = ROT(J1,1,L) / ROT(3,4,L) 8669 C 8670 C TRANSFORM AND CORRECT POINT IN EACH MODEL 8671 1620 DO 1621 J = 1,3 8672 1621 TSUM(J) = 0. 8673 NSUM = NN 8674 IF (NL .EQ. 0) GO TO 1623 8675 DO 1622 JL = 1,NL 8676 IF (LTER(1) .NE. LAX(1,JL)) GO TO 1622 8677 IF (LAX(2,JL) .GT. 0) NSUM = NN - 1 8678 GO TO 1624 8679 1622 CONTINUE 8680 1623 JL = 0 8681 1624 IL = 1 8682 DO 1630 L = 1,NN 8683 IF (L.GT.1 .AND. ILK.EQ.10) IL = 3 8684 J2 = LT(1,L) 8685 CALL MTRAN (LT,ROT,C,TT, NN,NW,NY, L) 8686 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1628 8687 DO 1625 J = IL,3 8688 1625 COR(J) = 0. 8689 DO 1626 J1 = 1,N3 8690 W = Z(J1,NZ) 8691 DO 1626 J = IL,3 8692 1626 COR(J) = COR(J) + C(J,J1) * W 8693 DO 1627 J = IL,3 8694 1627 TT(J) = TT(J) - COR(J) 8695 1628 DO 1630 J = IL,3 8696 TT1 = TT(J) + LORI(J+1,J2) 8697 RES(J,L) = TT1 - LTER(J+1) 8698 IF (JL .EQ. 0) GO TO 1629 8699 IF (L .NE. LAX(2,JL)) GO TO 1629 8700 IF (.NOT. (J.LT.3 .AND. L.EQ.1 .AND. ILK.EQ.10)) GO TO 1630 8701 1629 TSUM(J) = TSUM(J) + TT1 8702 1630 CONTINUE 8703 C 8704 C CORRECT APPROXIMATE COORDINATES, COMPUTE RESIDUALS 8705 J1 = 1 8706 J2 = 3 8707 JG = ILK 8708 IF (JG .GT. 6) JG = 0 8709 IF (JG .GT. 3) JG = JG - 3 8710 IF (JG .LE. 0) GO TO 1634 8711 GO TO (1631, 1632, 1632), JG 8712 1631 J2 = 2 8713 GO TO 1634 8714 1632 J1 = J1 + JG 8715 1634 DO 1636 J = IL,3 8716 IF (NSUM .GT. 1) TSUM(J) = TSUM(J) / NSUM 8717 IF (LIST.LE.0 .AND. NEWIT.GT.1) GO TO 1636 8718 COR(J) = TSUM(J) - LTER(J+1) 8719 DO 1635 L = 1,NN 8720 1635 RES(J,L) = RES(J,L) - COR(J) 8721 1636 CONTINUE 8722 DO 1638 J = 1,3 8723 IF (J.LT.J1 .OR. J.GT.J2) GO TO 1638 8724 ROUND = .5 8725 IF (TSUM(J) .LT. 0.) ROUND = - ROUND 8726 JTER = TSUM(J) + ROUND 8727 LTER(J+1) = JTER 8728 1638 CONTINUE 8729 WRITE (JD31) LTER, ILK, NN, ((LT(J1,J2), J1=1,4), J2=1,NN) 8730 IF (LIST.LE.0 .AND. NEWIT.GT.1) GO TO 1610 8731 C 8732 C LIST COORDINATES AND RESIDUALS AT TERRAIN SCALE 8733 DO 1641 J = 1,3 8734 1641 TER(J) = LTER(J+1) * 1D-2 8735 IF (NN .GT. 1) GO TO 1642 8736 WRITE (IPR,51) ILK, LTER1, TER 8737 GO TO 1645 8738 1642 IF (ILK.EQ.10 .OR. ILK.EQ.12) GO TO 1643 8739 WRITE (IPR,51) ILK, LTER1, TER, ((RES(J1,J2),J1=1,3),J2=1,NN) 8740 GO TO 1645 8741 1643 WRITE (IPR,52) ILK, LTER1, TER, (RES(3,J2),J2=1,NN) 8742 GO TO 1610 8743 1645 IF (JG .LE. 0) GO TO 1650 8744 GO TO (1649, 1648, 1648), JG 8745 1648 WRITE (IPR,53) (COR(J), J=1,JG) 8746 GO TO 1650 8747 1649 WRITE (IPR,54) COR(3) 8748 C 8749 1650 IF (ILK .LE. 0) GO TO 1610 8750 IF (NL .EQ. 0) GO TO 1660 8751 DO 1651 J = 1,NL 8752 IF (LTER(1) .EQ. LAX(1,J)) GO TO 1610 8753 1651 CONTINUE 8754 C 8755 C COMPUTE ROOT MEAN SQUARE VALUES 8756 1660 CALL MOD5A (COR, RES, NIM, NN) 8757 IF (NN) 1602, 1690, 1610 8758 C 8759 1680 IF (LIST.LE.0 .AND. NEWIT.GT.1) GO TO 1690 8760 NN = 0 8761 GO TO 1660 8762 C 8763 1690 REWIND JD32 8764 END FILE JD31 8765 REWIND JD31 8766 WRITE (IPR,60) 8767 ITR = ITR + 1 8768 NEWIT = NEWIT - 1 8769 RETURN 8770 END 8771 FTN SUBM5A SUBROUTINE MOD5A (COR, RES, NIM, NN) 8801 C 8802 C ROOT MEAN SQUARE VALUES 8803 C 8804 COMMON FLP(11), INT1(10),IPR, 8805 1 NL,LAX(2,40),INT2(5),ILK,INT3(5),KAUXL,INT4(2),N 8806 DIMENSION COR(3), RES(3,NIM), M(9), E(4,9), A(4,7) 8807 DATA A /'CONT','ROL ','POIN','TS ','AUXI','L. P','OINT','S ', 8808 1 'CHEC','K PO','INTS',' ','SPEC','IAL ','HEIG','TS ', 8809 2 'CONT','ROL ','POIN','TS ','OTHE','R PO','INTS',' ', 8810 3 'PROJ','. CE','NTRE','S '/ 8811 55 FORMAT (1H0 T10,'ROOT MEAN SQUARE VALUES ME MN MPLA 8812 1N' T78,'MH' / 1H T10,'RESIDUALS') 8813 56 FORMAT (1H0 T6,4A4, T60,I10,F10.0, T82,'CENTIM.' T20,I10,3F10.0) 8814 57 FORMAT (1H0 T10,'DEVIATIONS FROM MEAN') 8815 58 FORMAT (1H0 T6,4A4, T20,I10,3F10.0, T70,F10.0, T82,'CENTIM.') 8816 IF (NN) 1600, 1680, 1660 8817 C 8818 1600 DO 1601 J = 1,9 8819 M(J) = 0 8820 DO 1601 J1 = 1,4 8821 1601 E(J1,J) = 0. 8822 GO TO 1695 8823 C CONTRIBUTIONS 8824 C 8825 1660 IF (ILK .LE. 7) GO TO 1661 8826 IF (ILK - 9) 1671, 1670, 1670 8827 C CONTRIBUTIONS OF RESIDUALS 8828 1661 JG = ILK 8829 J2 = 3 8830 GO TO (1668, 1664, 1663, 1667, 1664, 1662, 1667), ILK 8831 1662 J2 = 6 8832 1663 JG = JG - 1 8833 1664 M(JG) = M(JG) + 1 8834 DO 1665 J = 1,2 8835 1665 E(J,JG) = E(J,JG) + COR(J)**2 8836 GO TO (1668, 1672, 1666, 1667, 1671, 1666), ILK 8837 1666 JG = JG - 1 8838 GO TO 1668 8839 1667 J2 = 6 8840 1668 M(JG) = M(JG) + 1 8841 E(1,JG) = E(1,JG) + COR(3)**2 8842 GO TO 1672 8843 C CONTRIBUTIONS OF DEVIATIONS FROM MEAN 8844 1670 J2 = 9 8845 GO TO 1672 8846 1671 J2 = 6 8847 1672 IF (NN .EQ. 1) GO TO 1695 8848 M(J2) = M(J2) + NN 8849 DO 1673 L = 1,NN 8850 DO 1673 J = 1,3 8851 1673 E(J,J2) = E(J,J2) + RES(J,L)**2 8852 GO TO 1695 8853 C 8854 C RMS VALUES 8855 1680 DO 1688 J = 1,9 8856 IF (M(J) .EQ. 0) GO TO 1688 8857 J1 = 1 8858 GO TO (1684, 1682, 1681, 1684, 1682, 1681, 1684, 1682, 1681), J 8859 1681 E(4,J) = E(3,J) 8860 J1 = 4 8861 GO TO 1683 8862 1682 J1 = 3 8863 1683 E(3,J) = E(1,J) + E(2,J) 8864 1684 DO 1687 J2 = 1,J1 8865 1687 E(J2,J) = (E(J2,J) / M(J))**.5 8866 1688 CONTINUE 8867 WRITE (IPR,55) 8868 J2 = 0 8869 DO 1691 J = 1,4,3 8870 J2 = J2 + 1 8871 IF (KAUXL.LT.0 .AND. J.EQ.4) J2 = 3 8872 IF (M(J)+M(J+1) .EQ. 0) GO TO 1691 8873 WRITE (IPR,56) (A(J1,J2),J1=1,4), M(J),E(1,J), M(J+1),(E(J1,J+1), 8874 1 J1=1,3) 8875 1691 CONTINUE 8876 IF (M(7) .NE. 0) WRITE (IPR,56) (A(J1,4),J1=1,4), M(7),E(1,7) 8877 IF (N .EQ. 1) GO TO 1695 8878 WRITE (IPR,57) 8879 J2 = 4 8880 DO 1692 J = 3,9,3 8881 J2 = J2 + 1 8882 IF (M(J) .EQ. 0) GO TO 1692 8883 WRITE (IPR,58) (A(J1,J2),J1=1,4), M(J),(E(J1,J),J1=1,4) 8884 1692 CONTINUE 8885 1695 RETURN 8886 END 8887 FTN FTRAN SUBROUTINE MTRAN (LT,ROT,C,TT, NN,NW,NZ, MM) 9001 C 9002 C TRANSFORMATION, AND TERMS FOR ADDED PARAMETERS 9003 C 9004 COMMON WT(10),CURV, INT1(2),N3, INT2(2),NEWIT, INT3(6),LAX(2,40), 9005 1 INT4(6),K, INT5(2),INC 9006 DIMENSION LT(4,NN), C(3,NZ), T(3) 9007 DOUBLE PRECISION ROT(3,4,NW), TT(3) 9008 C 9009 M = LT(1,MM) - K + 1 + INC 9010 DO 1212 J = 1,3 9011 1212 TT(J) = ROT(J,1,M)*LT(2,MM) + ROT(J,2,M)*LT(3,MM) + 9012 1 ROT(J,3,M)*LT(4,MM) 9013 IF (NEWIT.GT.1 .OR. CURV.EQ.0.) GO TO 1214 9014 C APPLY EARTH CURVATURE OR/AND REFRACTION CORRECTIONS 9015 DO 1213 J = 1,2 9016 1213 TT(J) = TT(J) - TT(J) * TT(3) * CURV 9017 TT(3) = TT(3) + (TT(1)**2+TT(2)**2) * (.5 * CURV) 9018 1214 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1219 9019 C 9020 C TERMS FOR ADDED PARAMETERS IN ADJUSTMENT OF MODELS 9021 DO 1215 J = 1,3 9022 1215 T(J) = TT(J) * 1D-5 9023 DO 1217 J = 1,2 9024 DO 1217 J2 = 1,N3 9025 1217 C(J,J2) = 0. 9026 COS = ROT(1,4,M) 9027 SIN = ROT(2,4,M) 9028 U = T(1) * COS + T(2) * SIN 9029 V = T(2) * COS - T(1) * SIN 9030 Z = T(3) 9031 C 9032 C COEFFICIENTS OF ADDED PARAMETERS AS FUNCTIONS OF U AND V 9033 GO TO (11,12, 21,22, 31,32), N3 9034 C CORRECTIONS FOR LONGITUDINAL HEIGHT CURVATURE AND TORSION 9035 32 C(2,6) = -U * Z 9036 C(3,6) = U * V 9037 31 C(1,5) = -2. * U * Z 9038 C(3,5) = U * U 9039 C 9040 C CORRECTIONS FOR CHANGING SCALE AND AZIMUTH 9041 22 C(1,4) = U * V 9042 21 C(2,3) = U * V 9043 C 9044 C AFFINE CORRECTIONS 9045 C DIFFERENTIAL AZIMUTH CORRECTION 9046 12 C(2,2) = U 9047 C DIFFERENTIAL SCALE CORRECTION 9048 11 C(2,1) = V 9049 C 9050 C CONVERSION FROM CORRECTIONS DU AND DV TO CORRECTIONS DX AND DY 9051 DO 1218 J = 1,N3 9052 TEM = C(1,J) * COS - C(2,J) * SIN 9053 C(2,J) = C(1,J) * SIN + C(2,J) * COS 9054 1218 C(1,J) = TEM 9055 1219 RETURN 9056 END 9057