C ADJUSTMENT OF BUNDLES, JUNE 1978 G.H.S. 3001 C 3002 C ******* ADJUSTMENT ******* MAIN ROUTINE ******* 3003 C 3004 C *************************************************************** 3005 C * SET DIMENSIONS AND 8 CONSTANTS AS FOLLOWS, OR LARGER * 3006 C * MPH NUMBER OF PHOTOGRAPHS IN BLOCK * 3007 C * MIM MAXIMUM NUMBER OF IMAGES OF ONE TERRAIN POINT * 3008 C * MW BANDWIDTH * 3009 C * MY NUMBER OF ADDED PARAMETERS, BUT NOT SMALLER THAN 1 * 3010 C * MZ NUMBER OF ADDED PARAMETERS + 1 * 3011 C * MP 6 * MIM + MZ * 3012 C * MQ (MW * (MW + 1)) / 2 * 3013 C * MH MAXIMUM NUMBER OF SPECIAL HEIGHT CONTROL POINTS * 3014 C * LAIR(8,MPH),V(6,MPH),LT(3,MIM),RES(2,MIM),R(3,3,MW), * 3015 C * C(2,MZ),T(3,MP),LSPH(2,MH) * 3016 C * P(3,MP),Q(6,6,MQ),S(6,MZ,MW),Z(MY,MZ),ROT(3,3,MW) * 3017 C *************************************************************** 3018 C 3019 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,LIST, JD31,JD32,JD33,IPR,LF, 3020 1 NL,LAX(2,40),NOT(11),K,KK,KDIF,INC,KKMAX, 3021 2 N,NIM,NW,NZ,NY,NP,NQ, DOUBP 3022 DOUBLE PRECISION DOUBP(16) 3023 EQUIVALENCE (P(1,1),V(1,1)) 3024 C *************************************************************** 3025 C DIMENSIONS FOR 2 EXAMPLES BUNDLE C ************************* DIMENSION LAIR(8,29),V(6,29),LT(3,9),RES(2,9),R(3,3,9), 1 C(2,7),T(3,61),LSPH(2,1) DOUBLE PRECISION P(3,61),Q(6,6,45),S(6,7,9),Z(6,7),ROT(3,3,9) DATA MPH,MIM,MW, MY,MZ, MP,MQ, MH / 29,9,9, 6,7, 61,45, 1 / C *************************************************************** 3032 C 3033 N2 = 6 3034 N4 = 2 3035 N5 = 3 3036 ITR = 0 3037 100 CALL BUND1 (LAIR,MPH,LSPH,MH, Q,MQ,S,Z,C,MY,MZ,MW, MIM,MP) 3038 C 3039 C FORM THE NORMAL EQUATIONS, ELIMINATE THE TERRAIN POINT CORRECTIONS 3040 200 CALL BUND2 (LAIR, LSPH,MH, R,ROT, LT,RES, C,S,Z, P,T, Q) 3041 C 3042 C SOLVE THE REDUCED NORMAL EQUATIONS, CORRECT THE ORIENTATIONS 3043 CALL BUND3 (Q, S, Z) 3044 IF (K .GT. 0) GO TO 200 3045 CALL BUND4 (Q, S, Z, LAIR, V) 3046 C 3047 C CORRECT THE TERRAIN POINT COORDINATES 3048 CALL BUND5 (LAIR, V, ROT, LT, RES, Z, C, T) 3049 IF (NEWIT .EQ. 0) STOP 77 3050 GO TO 100 3051 END 3052 SUBROUTINE BUND1 (LAIR,MPH,LSPH,MH, Q,MQ,S,Z,C,MY,MZ,MW, MIM,MP) 3101 C 3102 C READ CONTROL CARDS AND RECORD WITH APPROX. ORIENTATION PARAMETERS 3103 C 3104 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,LIST, JD31,JD32,JD33,IPR,LF, 3105 1 NL,LAX(2,40),LTER(4),KECC,KREF, NH,INT(8), KKMAX, 3106 2 N,NIM,NW,NZ,NY,NP,NQ 3107 DIMENSION A(19), LAIR(8,MPH), LSPH(2,MH), C(N4,MZ) 3108 DOUBLE PRECISION Q(N2,N2,MQ), S(N2,MZ,MW), Z(MY,MZ) 3109 1 FORMAT (19A4, I4 / 6I4) 3110 4 FORMAT (I10, 20X, I10) 3111 8 FORMAT (10I8) 3112 10 FORMAT (1H1 19A4 //I6,14H ITERATION(S),I3,19H ADDED PARAMETER(S)) 3113 11 FORMAT (30H0 EARTH CURVATURE CORRECTION) 3114 12 FORMAT (25H0 REFRACTION CORRECTION) 3115 13 FORMAT ( 9H0 RELAX (1H T10, 10I8)) 3116 14 FORMAT (1H0 I8,' SPECIAL HEIGHTS') 3117 15 FORMAT (56H0 INITIAL TERRAIN COORDINATES AND PHOTOGRAPH RESIDUAL 3118 1S /) 3119 21 FORMAT (13H0 MORE THAN I3, 15H RELAXED POINTS) 3120 22 FORMAT (13H0 MORE THAN I4, 16H SPECIAL HEIGHTS) 3121 23 FORMAT (19H0 NO DATA IN FILE I3) 3122 24 FORMAT (55H0 THE VARIABLE DIMENSIONS MPH MIM MW MY MZ MP MQ 3123 1 / 17H ARE PRESENTLY T28, 7I4 3124 2 / 20H BUT MUST BE .GE. T28, 7I4) 3125 IF (ITR .GT. 0) GO TO 1025 3126 C DEFINE INPUT-OUTPUT DEVICES 3127 IRD = 5 3128 IPR = 6 3129 JD31 = 31 3130 JD32 = 32 3131 JD33 = 33 3132 C 3133 C READ THE CONTROL CARDS 3134 N1 = 3 3135 C NLMAX IS MAXIMUM NUMBER OF RELAXED POINTS, SEE ARRAY LAX 3136 NLMAX = 40 3137 READ (IRD,1) A, LIST, NEWIT, N3, KECC,KREF, NL, NH 3138 WRITE (IPR,10) A, NEWIT, N3 3139 IF (KECC .GE. 0) WRITE (IPR,11) 3140 IF (KREF .GE. 0) WRITE (IPR,12) 3141 C LIST OF POINTS TO BE RELAXED 3142 IF (NL .LE. 0) GO TO 1010 3143 IF (NL .GT. NLMAX) GO TO 1051 3144 READ (IRD,8) ((LAX(J1,J2), J1=1,2), J2=1,NL) 3145 DO 1002 J = 1,NL 3146 1002 LAX(1,J) = LAX(1,J) - 1000000 * (LAX(1,J)/1000000) 3147 WRITE (IPR,13) ((LAX(J1,J2), J1=1,2), J2=1,NL) 3148 C SPECIAL HEIGHTS 3149 1010 IF (NH .LE. 0) GO TO 1020 3150 DO 1011 J2 = 1,NH 3151 IF (J2 .GT. MH) GO TO 1052 3152 READ (IRD,4, END=1012) (LSPH(J,J2), J=1,2) 3153 IF (LSPH(1,J2) .EQ. 0) GO TO 1012 3154 1011 CONTINUE 3155 GO TO 1013 3156 1012 NH = J2 - 1 3157 1013 WRITE (IPR,14) NH 3158 C 3159 C READ RECORD WITH ORIENTATION PARAMETERS 3160 1020 READ (JD31,END=1053)ITER,N,NIM,NW,LF,((LAIR(J1,J2),J1=1,8),J2=1,N) 3161 NZ = N3 + 1 3162 NY = N3 3163 IF (NY .EQ. 0) NY = 1 3164 NP = N2 * NIM + NZ 3165 NQ = (NW * (NW + 1)) / 2 3166 IF (N.GT.MPH.OR.NIM.GT.MIM.OR.NW.GT.MW) GO TO 1054 3167 IF (NY.GT.MY.OR.NZ.GT.MZ.OR.NP.GT.MP.OR.NQ.GT.MQ) GO TO 1054 3168 KKMAX = N - NW + 1 3169 GO TO 1030 3170 1025 READ (JD31) 3171 C 3172 C ZERO THE NORMAL EQUATIONS 3173 1030 DO 1031 J = 1,NQ 3174 DO 1031 J2 = 1,N2 3175 DO 1031 J1 = 1,N2 3176 1031 Q(J1,J2,J) = 0. 3177 C NOTE: DO NOT REPLACE MZ AND MY BY NZ AND NY 3178 DO 1032 J = 1,NW 3179 DO 1032 J2 = 1,MZ 3180 DO 1032 J1 = 1,N2 3181 1032 S(J1,J2,J) = 0. 3182 IF (NZ .EQ. 1) GO TO 1040 3183 DO 1034 J2 = 1,NZ 3184 DO 1033 J = 1,2 3185 1033 C(J,J2) = 0. 3186 DO 1034 J1 = 1,MY 3187 1034 Z(J1,J2) = 0. 3188 1040 IF ((LIST.GT.0 .AND. ITR.EQ.0) .OR. NEWIT.LE.0) WRITE (IPR,15) 3189 DO 1041 J = 5,8 3190 1041 INT(J) = 0 3191 RETURN 3192 C 3193 1051 WRITE (IPR,21) NLMAX 3194 GO TO 1060 3195 1052 WRITE (IPR,22) MH 3196 GO TO 1060 3197 1053 WRITE (IPR,23) JD31 3198 GO TO 1060 3199 1054 WRITE (IPR,24) MPH,MIM,MW,MY,MZ,MP,MQ, N,NIM,NW,NY,NZ,NP,NQ 3200 1060 STOP 1 3201 END 3202 SUBROUTINE BUND2 (LAIR, LSPH,MH, R,ROT, LT,RES, C,S,Z, P,T, Q) 3301 C 3302 C FORM NORMAL EQUATIONS AND ELIMINATE TERRAIN POINT CORRECTIONS 3303 C 3304 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,LIST, JD31,JD32,JD33,IPR,LF, 3305 1 NL,LAX(2,40),LTER(4),NOT(2),NH,JL,NUM,JGR,NN,K,KK,KDIF,INC,KKMAX, 3306 2 N,NIM,NW,NZ,NY,NP,NQ, ABCD 3307 DIMENSION LAIR(8,N),LT(N5,NIM),R(3,3,NW),C(N4,NZ),T(3,NP), 3308 1 RES(N4,NIM),LSPH(2,NH) 3309 DOUBLE PRECISION P(3,NP),Q(N2,N2,NQ),S(N2,NZ,NW),Z(NY,NZ), 3310 1 ROT(3,3,NW), ABCD(16),AR,BR,CR,DR 3311 EQUIVALENCE (AR,ABCD(1)),(BR,ABCD(2)),(CR,ABCD(3)),(DR,ABCD(4)) 3312 IF (KK .GT. 0) GO TO 1110 3313 C 3314 C READ A RECORD, FORM THE NN ORIENTATION MATRICES 3315 1100 READ (JD31,END=1140) LTER, NN, ((LT(J,J1),J=1,N5),J1=1,NN) 3316 IF (K .EQ. LT(1,1)) GO TO 1120 3317 K = LT(1,1) 3318 IF (K .GT. KKMAX) INC = K - KKMAX 3319 IF (KK .EQ. KKMAX) GO TO 1120 3320 IF (KK .GT. 0) GO TO 1102 3321 KK = K 3322 L1 = 1 3323 GO TO 1115 3324 C 3325 C GO TO BUND3 TO ELIMINATE THE PARAMETERS OF ONE OR MORE PHOTOGRAPHS 3326 1102 KDIF = K - KK 3327 IF (K .GT. KKMAX) KDIF = KKMAX - KK 3328 IF (NEWIT) 1110, 1110, 1141 3329 C 3330 C UPON RETURN FROM BUND3, SHIFT THE AVAILABLE ROTATION MATRICES 3331 1110 KK = KK + KDIF 3332 M1 = KDIF + 1 3333 DO 1112 J2 = M1,NW 3334 M2 = J2 - KDIF 3335 DO 1112 J1 = 1,3 3336 DO 1112 J = 1,3 3337 1112 ROT(J,J1,M2) = ROT(J,J1,J2) 3338 L1 = NW - KDIF + 1 3339 C 3340 C FORM THE MISSING MATRICES 3341 1115 DO 1118 L = L1,NW 3342 J = L + KK - 1 3343 DO 1116 J1 = 1,4 3344 1116 ABCD(J1) = LAIR(J1+4,J) * 1D-7 3345 ROT(2,3,L) = 2D0*(BR*CR - AR*DR) 3346 ROT(3,1,L) = 2D0*(CR*AR - BR*DR) 3347 ROT(1,2,L) = 2D0*(AR*BR - CR*DR) 3348 ROT(3,2,L) = 2D0*(BR*CR + AR*DR) 3349 ROT(1,3,L) = 2D0*(CR*AR + BR*DR) 3350 ROT(2,1,L) = 2D0*(AR*BR + CR*DR) 3351 DO 1117 J1 = 1,4 3352 1117 ABCD(J1) = ABCD(J1)**2 3353 ROT(1,1,L) = DR + AR - BR - CR 3354 ROT(2,2,L) = DR - AR + BR - CR 3355 1118 ROT(3,3,L) = DR - AR - BR + CR 3356 DO 1119 J = 1,NW 3357 DO 1119 J2 = 1,3 3358 DO 1119 J1 = 1,3 3359 1119 R(J1,J2,J) = ROT(J1,J2,J) 3360 C 3361 C CHECK POINT AGAINST LIST OF SPECIAL HEIGHTS 3362 1120 JGR = LTER(1) / 1000000 3363 NUM = LTER(1) - 1000000 * JGR 3364 IF (ITR .GT. 0) GO TO 1124 3365 C ERASE ANY PREVIOUS SPECIAL HEIGHTS, INSERT NEW ONES 3366 IF (JGR .LT. 7) GO TO 1121 3367 LTER(1) = NUM 3368 JGR = 0 3369 1121 IF (NH .EQ. 0) GO TO 1124 3370 DO 1122 J = 1,NH 3371 IF (NUM .NE. LSPH(1,J)) GO TO 1122 3372 LTER(4) = LSPH(2,J) 3373 LTER(1) = NUM + 7000000 3374 JGR = 7 3375 GO TO 1124 3376 1122 CONTINUE 3377 C CHECK POINT AGAINST LIST OF RELAXED POINTS 3378 1124 IF (NL .EQ. 0) GO TO 1126 3379 DO 1125 JL = 1,NL 3380 IF (NUM .EQ. LAX(1,JL)) GO TO 1127 3381 1125 CONTINUE 3382 1126 JL = 0 3383 GO TO 1130 3384 1127 IF (LAX(2,JL).EQ.-1 .OR. NN.EQ.1) GO TO 1100 3385 IF (LAX(2,JL) .GT. 0) GO TO 1128 3386 LTER(1) = NUM 3387 JGR = 0 3388 GO TO 1130 3389 1128 IF (NN.EQ.2 .AND. JGR.NE.3) GO TO 1100 3390 C 3391 C FORM CORR EQ AND CONTRIBUTIONS TO NORMAL EQ, 3392 1130 ND = NN * N2 + NZ 3393 CALL FORM2 (LAIR, LT,RES,NN, R,ROT, S,Z, P,ND, Q, C) 3394 IF (NEWIT .LE. 0) GO TO 1100 3395 C 3396 C ELIMINATE THE COORDINATE CORRECTIONS 3397 CALL SOLV2 (P,T,ND, Q, LT,NN, S, Z) 3398 WRITE (JD32) LTER, NN, ((LT(J,J1), J=1,N5), J1=1,NN), 3399 1 ND, ((T(J,J1), J=1,N1), J1=1,ND) 3400 GO TO 1100 3401 C 3402 1140 IF (NEWIT .LE. 0) STOP 77 3403 K = 0 3404 KDIF = NW 3405 1141 RETURN 3406 END 3407 SUBROUTINE FORM2 (LAIR, LT,RES,NN, R,ROT, S,Z, P,ND, Q, C) 3501 C 3502 C CONTRIBUTION OF ONE TERRAIN POINT TO THE NORMAL EQS. 3503 C 3504 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,LIST, JD31,JD32,JD33,IPR,LF, 3505 1 NL,LAX(2,40),LTER(4),KECC,KREF,NIL,JL,NOT(3),K,KK,NT1,INC,NOR, 3506 2 N,NIM,NW,NZ,NY,NP,NQ, U,T,WW 3507 DIMENSION LAIR(8,N), LT(N5,NN) 3508 DIMENSION A(2,3), B(2,3), C(N4,NZ), R(3,3,NW), RES(N4,NN), TSP(3) 3509 DOUBLE PRECISION U(3,3), P(3,ND),Q(N2,N2,NQ),S(N2,NZ,NW),Z(NY,NZ) 3510 DOUBLE PRECISION ROT(3,3,NW), T(6), WW 3511 28 FORMAT (1H 3I10,I8, (1H T40, 5(F9.0,F7.0))) 3512 C 3513 F = LF 3514 ROOTW = .001 3515 JR = 0 3516 J2 = ND - N3 3517 DO 1202 J1 = 1,N1 3518 DO 1201 J = 1,N1 3519 1201 U(J1,J) = 0. 3520 DO 1202 J = J2,ND 3521 1202 P(J1,J) = 0. 3522 C 3523 C FORM TWO CORR EQ FOR EACH IMAGE POINT 3524 DO 1259 MM = 1,NN 3525 J2 = LT(1,MM) 3526 M = J2 - K + 1 + INC 3527 M1 = (M-1) * NW - (M * (M-3)) / 2 3528 DO 1211 J1 = 2,4 3529 T(J1+2) = LTER(J1) - LAIR(J1,J2) 3530 1211 TSP(J1-1) = T(J1+2) 3531 IF (NEWIT.GT.1 .OR. (KECC.LT.0 .AND. KREF.LT.0)) GO TO 1213 3532 DZ = 0. 3533 W = TSP(1)**2 + TSP(2)**2 3534 IF (KECC .GE. 0) DZ = W / 12756E5 3535 IF (KREF .LT. 0) GO TO 1212 3536 D = -TSP(3) * 1E-5 3537 H = LAIR(4,J2) * 1E-5 3538 DZ = DZ + (TSP(3) + W / TSP(3)) * D*(132.92+3.2460*D-H*(9.8091- 3539 1 .22737*H)-D*H*(.09693-.0030841*D+.0042192*H)) * 1E-7 3540 1212 T(6) = T(6) - DZ 3541 1213 DO 1214 J = 1,3 3542 1214 T(J) = ROT(1,J,M)*T(4) + ROT(2,J,M)*T(5) + ROT(3,J,M)*T(6) 3543 T(3) = -1D0 / T(3) 3544 JR = JR + 1 3545 DO 1215 J = 1,2 3546 C(J,NZ) = LT(J+1,MM) * 1D-6 - T(J) * T(3) 3547 1215 RES(J,JR) = -C(J,NZ) * F 3548 IF (NEWIT .LE. 0) GO TO 1259 3549 IF (JL .EQ. 0) GO TO 1216 3550 IF (MM .EQ. LAX(2,JL)) GO TO 1250 3551 1216 UL = LT(2,MM) * 1E-6 3552 VL = LT(3,MM) * 1E-6 3553 W = T(3) * 1D6 3554 DO 1217 J = 1,3 3555 A(1,J) = (R(J,1,M) + R(J,3,M) * UL) * W 3556 1217 A(2,J) = (R(J,2,M) + R(J,3,M) * VL) * W 3557 B(2,2) = UL * VL 3558 B(1,1) = -B(2,2) 3559 B(2,1) = -1. - VL**2 3560 B(1,2) = 1. + UL**2 3561 B(1,3) = VL 3562 B(2,3) = -UL 3563 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1218 3564 CALL ADPAR (UL, VL, C, NZ, N3) 3565 C TO WEIGHT THE CORRECTION EQUATIONS AS A FUNCTION OF DISTANCE FROM 3566 C IMAGE POINT TO PRINCIPAL POINT, REPLACE THE NEXT 3 STATEMENTS BY: 3567 C1218 ROOTW = 1. / (1. + 7. * (UL**2 + VL**2)) 3568 C IF (JL .EQ. 0) GO TO 1220 3569 C IF (LAX(2,JL) .EQ. 0) ROOTW = .001 3570 C1220 DO 1222 J1 = 1,2 3571 1218 IF (JL .EQ. 0) GO TO 1240 3572 IF (LAX(2,JL) .GT. 0) GO TO 1240 3573 DO 1222 J1 = 1,2 3574 DO 1221 J2 = 1,3 3575 A(J1,J2) = A(J1,J2) * ROOTW 3576 1221 B(J1,J2) = B(J1,J2) * ROOTW 3577 DO 1222 J2 = 1,NZ 3578 1222 C(J1,J2) = C(J1,J2) * ROOTW 3579 C 3580 C FORM CONTRIBUTION OF THE CORR EQ TO THE NORMAL EQ, 3581 C FOR ALL POINTS IN THE SAME WAY 3582 C 3583 C ON-DIAGONAL BAND 3584 1240 J2 = N2 * (MM-1) 3585 DO 1241 J = 1,3 3586 JP = J2 + J 3587 DO 1241 L = 1,J 3588 P(L,JP) = -A(1,L) * A(1,J) - A(2,L) * A(2,J) 3589 U(L,J) = U(L,J) - P(L,JP) 3590 Q(L,J,M1) = Q(L,J,M1) - P(L,JP) 3591 1241 Q(L+3,J+3,M1) = Q(L+3,J+3,M1) + (B(1,L)*B(1,J) + B(2,L)*B(2,J)) 3592 P(2,J2+1) = P(1,J2+2) 3593 P(3,J2+1) = P(1,J2+3) 3594 P(3,J2+2) = P(2,J2+3) 3595 DO 1242 J = 1,3 3596 JQ = N1 + J 3597 JP = JQ + J2 3598 DO 1242 L = 1,3 3599 P(L,JP) = A(1,L) * B(1,J) + A(2,L) * B(2,J) 3600 1242 Q(L,JQ,M1) = Q(L,JQ,M1) - P(L,JP) 3601 C VERTICAL BAND 3602 JP = NN * N2 3603 DO 1245 J = 1,NZ 3604 JP = JP + 1 3605 DO 1243 L = 1,3 3606 WW = A(1,L) * C(1,J) + A(2,L) * C(2,J) 3607 P(L,JP) = P(L,JP) + WW 3608 S(L,J,M) = S(L,J,M) - WW 3609 1243 S(L+3,J,M) = S(L+3,J,M) + (B(1,L) * C(1,J) + B(2,L) * C(2,J)) 3610 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1245 3611 DO 1244 L = 1,J 3612 IF (L .EQ. NZ) GO TO 1245 3613 1244 Z(L,J) = Z(L,J) + (C(1,L)*C(1,J) + C(2,L)*C(2,J)) 3614 1245 CONTINUE 3615 GO TO 1259 3616 C 3617 1250 J1 = N2 * (MM-1) + 1 3618 J2 = J1 + 5 3619 DO 1251 J = J1,J2 3620 DO 1251 L = 1,3 3621 1251 P(L,J) = 0. 3622 1259 CONTINUE 3623 C 3624 IF ((LIST.LE.0 .OR. ITR.GT.0) .AND. NEWIT.GT.0) GO TO 1261 3625 WRITE (IPR,28) LTER, RES 3626 1261 RETURN 3627 END 3628 SUBROUTINE SOLV2 (P,T,ND, Q, LT,NN, S, Z) 3701 C 3702 C ELIMINATE THE COORDINATE CORRECTIONS FOR ONE TERRAIN POINT 3703 C 3704 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,NIL(4), IPR,LF, 3705 1 NL,LAX(2,40),NOT(8),NUM,JGR,NON(4),INC,NOR(3), 3706 2 NW,NZ,NY,NP,NQ, U,REC,FAC,WTEN,WTH,WSPH 3707 DIMENSION LT(N5,NN),T(3,ND) 3708 DOUBLE PRECISION U(3,3), P(3,ND),Q(N2,N2,NQ),S(N2,NZ,NW),Z(NY,NZ) 3709 DOUBLE PRECISION REC,FAC, W(6,3),V(3), WTEN,WTH,WSPH 3710 EQUIVALENCE (W(1,1),V(1)) 3711 29 FORMAT (34H BREAKDOWN OF SOLUTION AT POINT I7) 3712 NA = NN * N2 3713 C 3714 C WEIGHT GROUND CONTROL 3715 WTEN = 25D0 3716 WTH = 25D0 3717 WSPH = WTH ** .25 3718 IF (JGR .EQ. 0) GO TO 1265 3719 GO TO (1261, 1262, 1261, 1265, 1265, 1265, 1263), JGR 3720 1261 U(3,3) = U(3,3) * WTH 3721 IF (JGR .EQ. 1) GO TO 1265 3722 1262 U(1,1) = U(1,1) * WTEN 3723 U(2,2) = U(2,2) * WTEN 3724 GO TO 1265 3725 1263 U(3,3) = U(3,3) * WSPH 3726 C 3727 C INVERT ON-DIAGONAL SUBMATRIX IN 3 BY 3 ARRAY 3728 1265 DO 411 J2 = 2,N1 3729 DO 411 J = J2,N1 3730 411 U(J,J2-1) = 0. 3731 DO 449 J = 1,N1 3732 IF (U(J,J) .GT. 0) GO TO 425 3733 WRITE (IPR,29) NUM 3734 STOP 2 3735 425 REC = 1D0 / U(J,J) 3736 U(J,J) = 1D0 3737 DO 439 J2 = 1,N1 3738 IF (J2 - J) 435, 439, 431 3739 C ELIMINATION 3740 431 FAC = U(J,J2) * REC 3741 DO 432 L = 1,N1 3742 IF (L.GT.J .AND. L.LT.J2) GO TO 432 3743 U(J2,L) = U(J2,L) - FAC * U(J,L) 3744 432 CONTINUE 3745 U(J,J2) = FAC 3746 GO TO 439 3747 C BACK SUBSTITUTION 3748 435 FAC = U(J2,J) * REC 3749 DO 436 L = 1,N1 3750 IF (L.GT.J2 .AND. L.LE.J) GO TO 436 3751 U(J2,L) = U(J2,L) - FAC * U(J,L) 3752 436 CONTINUE 3753 439 CONTINUE 3754 DO 441 J2 = 1,J 3755 IF (J2 .EQ. J) GO TO 442 3756 441 U(J,J2) = U(J,J2) * REC 3757 442 U(J,J) = REC 3758 449 CONTINUE 3759 DO 451 J2 = 2,N1 3760 J1 = J2 - 1 3761 DO 451 J = 1,J1 3762 451 U(J,J2) = U(J2,J) 3763 C 3764 C REDUCE EACH FOLLOWING ROW IN TURN, EXCEPT THE LAST ONE 3765 DO 1279 L2 = 1,NN 3766 M = LT(1,L2) - LT(1,1) + 1 + INC 3767 M1 = (M-1) * NW - (M * (M-3)) / 2 3768 JJ = (L2 - 1) * N2 3769 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 3770 DO 1271 J = 1,N1 3771 DO 1271 J1 = 1,N2 3772 J2 = JJ + J1 3773 1271 W(J1,J) = P(1,J2)*U(1,J) + P(2,J2)*U(2,J) + P(3,J2)*U(3,J) 3774 C THEN, REDUCE SUBMATRICES IN DIAGONAL BAND 3775 DO 1273 L3 = L2,NN 3776 M2 = M1 + LT(1,L3) - LT(1,L2) 3777 J2 = N2 3778 DO 1272 J1 = 1,N2 3779 IF (L3 .EQ. L2) J2 = J1 3780 L = (L3-1) * N2 + J1 3781 DO 1272 J = 1,J2 3782 1272 Q(J,J1,M2) = Q(J,J1,M2) -W(J,1)*P(1,L)-W(J,2)*P(2,L)-W(J,3)*P(3,L) 3783 1273 CONTINUE 3784 C THEN, REDUCE SUBMATRIX AND VECTOR IN VERTICAL BAND 3785 DO 1275 L = 1,NZ 3786 J2 = NA + L 3787 DO 1275 J = 1,N2 3788 1275 S(J,L,M) = S(J,L,M) -W(J,1)*P(1,J2)-W(J,2)*P(2,J2)-W(J,3)*P(3,J2) 3789 C FINALLY, STORE MULTIPLIER IN PIVOTAL ROW 3790 DO 1277 J = 1,N2 3791 J2 = JJ + J 3792 DO 1277 L = 1,N1 3793 1277 T(L,J2) = W(J,L) 3794 1279 CONTINUE 3795 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1290 3796 C 3797 C REDUCE THE LAST ROW 3798 DO 1283 L = 1,N3 3799 L2 = NA + L 3800 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 3801 DO 1281 J = 1,N1 3802 1281 V(J) = P(1,L2) * U(1,J) + P(2,L2) * U(2,J) + P(3,L2) * U(3,J) 3803 C THEN, REDUCE SUBMATRIX OF ADDED PARAMETERS 3804 C AND REDUCE ADJOINING VECTOR OF SECOND PARTS 3805 DO 1282 J = L,NZ 3806 J2 = NA + J 3807 1282 Z(L,J) = Z(L,J) - V(1)*P(1,J2) - V(2)*P(2,J2) - V(3)*P(3,J2) 3808 C STORE LAST MULTIPLIER IN PIVOTAL ROW 3809 DO 1283 J = 1,N1 3810 1283 T(J,L2) = V(J) 3811 C 3812 C REDUCE VECTOR OF SECOND PARTS IN PIVOTAL ROW 3813 1290 DO 1291 J = 1,N1 3814 1291 T(J,ND) = U(J,1) * P(1,ND) + U(J,2) * P(2,ND) + U(J,3) * P(3,ND) 3815 RETURN 3816 END 3817 SUBROUTINE BUND3 (Q, S, Z) 3901 C 3902 C ELIMINATE THE ORIENTATION PARAMETERS, SOLVE FOR ADDED PARAMETERS 3903 C 3904 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,LIST, JD31,JD32,JD33,IPR,LF, 3905 1 NL,LAX(2,40),NIL(11),K,KK,KDIF,NOT(4), NW,NZ,NY,NP,NQ 3906 DOUBLE PRECISION Q(N2,N2,NQ),S(N2,NZ,NW),Z(NY,NZ), 3907 1 U(6,6),W(6,6),T(6) 3908 EQUIVALENCE (W(1,1),T(1)) 3909 C 3910 C ELIMINATE THE PARAMETERS OF ONE OR MORE PHOTOGRAPHS 3911 KOUNT = KK 3912 DO 1340 LL = 1,KDIF 3913 L1 = (LL-1) * NW - ((LL-1) * (LL-2)) / 2 + 1 3914 C INVERT PIVOTAL SUBMATRIX 3915 DO 1301 J2 = 1,N2 3916 DO 1301 J1 = 1,J2 3917 1301 U(J1,J2) = Q(J1,J2,L1) 3918 CALL SOLV3 (U,N2,N2,KOUNT,IPR) 3919 IF (NW .EQ. 1) GO TO 1320 3920 C 3921 C REDUCE EACH FOLLOWING ROW IN TURN, EXCEPT THE LAST ONE 3922 DO 1317 L2 = 2,NW 3923 MM = LL + L2 - 1 3924 M1 = (MM-1) * NW - ((MM-1) * (MM-2)) / 2 + 1 3925 IF (MM .GT. NW) GO TO 1320 3926 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 3927 K1 = L1 + L2 - 1 3928 DO 1311 J2 = 1,N2 3929 DO 1311 J1 = 1,N2 3930 W(J1,J2) = 0. 3931 DO 1311 M = 1,N2 3932 1311 W(J1,J2) = W(J1,J2) + Q(M,J1,K1) * U(M,J2) 3933 C THEN, REDUCE SUBMATRICES IN ON-DIAG BAND 3934 L4 = M1 + NW - MM 3935 DO 1312 L3 = M1,L4 3936 L5 = K1 + L3 - M1 3937 J3 = N2 3938 DO 1312 J2 = 1,N2 3939 IF (L3 .EQ. M1) J3 = J2 3940 DO 1312 J1 = 1,J3 3941 DO 1312 M = 1,N2 3942 1312 Q(J1,J2,L3) = Q(J1,J2,L3) - W(J1,M) * Q(M,J2,L5) 3943 C THEN, REDUCE SUBMATRIX AND VECTOR IN VERTICAL BAND 3944 DO 1316 L = 1,NZ 3945 DO 1316 J = 1,N2 3946 DO 1316 M = 1,N2 3947 1316 S(J,L,MM) = S(J,L,MM) - W(J,M) * S(M,L,LL) 3948 C FINALLY, STORE MULTIPLIER IN PIVOTAL ROW 3949 DO 1317 J = 1,N2 3950 DO 1317 L = 1,N2 3951 1317 Q(L,J,K1) = W(J,L) 3952 C 3953 C REDUCE THE LAST ROW 3954 1320 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1330 3955 DO 1323 L = 1,N3 3956 C FIRST, FORM (TRANSPOSE OF OFF-DIAG MATRIX) * INVERSE 3957 DO 1321 J = 1,N2 3958 T(J) = 0. 3959 DO 1321 M = 1,N2 3960 1321 T(J) = T(J) + S(M,L,LL) * U(M,J) 3961 C THEN, REDUCE SUBMATRIX OF ADDED PARAMETERS 3962 C AND REDUCE ADJOINING VECTOR OF SECOND PARTS 3963 DO 1322 J = L,NZ 3964 DO 1322 M = 1,N2 3965 1322 Z(L,J) = Z(L,J) - T(M) * S(M,J,LL) 3966 C FINALLY, STORE MULTIPLIER IN PIVOTAL ROW 3967 DO 1323 J = 1,N2 3968 1323 S(J,L,LL) = T(J) 3969 C 3970 C REDUCE VECTOR OF SECOND PARTS IN PIVOTAL ROW 3971 1330 DO 1331 L = 1,N2 3972 T(L) = 0. 3973 DO 1331 J = 1,N2 3974 1331 T(L) = T(L) + U(L,J) * S(J,NZ,LL) 3975 DO 1332 J = 1,N2 3976 1332 S(J,NZ,LL) = T(J) 3977 KOUNT = KOUNT +1 3978 1340 CONTINUE 3979 IF (K .EQ. 0) GO TO 1390 3980 C 3981 C STORE AND SHIFT REDUCED EQUATIONS 3982 KD1 = KDIF + 1 3983 NW1 = NW + 1 3984 DO 1369 L = 1,NW1 3985 IF (L .EQ. 1) GO TO 1368 3986 IF (L .EQ. NW1) GO TO 1364 3987 L3 = NW1 + NW1 - L 3988 L1 = ((L-1) * L3) / 2 + 1 3989 L2 = L1 + NW - L 3990 L4 = 1 3991 IF (L .GT. KD1) L4 = ((L-KD1) * (L3 + KDIF)) / 2 + 1 3992 DO 1362 L3 = L1,L2 3993 DO 1361 J2 = 1,N2 3994 DO 1361 J1 = 1,N2 3995 1361 Q(J1,J2,L4) = Q(J1,J2,L3) 3996 1362 L4 = L4 + 1 3997 L5 = L4 3998 IF (L .GT. KD1) L5 = L4 + KDIF - 1 3999 GO TO 1365 4000 1364 L4 = L5 + 1 4001 L5 = NQ 4002 1365 DO 1366 L3 = L4,L5 4003 DO 1366 J2 = 1,N2 4004 DO 1366 J1 = 1,N2 4005 1366 Q(J1,J2,L3) = 0. 4006 IF (L .GT. KDIF) GO TO 1369 4007 1368 WRITE (JD33) (((Q(J1,J2,J), J1=1,N2), J2=1,N2), J=2,NW), 4008 1 ((S(J1,J2,L), J1=1,N2), J2=1,NZ) 4009 1369 CONTINUE 4010 DO 1379 L4 = 1,NW 4011 IF (L4 .GT. NW-KDIF) GO TO 1372 4012 L3 = L4 + KDIF 4013 DO 1371 J2 = 1,NZ 4014 DO 1371 J1 = 1,N2 4015 1371 S(J1,J2,L4) = S(J1,J2,L3) 4016 GO TO 1379 4017 1372 DO 1373 J2 = 1,NZ 4018 DO 1373 J1 = 1,N2 4019 1373 S(J1,J2,L4) = 0. 4020 1379 CONTINUE 4021 RETURN 4022 C 4023 1390 END FILE JD32 4024 REWIND JD32 4025 REWIND JD31 4026 IF (NEWIT.GT.1 .OR. N3.LE.1) GO TO 1391 4027 C SOLVE FOR THE ADDED PARAMETERS 4028 CALL SOLV3 (Z,NY,NZ,0,IPR) 4029 1391 RETURN 4030 END 4031 SUBROUTINE SOLV3 (S,N,N1,KOUNT,IPR) 4101 C 4102 C INVERT, OR INVERT AND SOLVE IN N BY N+1 ARRAY 4103 C 4104 DOUBLE PRECISION S(N,N1), REC, FAC, TEST 4105 31 FORMAT (41H BREAKDOWN OF SOLUTION AT PHOTOGRAPH # I5) 4106 32 FORMAT (37H ILL-CONDITIONED ADDED PARAMETER # I3, 16H EQUATED 4107 1TO ZERO) 4108 C 4109 DO 1411 K = 2,N 4110 DO 1411 J = K,N 4111 1411 S(J,K-1) = 0. 4112 TEST = 0. 4113 IF (N1 .EQ. N) GO TO 1420 4114 DO 1412 J = 1,N 4115 1412 TEST = TEST + S(J,J) 4116 TEST = 1D-8 * TEST / N 4117 C 4118 1420 DO 1449 J = 1,N 4119 IF (S(J,J) .GT. TEST) GO TO 1430 4120 IF (N1 .GT. N) GO TO 1421 4121 WRITE (IPR,31) KOUNT 4122 STOP 2 4123 1421 WRITE (IPR,32) J 4124 DO 1423 K = 1,J 4125 1423 S(K,J) = 0. 4126 DO 1424 K = 1,N1 4127 1424 S(J,K) = 0. 4128 S(J,J) = 1D0 4129 GO TO 1449 4130 C 4131 1430 J1 = J + 1 4132 REC = 1D0 / S(J,J) 4133 S(J,J) = 1D0 4134 DO 1439 K = 1,N1 4135 IF (K-J) 1435, 1439, 1431 4136 C ELIMINATION 4137 1431 FAC = S(J,K) * REC 4138 IF (K .GT. N) GO TO 1434 4139 DO 1432 L = 1,J 4140 1432 S(K,L) = S(K,L) - FAC * S(J,L) 4141 DO 1433 L = K,N1 4142 1433 S(K,L) = S(K,L) - FAC * S(J,L) 4143 1434 S(J,K) = FAC 4144 GO TO 1439 4145 C BACK SUBSTITUTION 4146 1435 FAC = S(K,J) * REC 4147 DO 1436 L = 1,K 4148 1436 S(K,L) = S(K,L) - FAC * S(J,L) 4149 IF (J1 .GT. N1) GO TO 1439 4150 DO 1437 L = J1,N1 4151 1437 S(K,L) = S(K,L) - FAC * S(J,L) 4152 1439 CONTINUE 4153 DO 1441 K = 1,J 4154 IF (K .EQ. J) GO TO 1442 4155 1441 S(J,K) = S(J,K) * REC 4156 1442 S(J,J) = REC 4157 1449 CONTINUE 4158 C 4159 DO 1451 K = 2,N 4160 J1 = K - 1 4161 DO 1451 J = 1,J1 4162 1451 S(J,K) = S(K,J) 4163 RETURN 4164 END 4165 SUBROUTINE BUND4 (Q, S, Z, LAIR, V) 4301 C 4302 C SOLVE FOR THE ORIENTATION CORRECTIONS 4303 C 4304 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,LIST, JD31,JD32,JD33,IPR,LF, 4305 1 NL,LAX(2,40),NOT(16), N,NIM,NW,NZ,NY,NP,NQ 4306 DIMENSION LAIR(8,N), V(N2,N) 4307 DOUBLE PRECISION Q(N2,N2,NQ),S(N2,NZ,NW),Z(NY,NZ), W(12), SD,SDEV 4308 41 FORMAT (1H1) 4309 42 FORMAT (24H0 RESULTS OF ITERATION I2 /) 4310 43 FORMAT (48H ADDED PAR, ST. DEV, CORRELATION MATRIX) 4311 44 FORMAT (1H 2D12.3, (1H T29, 12F7.3)) 4312 45 FORMAT (26H ORIENTATION PARAMETERS // (1H 8I10)) 4313 IF (LIST .GT. 0) WRITE (IPR,41) 4314 ITER = ITER + 1 4315 WRITE (IPR,42) ITER 4316 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1510 4317 C 4318 C FIRST, COMPUTE VARIANCE AND CORRELATION OF ADDED PARAMETERS 4319 WRITE (IPR,43) 4320 IF (N3 .GT. 1) GO TO 1501 4321 Z(1,1) = 1D0 / Z(1,1) 4322 Z(1,2) = Z(1,2) * Z(1,1) 4323 GO TO 1505 4324 1501 DO 1502 J2 = 2,NY 4325 J = J2 - 1 4326 DO 1502 J1 = 1,J 4327 1502 Z(J1,J2) = Z(J1,J2) / (Z(J1,J1) * Z(J2,J2)) ** .5D0 4328 C ASSUME: STANDARD DEVIATION OF MEASURED PHOT. COORD. = 2 MICROM. 4329 1505 SD = 2. / LF 4330 DO 1507 J1 = 1,NY 4331 SDEV = Z(J1,J1)**.5D0 * SD 4332 Z(J1,J1) = 1D0 4333 1507 WRITE (IPR,44) Z(J1,NZ), SDEV, (Z(J,J1), J=1,J1) 4334 WRITE (IPR,44) 4335 C 4336 C BACK SUBSTITUTION FOR ORIENTATION CORRECTIONS 4337 1510 DO 1529 JJ = 1,N 4338 MV = N + 1 - JJ 4339 IF (JJ .GT. NW) GO TO 1512 4340 M = NW + 1 - JJ 4341 IF (JJ .EQ. 1) GO TO 1525 4342 M1 = (M-1) * NW - (M * (M-3)) / 2 4343 L2 = JJ 4344 C NOTE: L2.LE.NW, M.GE.1 4345 GO TO 1520 4346 C SHIFT VECTORS AND READ A ROW OF STORED SUBMATRICES 4347 1512 DO 1513 J1 = 2,NW 4348 J = NW - J1 + 2 4349 DO 1513 L = 1,N2 4350 1513 S(L,NZ,J) = S(L,NZ,J-1) 4351 BACKSPACE JD33 4352 READ (JD33) (((Q(J1,J2,J), J1=1,N2), J2=1,N2), J=2,NW), 4353 1 ((S(J1,J2,1), J1=1,N2), J2=1,NZ) 4354 BACKSPACE JD33 4355 C CORRECTIONS FOR ONE PHOTOGRAPH 4356 1520 DO 1522 L3 = 2,L2 4357 L4 = NW - L2 + L3 4358 L7 = M1 + L3 - 1 4359 DO 1521 L = 1,N2 4360 DO 1521 J = 1,N2 4361 1521 S(L,NZ,M) = S(L,NZ,M) - Q(L,J,L7) * S(J,NZ,L4) 4362 1522 CONTINUE 4363 1525 IF (N3 .EQ. 0) GO TO 1527 4364 DO 1526 L = 1,N2 4365 DO 1526 J = 1,NY 4366 1526 S(L,NZ,M) = S(L,NZ,M) - S(L,J,M) * Z(J,NZ) 4367 1527 DO 1528 L = 1,N2 4368 1528 V(L,MV) = S(L,NZ,M) 4369 1529 CONTINUE 4370 C 4371 C CORRECT THE ORIENTATION PARAMETERS 4372 DO 1554 J2 = 1,N 4373 DO 1551 J = 1,3 4374 COR = V(J,J2) * 1D6 + .5D0 4375 IF (COR .LT. .5) COR = COR - 1. 4376 KOR = COR 4377 1551 LAIR(J+1,J2) = LAIR(J+1,J2) + KOR 4378 DO 1552 J = 1,4 4379 1552 W(J) = LAIR(J+4,J2) * 1D-7 4380 DO 1553 J = 4,6 4381 1553 W(J+1) = V(J,J2) * .5 4382 W( 9) = W(4)*W(5) - W(3)*W(6) + W(2)*W(7) + W(1) 4383 W(10) = W(3)*W(5) + W(4)*W(6) - W(1)*W(7) + W(2) 4384 W(11) = -W(2)*W(5) + W(1)*W(6) + W(4)*W(7) + W(3) 4385 W(12) = -W(1)*W(5) - W(2)*W(6) - W(3)*W(7) + W(4) 4386 W(8) = 1D7 / (W(9)**2 + W(10)**2 + W(11)**2 + W(12)**2)**.5 4387 DO 1554 J = 1,4 4388 W(J) = W(J+8) * W(8) + .5D0 4389 IF (W(J) .LT. .5D0) W(J) = W(J) - 1D0 4390 1554 LAIR(J+4,J2) = W(J) 4391 WRITE (JD31) ITER,N,NIM,NW, LF, ((LAIR(J1,J2), J1=1,8), J2=1,N) 4392 WRITE (IPR,45) ((LAIR(J1,J2), J1=1,8), J2=1,N) 4393 RETURN 4394 END 4395 SUBROUTINE BUND5 (LAIR, V, ROT, LT, RES, Z, C, P) 4501 C 4502 C CORRECT THE APPROXIMATE COORDINATES OF THE TERRAIN POINTS 4503 C 4504 COMMON N1,N2,N3,N4,N5, ITER,ITR,NEWIT,LIST, JD31,JD32,JD33,IPR,LF, 4505 1 NL,LAX(2,40),LTER(4),KECC,KREF,NOT(10), 4506 2 N,NIM,NW,NZ,NY,NP,NQ, ABCD,T 4507 DIMENSION LAIR(8,N),LT(N5,NIM),P(3,NP),C(N4,NY),RES(N4,NIM),COR(3) 4508 DIMENSION V(N2,N), KOR(3),JTER(4),TSP(3), M(7),E(3,7),A(4,5) 4509 DOUBLE PRECISION ROT(3,3,NW), T(6), Z(NY,NZ), ABCD(4),AR,BR,CR,DR 4510 EQUIVALENCE (AR,ABCD(1)),(BR,ABCD(2)),(CR,ABCD(3)),(DR,ABCD(4)) 4511 DATA A /'CONT','ROL ','POIN','TS ','CHEC','K PO','INTS',' ', 4512 1 'SPEC','IAL ','HEIG','TS ', 4513 2 'CONT','R PN','T IM','AGES','ALL ','OTHE','R IM','AGES'/ 4514 51 FORMAT (48H0 TERRAIN COORDINATES AND PHOTOGRAPH RESIDUALS) 4515 52 FORMAT (1H 3I10,I8, (1H T40, 5(F9.0,F7.0))) 4516 53 FORMAT (1H 11X, 2I10, I8) 4517 54 FORMAT (1H 31X, I8) 4518 55 FORMAT (1H0 T35,'ME,MU MN,MV MPLAN' T78,'MH') 4519 56 FORMAT (1H0 4A4, T60, I10, F10.1, T82,'CENTIM.' T20, I10, 3F10.1) 4520 57 FORMAT (1H0 4A4, T20, I10, 3F10.1, T82, 'MICROM.') 4521 IF (LIST.NE.0 .OR. NEWIT.EQ.1) WRITE (IPR,51) 4522 F = LF 4523 K = 0 4524 KK = 0 4525 KKMAX = N - NW + 1 4526 INC = 0 4527 JA = 1 4528 C 4529 C READ EACH RECORD 4530 1610 READ (JD32, END=1655) LTER, NN, ((LT(J,J1), J=1,N5), J1=1,NN), 4531 1 ND, ((P(J,J1), J=1,N1), J1=1,ND) 4532 NA = NN * N2 4533 IF (K .EQ. LT(1,1)) GO TO 1620 4534 IF (LIST.NE.0 .OR. NEWIT.EQ.1) WRITE (IPR,52) 4535 K = LT(1,1) 4536 IF (K .GT. KKMAX) INC = K - KKMAX 4537 IF (KK .EQ. KKMAX) GO TO 1620 4538 IF (KK .GT. 0) GO TO 1611 4539 KK = K 4540 L1 = 1 4541 GO TO 1615 4542 C 4543 C SHIFT THE AVAILABLE ROTATION MATRICES 4544 1611 KDIF = K - KK 4545 IF (K .GT. KKMAX) KDIF = KKMAX - KK 4546 KK = KK + KDIF 4547 M1 = KDIF + 1 4548 DO 1612 J2 = M1,NW 4549 M2 = J2 - KDIF 4550 DO 1612 J1 = 1,3 4551 DO 1612 J = 1,3 4552 1612 ROT(J,J1,M2) = ROT(J,J1,J2) 4553 L1 = NW - KDIF + 1 4554 C FORM THE MISSING MATRICES 4555 1615 DO 1618 L = L1,NW 4556 J = L + KK - 1 4557 DO 1616 J1 = 1,4 4558 1616 ABCD(J1) = LAIR(J1+4,J) * 1D-7 4559 ROT(2,3,L) = 2D0*(BR*CR - AR*DR) 4560 ROT(3,1,L) = 2D0*(CR*AR - BR*DR) 4561 ROT(1,2,L) = 2D0*(AR*BR - CR*DR) 4562 ROT(3,2,L) = 2D0*(BR*CR + AR*DR) 4563 ROT(1,3,L) = 2D0*(CR*AR + BR*DR) 4564 ROT(2,1,L) = 2D0*(AR*BR + CR*DR) 4565 DO 1617 J1 = 1,4 4566 1617 ABCD(J1) = ABCD(J1)**2 4567 ROT(1,1,L) = DR + AR - BR - CR 4568 ROT(2,2,L) = DR - AR + BR - CR 4569 1618 ROT(3,3,L) = DR - AR - BR + CR 4570 C 4571 C COMPUTE CORRECTIONS 4572 1620 DO 1621 J = 1,3 4573 JTER(J+1) = LTER(J+1) 4574 1621 COR(J) = P(J,ND) 4575 DO 1622 L1 = 1,NN 4576 L2 = LT(1,L1) 4577 L = (L1 - 1) * N2 4578 DO 1622 J = 1,N2 4579 J2 = L + J 4580 W = V(J,L2) 4581 DO 1622 J1 = 1,3 4582 1622 COR(J1) = COR(J1) - P(J1,J2) * W 4583 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1630 4584 DO 1623 J = 1,N3 4585 J2 = NA + J 4586 W = Z(J,NZ) 4587 DO 1623 J1 = 1,3 4588 1623 COR(J1) = COR(J1) - P(J1,J2) * W 4589 C 4590 C UPDATE UNKNOWN TERRAIN COORDINATES 4591 1630 JGR = LTER(1) / 1000000 4592 NUM = LTER(1) - 1000000 * JGR 4593 J1 = 1 4594 J2 = 3 4595 JG = JGR 4596 IF (JGR .EQ. 0) GO TO 1634 4597 IF (JG .GT. 3) JG = JG - 3 4598 GO TO (1631, 1632, 1633, 1631), JG 4599 1631 J2 = 2 4600 GO TO 1634 4601 1632 J1 = 3 4602 GO TO 1634 4603 1633 J1 = 4 4604 1634 DO 1636 J = 1,3 4605 ROUND = .5 4606 IF (COR(J) .LT. 0) ROUND = -ROUND 4607 KOR(J) = COR(J) * 1E6 + ROUND 4608 IF (J.LT.J1 .OR. J.GT.J2) GO TO 1636 4609 LTER(J+1) = LTER(J+1) + KOR(J) 4610 1636 CONTINUE 4611 WRITE (JD31) LTER, NN, ((LT(J1,J2), J1=1,N5), J2=1,NN) 4612 IF (LIST.EQ.0 .AND. NEWIT.GT.1) GO TO 1610 4613 C 4614 C COMPUTE RESIDUALS AT PHOTOGRAPH SCALE AND LIST 4615 DO 1647 L = 1,NN 4616 J2 = LT(1,L) 4617 DO 1641 J = 2,4 4618 T(J+2) = (JTER(J) - LAIR(J,J2)) + COR(J-1) * 1D6 4619 1641 TSP(J-1) = T(J+2) 4620 IF (NEWIT.GT.1 .OR. (KECC.LT.0 .AND. KREF.LT.0)) GO TO 1643 4621 DZ = 0. 4622 W = TSP(1)**2 + TSP(2)**2 4623 IF (KECC .GE. 0) DZ = W / 12756E5 4624 IF (KREF .LT. 0) GO TO 1642 4625 D = -TSP(3) * 1E-5 4626 H = LAIR(4,J2) * 1E-5 4627 DZ = DZ + (TSP(3) + W / TSP(3)) * D*(132.92+3.2460*D-H*(9.8091- 4628 1 .22737*H)-D*H*(.09693-.0030841*D+.0042192*H)) * 1E-7 4629 1642 T(6) = T(6) - DZ 4630 1643 L1 = J2 - K + 1 + INC 4631 DO 1644 J = 1,3 4632 1644 T(J) = ROT(1,J,L1)*T(4) + ROT(2,J,L1)*T(5) + ROT(3,J,L1)*T(6) 4633 C UNIT OF RES EQUALS L.S.D. OF PUNCHED FOCAL LENGTH 4634 DO 1645 J = 1,2 4635 W = -T(J) / T(3) - LT(J+1,L) * 1D-6 4636 1645 RES(J,L) = W * F 4637 IF (NEWIT.GT.1 .OR. N3.EQ.0) GO TO 1647 4638 XL = LT(2,L) * 1E-6 4639 YL = LT(3,L) * 1E-6 4640 CALL ADPAR (XL, YL, C, NZ, N3) 4641 DO 1646 J = 1,N3 4642 W = Z(J,NZ) 4643 DO 1646 J1 = 1,2 4644 1646 RES(J1,L) = RES(J1,L) + C(J1,J) * W * F 4645 1647 CONTINUE 4646 WRITE (IPR,52) LTER, ((RES(J1,J2),J1=1,N4),J2=1,NN) 4647 IF (JGR .EQ. 0) GO TO 1650 4648 GO TO (1649, 1648, 1648, 1649), JG 4649 1648 WRITE (IPR,53) (KOR(J), J=1,JG) 4650 GO TO 1650 4651 1649 WRITE (IPR,54) KOR(3) 4652 C 4653 1650 IF (NL .EQ. 0) GO TO 1660 4654 DO 1651 J = 1,NL 4655 IF (NUM .EQ. LAX(1,J)) GO TO 1610 4656 1651 CONTINUE 4657 GO TO 1660 4658 1655 IF (LIST.EQ.0 .AND. NEWIT.GT.1) GO TO 1690 4659 JA = 3 4660 C 4661 C COMPUTE ROOT MEAN SQUARE VALUES 4662 1660 IF (JA - 2) 1661, 1670, 1680 4663 1661 DO 1662 J = 1,7 4664 M(J) = 0 4665 DO 1662 J1 = 1,3 4666 1662 E(J1,J) = 0. 4667 JA = 2 4668 C CONTRIBUTIONS 4669 1670 IF (JGR .EQ. 0) GO TO 1679 4670 JG = JGR 4671 J2 = 3 4672 GO TO (1677, 1673, 1672, 1676, 1673, 1671, 1677), JGR 4673 1671 J2 = 6 4674 1672 JG = JG - 1 4675 1673 M(JG) = M(JG) + 1 4676 DO 1674 J = 1,2 4677 1674 E(J,JG) = E(J,JG) + COR(J)**2*1E12 4678 GO TO (1681, 1681, 1675, 1681, 1679, 1675), JGR 4679 1675 JG = JG - 1 4680 GO TO 1677 4681 1676 J2 = 6 4682 1677 M(JG) = M(JG) + 1 4683 E(1,JG) = E(1,JG) + COR(3)**2*1E12 4684 GO TO 1681 4685 1679 J2 = 6 4686 1681 M(J2) = M(J2) + NN 4687 DO 1682 L = 1,NN 4688 DO 1682 J = 1,2 4689 1682 E(J,J2) = E(J,J2) + RES(J,L)**2 4690 GO TO 1610 4691 C RMS ERRORS 4692 1680 DO 1688 J = 1,7 4693 IF (M(J) .EQ. 0) GO TO 1688 4694 J1 = 1 4695 GO TO (1686, 1685, 1685, 1686, 1685, 1685, 1686), J 4696 1685 J1 = 3 4697 E(3,J) = E(1,J) + E(2,J) 4698 1686 DO 1687 J2 = 1,J1 4699 1687 E(J2,J) = (E(J2,J) / M(J))**.5 4700 1688 CONTINUE 4701 WRITE (IPR,55) 4702 J2 = 0 4703 DO 1691 J = 1,4,3 4704 J2 = J2 + 1 4705 IF (M(J)+M(J+1) .EQ. 0) GO TO 1691 4706 WRITE (IPR,56) (A(J1,J2),J1=1,4), M(J),E(1,J), M(J+1),(E(J1,J+1), 4707 1 J1=1,3) 4708 1691 CONTINUE 4709 IF (M(7) .NE. 0) WRITE (IPR,56) (A(J1,3),J1=1,4), M(7),E(1,7) 4710 J2 = 3 4711 DO 1692 J = 3,6,3 4712 J2 = J2 + 1 4713 1692 WRITE (IPR,57) (A(J1,J2),J1=1,4), M(J),(E(J1,J),J1=1,3) 4714 C 4715 1690 REWIND JD32 4716 END FILE JD31 4717 REWIND JD31 4718 ITR = ITR + 1 4719 NEWIT = NEWIT - 1 4720 RETURN 4721 END 4722 SUBROUTINE ADPAR (U, V, C, NZ, N3) 4801 C 4802 C TERMS FOR ADDED PARAMETERS 4803 DIMENSION C(2,NZ) 4804 RR = U*U + V*V 4805 GO TO (11,12, 21,22, 31,32), N3 4806 C 4807 C TORSION CORRECTION 4808 32 UV = U * V 4809 C(1,6) = U * UV 4810 C(2,6) = V * UV 4811 C HEIGHT CURVATURE CORRECTION 4812 31 C(1,5) = U * RR 4813 C(2,5) = V * RR 4814 C 4815 C PSEUDO-CONFORMAL CORRECTIONS 4816 22 C(2,4) = RR 4817 21 C(1,3) = RR 4818 C 4819 C AFFINE CORRECTIONS 4820 C DIFFERENTIAL AZIMUTH CORRECTION 4821 12 C(2,2) = U 4822 C DIFFERENTIAL SCALE CORRECTION 4823 11 C(2,1) = V 4824 RETURN 4825 END 4826 C ADJUSTMENT OF BUNDLES, SORTING JUNE 1978 G.H.S. 1001 C 1002 C ***** SORT MONOCOMPARATOR DATA **** MAIN ROUTINE ***** 1003 C 1004 C ************************************************************** 1005 C * SET DIMENSIONS AND 7 CONSTANTS AS FOLLOWS, OR LARGER * 1006 C * MPH IS NUMBER OF PHOTOGRAPHS * 1007 C * MENH IS NUMBER OF CONTROL AND CHECK POINTS, PLUS ONE * 1008 C * MXY IS MAXIMUM NUMBER OF IMAGE POINTS IN ONE PHOTOGRAPH * 1009 C * MID IS NUMBER OF IMAGE POINTS * 1010 C * MIM IS MAXIMUM NUMBER OF IMAGES OF ONE TERRAIN POINT * 1011 C * MTY IS MAXIMUM NUMBER OF POINTS IN ARRAY LTIE (.LE. MXY) * 1012 C * MW IS BANDWIDTH * 1013 C * LAIR(8,MPH),LPH(2,MPH),LENH(4,MENH), * 1014 C * LXY(3,MXY),ID(MID),LTIE(3,MIM,MTY) * 1015 C * ROT(3,3,MW) * 1016 C ************************************************************** 1017 C 1018 COMMON INT(16), CLIST(161) 1019 EQUIVALENCE (INT(1),N), (INT(2),NIM), (INT(3),NW), (INT(5),NP), 1020 1 (LXY(1,1),ROT(1,1,1)), (ID(1),LENH(1,1)) 1021 C ************************************************************** 1022 C DIMENSIONS FOR ST FAITH TEST BLOCK MONO C ********************************** DIMENSION LAIR(8,14),LPH(2,14),LENH(4,225), 1 LXY(3, 75),ID( 700),LTIE(3,6,25) DOUBLE PRECISION ROT(3,3, 6) DATA MPH,MENH,MXY,MID,MIM,MTY,MW / 14,225, 75, 700,6,25,6 / C ************************************************************** 1029 C 1030 C READ PHOTOGRAPH SPECIFICATIONS, READ AIR STATIONS IN DESIRED ORDER 1031 C WITH APPROXIMATE VALUES OF ALL ORIENTATION PARAMETERS 1032 CALL COMM1 (LAIR,LPH,2,MPH) 1033 C 1034 C READ, REDUCE, AND SORT IMAGE COORDINATES 1035 CALL MONO2 (LXY,MXY, LPH,MPH, ID,MID) 1036 INT(3) = MW 1037 CALL MONO3 (LPH, ID, LXY,MXY, LTIE,MIM,MTY) 1038 C 1039 C READ GROUND CONTROL AND COMPUTE APPROXIMATE COORDINATES 1040 CALL COMM4 (LENH,MENH, LAIR, LTIE,MTY, ROT) 1041 STOP 77 1042 END 1043 SUBROUTINE MONO2 (LXY,MXY, LPH,MPH, ID,MID) 1201 C 1202 C READ COMPARATOR COORDINATES, COMPUTE REDUCED PHOTOGRAPH COORD. 1203 C CONSTRUCT ONE RECORD WITH REDUCED COORDINATES FOR EACH PHOTOGRAPH. 1204 C 1205 COMMON INT(6),JD31,JD32,JD33,IPR,JW11,JW12,JW13,LFOC,LSC(2), 1206 1 CLIST(160),DELR 1207 DOUBLE PRECISION S(2),FREC,T,W(2),ROUND 1208 DIMENSION L(3,3),LXY(3,MXY),LPH(2,MPH),ID(MID) 1209 EQUIVALENCE (INT(1),N), (INT(2),KODE),(INT(3),LENS),(INT(4),JTOT), 1210 1 (INT(5),NP) 1211 FREC = 1D0 / (1D-3 * LFOC) 1212 C 1213 8 FORMAT (10I8) 1214 31 FORMAT (15H0 PHOTOGRAPHS T27,'CODES' T38,'SCALE FACTORS' / 1215 1 T11, '#' T20, 'PH') 1216 32 FORMAT (1H 3I10, 3P2F10.5) 1217 33 FORMAT (1H0 I8, 12H PHOTOGRAPHS / I9, 13H IMAGE POINTS) 1218 40 FORMAT (27H0 NO MEASURED COORDINATES) 1219 41 FORMAT (23H0 MEASURED PHOTOGRAPH I8, 19H HAS NO AIR STATION) 1220 42 FORMAT (15H0 AIR STATION I8, 27H HAS NO MEASURED PHOTOGRAPH) 1221 43 FORMAT (16H0 INVALID CODE I3, 21H FOR ROTATION OF AXES) 1222 44 FORMAT (33H0 NO IMAGE POINTS IN PHOTOGRAPH I8) 1223 45 FORMAT (39H0 TOO MANY IMAGE POINTS IN PHOTOGRAPH I8) 1224 46 FORMAT (35H0 TOO MANY IMAGE POINTS IN BLOCK: I4) 1225 47 FORMAT (32H0 NO SPACE FOR NEXT PHOTOGRAPH) 1226 48 FORMAT (16H POINT NUMBER I8, 27H IS DUPLICATE IN PHOTOGRAPH I8) 1227 49 FORMAT ( 9H POINT I7, 14H IN PHOTOGRAPH I8, 1228 1 42H IS OUTSIDE RANGE OF LENS CORRECTION TABLE) 1229 C 1230 C FOR EACH PHOTOGRAPH, READ FIRST THE PRINCIPAL POINT 1231 C THEN THE MEASURED POINTS, 1 TO 3 PER CARD. 1232 C NO SEPARATION BETWEEN SETS OF CARDS FOR DIFFERENT PHOTOGRAPHS. 1233 C END OF DECK OR BLANK CARD ENDS READING 1234 WRITE (IPR,31) 1235 MISS = 0 1236 NPH = 1 1237 NP = 0 1238 JFOT = 0 1239 205 READ (JW12,8,END=240) JPHOT, ((L(J1,J2), J1=1,3), J2=1,3) 1240 IF (JPHOT .EQ. JFOT) GO TO 230 1241 IF (JFOT .NE. 0) GO TO 250 1242 C 1243 C PROCESS PRINCIPAL POINT CARD 1244 210 DO 211 JSEQ = 1,N 1245 IF (LPH(1,JSEQ) .NE. JPHOT) GO TO 211 1246 LPH(2,JSEQ) = NPH 1247 GO TO 212 1248 211 CONTINUE 1249 WRITE (IPR,41) JPHOT 1250 MISS = MISS + 1 1251 212 JFOT = JPHOT 1252 LPPX = L(2,1) 1253 LPPY = L(3,1) 1254 C PREPARE FOR CHANGE TO U,V COORDINATES 1255 IF (L(1,2) .NE. 0) KODE = L(1,2) 1256 L2 = 2 1257 L3 = 3 1258 JU = KODE / 10 1259 JV = KODE - 10 * JU 1260 IF (JU.GT.4 .OR. JV.GT.4) GO TO 293 1261 JU = JU + 1 1262 JV = JV + 1 1263 GO TO (221, 222, 221, 222, 221), JU 1264 221 GO TO (225, 293, 225, 293, 225), JV 1265 222 GO TO (293, 223, 293, 223, 293), JV 1266 223 L2 = 3 1267 L3 = 2 1268 C SET SCALE FACTORS 1269 225 DO 226 J = 1,2 1270 IF (L(J+1,2) .NE. 0) LSC(J) = L(J+1,2) 1271 226 S(J) = LSC(J) * 1D-8 1272 IF (JU.EQ.3 .OR. JU.EQ.4) S(1) = -S(1) 1273 IF (JV.EQ.2 .OR. JV.EQ.3) S(2) = -S(2) 1274 L1 = 0 1275 GO TO 205 1276 C 1277 C SQUEEZE IMAGE POINTS INTO ARRAY LXY, OMITTING DUPLICATES. 1278 C SHIFT ORIGIN AND, IF REQUIRED BY CODE, INTERCHANGE COORDINATES 1279 230 DO 233 J2 = 1,3 1280 IF (L(1,J2) .EQ. 0) GO TO 233 1281 IF (L1 .EQ. 0) GO TO 232 1282 DO 231 J = 1,L1 1283 IF (L(1,J2) .NE. LXY(1,J)) GO TO 231 1284 WRITE (IPR,48) L(1,J2), JFOT 1285 GO TO 233 1286 231 CONTINUE 1287 IF (L1 .EQ. MXY) GO TO 295 1288 232 L1 = L1 + 1 1289 LXY(1,L1) = L(1,J2) 1290 LXY(2,L1) = L(2,J2) - LPPX 1291 LXY(3,L1) = L(3,J2) - LPPY 1292 233 CONTINUE 1293 GO TO 205 1294 C 1295 240 IF (JFOT .EQ. 0) GO TO 290 1296 JPHOT = -1 1297 C 1298 C PROCESS THE POINTS IN ONE PHOTOGRAPH 1299 250 IF (L1 .EQ. 0) GO TO 294 1300 DO 259 J1 = 1,L1 1301 W(L2-1) = LXY(2,J1) * S(1) 1302 W(L3-1) = LXY(3,J1) * S(2) 1303 T = 0. 1304 IF (LENS .EQ. 0) GO TO 253 1305 T1 = W(1)**2 + W(2)**2 1306 T2 = T1**.5 + 1E-7 1307 J = T2 / DELR + 1. 1308 IF (J .LT. JTOT) GO TO 252 1309 WRITE (IPR,49) LXY(1,J1), JFOT 1310 GO TO 253 1311 252 T3 = DELR * J - T2 1312 T = (T3*CLIST(J) + (DELR-T3)*CLIST(J+1)) / T2 1313 253 DO 254 J = 1,2 1314 W(J) = (W(J) + W(J) * T) * FREC 1315 ROUND = .5D0 1316 IF (W(J) .LT. 0) ROUND = - ROUND 1317 254 LXY(J+1,J1) = W(J) * 1D6 + ROUND 1318 259 CONTINUE 1319 C 1320 C STORE THE DATA, ONE RECORD PER PHOTOGRAPH 1321 WRITE (JD33) L1, ((LXY(J1,J2), J1=1,3), J2=1,L1) 1322 WRITE (IPR,32) JSEQ, JFOT, KODE, S 1323 NP = NP + L1 1324 IF (JPHOT .LE. 0) GO TO 260 1325 NPH = NPH + 1 1326 IF (NPH .GT. MPH) GO TO 297 1327 GO TO 210 1328 C 1329 260 END FILE JD33 1330 REWIND JD33 1331 IF (NP .GT. MID) GO TO 296 1332 C 1333 C SORT PHOTOGRAPHS IN SEQUENCE FOR ADJUSTMENT 1334 NP = 0 1335 J = 1 1336 DO 279 J1 = 1,N 1337 J2 = LPH(2,J1) 1338 IF (J2 .NE. 0) GO TO 272 1339 WRITE (IPR,42) LPH(1,J1) 1340 MISS = MISS + 1 1341 GO TO 279 1342 272 IF (J2 - J) 273, 275, 274 1343 273 BACKSPACE JD33 1344 J = J - 1 1345 GO TO 272 1346 274 READ (JD33) 1347 J = J + 1 1348 GO TO 272 1349 275 READ (JD33) L1, ((LXY(L2,L3), L2 = 1,3), L3 = 1,L1) 1350 J = J + 1 1351 WRITE (JD32) J1, L1, ((LXY(L2,L3), L2 = 1,3), L3 = 1,L1) 1352 C STORE POINT NUMBERS IN ARRAY ID 1353 DO 278 J2 = 1,L1 1354 NP = NP + 1 1355 278 ID(NP) = LXY(1,J2) 1356 LPH(2,J1) = NP 1357 279 CONTINUE 1358 C 1359 IF (MISS .GT. 0) GO TO 299 1360 WRITE (IPR,33) N, NP 1361 END FILE JD32 1362 REWIND JD32 1363 REWIND JD33 1364 RETURN 1365 C 1366 C ERROR MESSAGES 1367 290 WRITE (IPR,40) 1368 GO TO 299 1369 293 WRITE (IPR,43) KODE 1370 GO TO 299 1371 294 WRITE (IPR,44) JFOT 1372 GO TO 299 1373 295 WRITE (IPR,45) JPHOT 1374 GO TO 299 1375 296 WRITE (IPR,46) NP 1376 GO TO 299 1377 297 WRITE (IPR,47) 1378 299 STOP 02 1379 END 1380 SUBROUTINE MONO3 (LPH, ID, LXY,MXY, LTIE,MIM,MTY) 1501 C 1502 C FOR EACH TERRAIN POINT, CODE INFORMATION ON EACH IMAGE POINT 1503 C AND PLACE IN ARRAY LTIE 1504 C 1505 COMMON N,NIM,NW,MTOT,NP,LIST, JD31,JD32,JD33,IPR 1506 DIMENSION LPH(2,N), ID(NP), LXY(3,MXY), LTIE(3,MIM,MTY),INT(6) 1507 EQUIVALENCE (INT(1),N) 1508 51 FORMAT (33H0 PHOTOGRAPHS WITH IMAGE POINTS /) 1509 52 FORMAT (1H I10, (1H T18, 12I8)) 1510 53 FORMAT (22H0 A MAXIMUM OF NIM = I3, 25H IMAGES PER TERRAIN POINT 1511 1 / 27H REQUIRED BANDWIDTH NW = I3) 1512 61 FORMAT (9H0 POINT I8, 23H, FIRST IN PHOTOGRAPH # I4, 1513 1 24H, OCCURS MORE THAN MIM = I3, 6H TIMES) 1514 62 FORMAT (9H0 POINT I8, 23H, FIRST IN PHOTOGRAPH # I4, 1515 1 26H, OVERFLOWS BANDWIDTH MW = I3) 1516 C 1517 MW = INT(3) 1518 DO 301 J = 2,4 1519 301 INT(J) = 0 1520 KK = 0 1521 JSUM = 0 1522 L3 = 1 1523 C 1524 C FOR EACH TERRAIN POINT IN EACH PHOTOGRAPH, 1525 DO 379 K = 1,N 1526 L1 = L3 1527 L2 = LPH(2,K) 1528 L3 = L2 + 1 1529 L4 = K + 2 * MW - 1 1530 IF (L4 .GT. N) L4 = N 1531 L4 = LPH(2,L4) 1532 MT = 0 1533 C 1534 C AT FIRST OCCURRENCE OF A TERRAIN POINT, STORE IN ARRAY LTIE 1535 C POINT NUMBER AND POINT SEQ.# 1536 DO 379 L = L1,L2 1537 IF (ID(L) .LE. 0) GO TO 329 1538 IF (MT .GT. 0) GO TO 311 1539 KMAX = K 1540 MZ = 1 1541 311 M1 = 1 1542 MT = MT + 1 1543 LTIE(1,1,MT) = ID(L) 1544 LTIE(2,1,MT) = L - (L1-1) 1545 IF (K .EQ. N) GO TO 326 1546 C 1547 C SEARCH FOR POINT NUMBER IN FOLLOWING PHOTOGRAPHS, 1548 C RECORD EACH OCCURRENCE BY PHOTOGRAPH SEQ.# AND POINT SEQ.# 1549 DO 325 J = L3,L4 1550 IF (ID(J) .NE. ID(L)) GO TO 325 1551 C FOUND AN IMAGE POINT IN ARRAY ID; NOW FIND ITS PHOTOGRAPH 1552 M1 = M1 + 1 1553 J1 = K + 1 1554 DO 321 J2 = J1,L4 1555 IF (LPH(2,J2) .GE. J) GO TO 322 1556 321 CONTINUE 1557 322 IF (J2 .GT. KMAX) KMAX = J2 1558 IF (M1.LE.MIM .AND. J2-K.LT.MW) GO TO 323 1559 IF (M1 .GT. MIM) WRITE (IPR,61) ID(L), K, MIM 1560 IF (J2-K.GE.MW) WRITE (IPR,62) ID(L), K, MW 1561 JSUM = JSUM + 1 1562 GO TO 324 1563 323 LTIE(1,M1,MT) = J2 1564 LTIE(2,M1,MT) = J - LPH(2,J2-1) 1565 324 ID(J) = -ID(J) 1566 325 CONTINUE 1567 C 1568 C COMPLETE THE PROCESSING OF A TERRAIN POINT 1569 326 ID(L) = -ID(L) 1570 IF (M1 .GT. MZ) MZ = M1 1571 C MARK UNUSED PART OF RECORD 1572 IF (M1 .LT. MIM) LTIE(1,M1+1,MT) = 7777 1573 IF (MT . EQ. MTY) GO TO 330 1574 329 IF (L.LT.L2 .OR. MT.EQ.0) GO TO 379 1575 C 1576 C INSERT REDUCED COORDINATES OF ALL TERRAIN POINTS IN ARRAY LTIE 1577 C 1578 C FIRST, UPDATE COUNTERS IN ARRAY INT 1579 C K AND KMAX ARE SEQUENCE NUMBERS OF FIRST AND LAST PHOTOGRAPH, 1580 C MT IS NUMBER OF TERRAIN POINTS COLLECTED IN LTIE 1581 330 IF (KMAX-K .GE. NW) NW = KMAX - K + 1 1582 IF (MZ .GT. NIM) NIM = MZ 1583 INT(4) = INT(4) + MT 1584 C READ ONLY NEEDED RECORDS INTO ARRAY LXY 1585 DO 349 J3 = K,KMAX 1586 DO 348 J2 = 1,MT 1587 IF (J3 .GT. K) GO TO 335 1588 J1 = 1 1589 GO TO 340 1590 335 DO 336 J1 = 2,MZ 1591 IF (LTIE(1,J1,J2) - J3) 336, 340, 348 1592 336 CONTINUE 1593 GO TO 348 1594 C READ NEEDED RECORD OF PHOTOGRAPH J3 INTO ARRAY LXY 1595 340 IF (J3 .EQ. KK) GO TO 345 1596 341 IF (KK+1 - J3) 342, 344, 343 1597 342 READ (JD32) 1598 KK = KK + 1 1599 GO TO 341 1600 343 BACKSPACE JD32 1601 KK = KK - 1 1602 GO TO 341 1603 344 READ (JD32) KK,JJ, ((LXY(J ,JJ1), J =1,3), JJ1=1,JJ) 1604 345 J = LTIE(2,J1,J2) 1605 DO 346 JJ = 2,3 1606 346 LTIE(JJ,J1,J2) = LXY(JJ,J) 1607 348 CONTINUE 1608 349 CONTINUE 1609 C 1610 C STORE ZEROS IN UNUSED LOCATIONS OF ARRAY LTIE 1611 DO 377 J2 = 1,MT 1612 DO 372 J1 = 2,MZ 1613 IF (LTIE(1,J1,J2) .NE. 7777) GO TO 372 1614 DO 371 J3 = J1,MZ 1615 DO 371 JJ = 1,3 1616 371 LTIE(JJ,J3,J2) = 0 1617 GO TO 377 1618 372 CONTINUE 1619 377 CONTINUE 1620 WRITE (JD33) K,MZ,MT, (((LTIE(J1,J2,J3),J1=1,3),J2=1,MZ),J3=1,MT) 1621 MT = 0 1622 379 CONTINUE 1623 C 1624 C OPTIONALLY, LIST ALL IMAGE POINTS IN EACH PHOTOGRAPH 1625 IF (LIST.LE.0 .AND. JSUM.EQ.0) GO TO 390 1626 DO 381 J = 1,NP 1627 381 ID(J) = -ID(J) 1628 WRITE (IPR,51) 1629 J2 = 0 1630 DO 382 L = 1,N 1631 J1 = J2 + 1 1632 J2 = LPH(2,L) 1633 WRITE (IPR,52) LPH(1,L), (ID(J), J=J1,J2) 1634 382 CONTINUE 1635 IF (JSUM .GT. 0) STOP 03 1636 C 1637 390 WRITE (IPR,53) NIM, NW 1638 END FILE JD33 1639 REWIND JD33 1640 REWIND JD32 1641 RETURN 1642 END 1643 C PUNCH ORIENTATION ELEMENTS AND TERRAIN COORDINATES C DIMENSION LAIR(8,300), LTER(4) 8 FORMAT (8I10) 9 FORMAT (19H0 NO DATA IN FILE I3) C C DEFINE INPUT AND OUTPUT DEVICES JD31 = 31 JD32 = 32 JD33 = 33 C C ORIENTATION ELEMENTS READ (JD31,END=208) ITER,N,NIM,NW,LF,((LAIR(J1,J2),J1=1,8),J2=1,N) DO 201 J2 = 1,N 201 WRITE (JD32,8) (LAIR(J1,J2), J1=1,8) C C TERRAIN COORDINATES 202 READ (JD31,END=209) (LTER(J), J=1,4) C OMIT TERRAIN POINTS WITH POINT NUMBER 7777 IF (LTER(1) .EQ. 7777) GO TO 202 LTER(1) = LTER(1) - (LTER(1) / 1000000) * 1000000 WRITE (JD33,8) (LTER(J), J=1,4) GO TO 202 C 208 WRITE (IPR,9) JD31 STOP 1 209 STOP 77 END C APPROXIMATE VALUES OF ORIENTATION PARAMETERS C C ALLOW A MAXIMUM OF 100 PHOTOGRAPHS PER FLIGHT LINE DOUBLE PRECISION B1, B2, COSI DIMENSION LAIR(8,100), LADD(4) MPH = 100 1 FORMAT (1H1) 3 FORMAT (8I10) 5 FORMAT (13H0 MORE THAN I4, 21H PHOTOGRAPHS IN STRIP) C C DEFINE INPUT AND OUTPUT DEVICES IRD = 5 IPR = 6 ICD = 7 WRITE (IPR,1) C C READ ALL AIR STATIONS, ONE FLIGHT LINE AT A TIME 110 DO 112 N = 1,MPH READ (IRD, 3, END=114) (LAIR(J,N), J=1,4) IF (LAIR(1,N) .LE. 0) GO TO 115 112 CONTINUE READ (IRD, 3, END=120) (LADD(J), J=1,4) IF (LADD(1) .GT. 0) GO TO 153 NXT = LADD(4) N = MPH GO TO 120 114 NXT = -1 GO TO 116 115 NXT = LAIR(4,N) 116 N = N - 1 IF (N-1) 142, 117, 120 117 LADD(2) = 1 LADD(3) = 0 GO TO 130 C C COMPUTE BASE 120 NM1 = N-1 DO 125 J = 2,4 125 LADD(J) = (LAIR(J,N) - LAIR(J,1)) / NM1 IF (N .EQ. 2) GO TO 130 IF (LAIR(4,2) .NE. 0) GO TO 130 C C INTERPOLATE E, N, H DO 126 J2 = 2,NM1 DO 126 J1 = 2,4 126 LAIR(J1,J2) = LAIR(J1,J2-1) + LADD(J1) C C COMPUTE A, B, C, D 130 GO TO (132,131,132,131), KODE 131 B1 = LADD(2) B2 = LADD(3) GO TO 133 132 B1 = LADD(3) B2 =-LADD(2) 133 IF (KODE .GT. 2) GO TO 134 B1 = -B1 B2 = -B2 134 COSI = B1 / (B1**2 + B2**2)**.5D0 LADD(3) = ((1D0 - COSI) * .5D0)**.5D0 * 1D7 LADD(4) = ((1D0 + COSI) * .5D0)**.5D0 * 1D7 IF (B2 .LT. 0) LADD(3) = -LADD(3) IF (LADD(4) .LT. 0) LADD(4) = -LADD(4) LADD(1) = 0 LADD(2) = 0 DO 135 J2 = 1,N DO 135 J1 = 1,4 135 LAIR(J1+4,J2) = LADD(J1) C C PUNCH ORIENTATION PARAMETERS DO 141 J2 = 1,N 141 WRITE (ICD,3) (LAIR(J1,J2), J1=1,8) 142 IF (NXT .LT. 0) GO TO 159 KODE = NXT IF (KODE .EQ. 0) KODE = 4 GO TO 110 C 153 WRITE (IPR,5) MPH 159 STOP END