C***********************************************************************00000010 C* *00000020 C* NRC34 AN INTRODUCTION TO ANALYTICAL STRIP TRIANGULATION *00000030 C* WITH A FORTRAN PROGRAM. G.H. SCHUT,PHOTOGRAMMETRIC RESERACH, *00000040 C* DIVISION OF APPLIED PHYSICS, NATIONAL RESEARCH COUNCIL OF *00000050 C* CANADA. FOR DESCRIPTION AND INSTRUCTIONS SEE REPORT NRC-9396 *00000060 C* DEC 1966 *00000070 C* *00000080 C***********************************************************************00000090 DOUBLE PRECISION DSQRT 00000100 DOUBLE PRECISION R(9),AR(12),AL(9),W(4),PP(4),T(9), 00000110 1 S(20),CLIST(162),SLIST(10),PHC(404) 00000120 DOUBLE PRECISION R1,R2,R3,R4,R5,R6,R7,R8,R9, AR1,AR2, 00000130 1 AR3,AR4,AR5,AR6,AR7,AR8,AR9, BX,BY,BZ, AL1,AL2,AL3, 00000140 2 AL4,AL5,AL6,AL7,AL8,AL9, U1,V1,U2,V2, PP1,PP2,PP3, 00000150 3 PP4, T1,T2,T3,T4,T5,T6,T7,T8,T9, A1,A2,A3,DBY,DBZ, 00000160 4 X1,Y1,Z1,X2,Y2,Z2,PX,PY,PZ, DELR,WANT, F,CX,CY,BX1, 00000170 5 CE,CR, FF,BY1 00000180 DIMENSION ID(101),LIST(10),LYST(10) 00000190 COMMON R1,R2,R3,R4,R5,R6,R7,R8,R9, AR1,AR2, 00000200 1 AR3,AR4,AR5,AR6,AR7,AR8,AR9, BX,BY,BZ, AL1,AL2,AL3, 00000210 2 AL4,AL5,AL6,AL7,AL8,AL9, U1,V1,U2,V2, PP1,PP2,PP3, 00000220 3 PP4, T1,T2,T3,T4,T5,T6,T7,T8,T9, S 00000230 EQUIVALENCE (R1,R(1)), (AR1,AR(1)), (AL1,AL(1)), 00000240 1 (U1,W(1)), (PP1,PP(1)), (T1,T(1)), (A1,S(16)), 00000250 2 (A2,S(17)), (A3,S(18)), (DBY,S(19)), (DBZ,S(20)) 00000260 1 FORMAT (I4,I5,F7.3,2F7.5,F7.0,F7.3,F7.4,I7,20X,I2) 00000270 2 FORMAT (I4, F5.1, 9F7.2, 6X, I2) 00000280 3 FORMAT (I4, I5, 4F7.3, I3) 00000290 4 FORMAT (1H,I4, 5X, 5F14.10) 00000300 5 FORMAT (I4, I5, 4I9) 00000310 55 FORMAT (1H, I4, I5, 4I9) 00000320 6 FORMAT (6H0ERROR,I2, 14H. EXIT AT,CARD2I6) 00000330 9 FORMAT (1H1) 00000340 9601 FORMAT(//13X,'CORRECTIONS'/13X,11('-')/) 9604 FORMAT(//12X,'P.C. COORDINATES'/12X,16('-')//6X,'PT.#',6X,'XS', *7X,'YS',7X,'ZS') 9606 FORMAT(/12X,'STRIP COORDINATES'/12X,17('-')//6X,'PT.#',6X,'XS', *7X,'YS',7X,'ZS',5X,'DISCREPANCIES(MIN. DISTANCE BETWEEN'/57X, *'CORRESPONDING RAYS)') 9610 FORMAT(2X,'ITR.#',6X,'A1',12X,'A2',12X,'A3',12X,'DBY',11X,'DBZ') 9603 FORMAT(//2X,'STRIP›',4X,'ELEMENTS OF ORIENTATION MATRIX'/2X, *'MODEL#',4X,30('-')) IRCD = 5 00000350 IWCD = 7 00000360 IPR = 6 00000370 C CALCULATE ACTUAL COMPUTATION TIME 00000380 94 FORMAT (7H0, TIME F5.1, 8H SECONDS) 00000390 C ITEM2=KLOCK(J) 00000400 C 00000410 1000 READ (IRCD,1) KK, KKK, F, CX, CY, BX1, CE, CR, NOCRD, K7 00000420 CR = CR / 1000. 00000430 CE =(CE / 12756. + CR) / (F*F) 00000440 WRITE (IPR,9) 00000450 KK = KK + 1 00000460 K8 = 1 00000470 IF (KK) 2901, 2901, 1001 00000480 1001 IF (KK-5) 1010, 1010, 2901 00000490 1010 K1 = 1 00000500 1011 K2 = K1 + 8 00000510 READ (IRCD,2) K3, T3, (CLIST(I),I=K1,K2), K5 00000520 IF (K1-1) 1012, 1012, 1014 00000530 1012 K6 = K3 00000540 DELR = T3 00000550 K8 = 2 00000560 IF (K6-162) 1014, 1014, 2902 00000570 1014 K7 = K7 + 1 00000580 K8 = 3 00000590 IF (K5-K7) 2902, 1015, 2902 00000600 1015 K1 = K1 + 9 00000610 IF (K6-K2) 1016, 1016, 1011 00000620 1016 DO 1017 I = 1,K6 00000630 1017 CLIST(I) = CLIST(I) / (DELR * 1000.) 00000640 K4MAX = 100 00000650 READ (IRCD,3) ID1, ID2, (PP(I),I=1,4), K6 00000660 1018 K4 = K6 00000670 MOD = 1 00000680 K3 = 1 00000690 PX = 200000. 00000700 PY = 400000. 00000710 PZ = 600000. 00000720 BY1 = 0. 00000730 1100 K = 1 00000740 KA = 1 00000750 K2 = 1 00000760 KK1 = 1 00000770 BY = BY1 00000780 BZ = 0. 00000790 WRITE(6,9601) WRITE(6,9610) DO 1101 I=1,10 00000810 1101 LIST(I) = -1 00000820 LI = 1 00000830 K5 = K4MAX 00000840 IF (K4 - K4MAX) 1102, 2030, 2030 00000850 1102 K5 = K4 00000860 K8 = 4 00000870 IF (K4-6) 2903, 2030, 2030 00000880 1110 K1 = 1 00000890 GO TO 2010 00000900 1111 DO 1112 I = 1,4 00000910 PHC(K2) = W(I) 00000920 1112 K2 = K2 + 1 00000930 IF (LI-10) 1113, 1113, 1120 00000940 1113 IF (K6 - 1) 1120, 1114, 1120 00000950 1114 LIST(LI) = ID(K1) 00000960 LI = LI + 1 00000970 KK1 = 2 00000980 1120 X2 = U2 00000990 Y2 = V2 00001000 Z2 = 1. 00001010 GO TO 2050 00001020 1121 K1 = K1 + 1 00001030 IF (K1-K5) 2010, 2010, 2100 00001040 1150 DO 1151 I=1,9 00001050 1151 AR(I) = R(I) 00001060 DO 1152 I=16,20 00001070 K8 = 30. * S(I) 00001080 IF (K8) 1200, 1152, 1200 00001090 1152 CONTINUE 00001100 GO TO 1300 00001110 1200 K = 2 00001120 1201 K2 = 1 00001130 GO TO 2030 00001140 1202 K1 = 0 00001150 1205 K1 = K1 + 1 00001160 DO 1203 I=1,4 00001170 W(I) = PHC(K2) 00001180 1203 K2 = K2 + 1 00001190 GO TO 2040 00001200 1204 IF (K1-K5) 1205, 1206, 1206 00001210 1206 GO TO (2999, 2100, 1302, 1602), K 00001220 1300 K = 3 00001230 KB = 1 00001240 GO TO 1201 00001250 1301 GO TO (1204, 1303), KB 00001260 1302 KB = 2 00001270 K1 = K4MAX + 1 00001280 K8 = K5 00001290 1303 K8 = K8 + 1 00001300 IF (K8-K4) 2010, 2010, 2100 00001310 1400 KA = 2 00001320 BY1 = BY 00001330 GO TO (1460, 1401), MOD 00001340 1401 KD = 1 00001350 BX = 1. 00001360 J = 0 00001370 GO TO (1410, 1430), KK2 00001380 1410 KKB = 3 00001390 GO TO (2999, 1411, 1412, 1414, 1413), KK 00001400 1411 KKB = 2 00001410 1412 KKA = 2 00001420 GO TO 1420 00001430 1413 KKB = 4 00001440 1414 KKA = 1 00001450 1420 M2 = KKA 00001460 1421 DO 1422 I = 1,4 00001470 K7 = 4 * M2 - 4 + I 00001480 1422 W(I) = PHC(K7) 00001490 GO TO 2040 00001500 1423 J = J + 1 00001510 SLIST(J) = SLIST(J) / T6 00001520 GO TO (1424, 1432), KK2 00001530 1424 M2 = M2 + 1 00001540 IF (M2-KKB) 1421, 1421, 1440 00001550 1430 M1 = 0 00001560 1435 M1 = M1 + 1 00001570 IF (LYST(M1)) 1440, 1431, 1431 00001580 1431 M2 = 0 00001590 1436 M2 = M2 + 1 00001600 IF (ID(M2)-LYST(M1)) 1432, 1421, 1432 00001610 1432 IF (M2-K5) 1436, 1437, 1437 00001620 1437 IF(M1-10)1435,1440,1440 00001630 1440 T1 = 0. 00001640 BX = 0. 00001650 DO 1443 I = 1,J 00001660 IF (SLIST(I)) 1443, 2999, 1442 00001670 1442 BX = BX + SLIST(I) 00001680 T1 = T1 + 1. 00001690 1443 CONTINUE 00001700 BX = BX / T1 00001710 T2 = 0. 00001720 DO 1455 I = 1,J 00001730 IF (SLIST(I)) 1455, 2999, 1451 00001740 1451 T3 = SLIST(I) - BX 00001750 IF (T3) 1452, 1455, 1453 00001760 1452 T3 = -T3 00001770 1453 IF (T3 - .9999999999 * T2) 1455, 1454, 1454 00001780 1454 K7 = I 00001790 T2 = T3 00001800 1455 CONTINUE 00001810 IF (T2 / BX - .0005) 1500, 1500, 1456 00001820 1456 SLIST(K7) = -SLIST(K7) 00001830 WRITE (IPR,55) K7 00001840 GO TO 1440 00001850 1460 BX = BX1 00001860 1500 BY = BX * BY 00001870 BZ = BX * BZ 00001880 GO TO (1510, 1501), MOD 00001890 1501 DO 1502 I = 1,9 00001900 1502 R(I) = AL(I) 00001910 J = 12 00001920 GO TO 2210 00001930 1510 WRITE(6,9603) DO 1511 I=1,3 1511 WRITE (IPR,4) ID1, AR(I), AR(I+3), AR(I+6) 00001950 WRITE (IPR,55) 00001960 WRITE(6,9604) K7 = 0 00001970 GO TO (1512, 1514), MOD 00001980 1512 IX = PX 00001990 IY = PY 00002000 IZ = PZ 00002010 WRITE (IPR,55) ID1, K7, IX, IY, IZ 00002020 IF (NOCRD) 1514, 1513, 1514 00002030 1513 WRITE (IWCD,5) ID1, K7, IX, IY, IZ 00002040 1514 IX = PX + BX 00002050 IY = PY + BY 00002060 IZ = PZ + BZ 00002070 WRITE (IPR,55) ID1, K7, IX, IY, IZ 00002080 WRITE(6,9606) IF (NOCRD) 1520, 1515, 1520 00002090 1515 WRITE (IWCD,5) ID1, K7, IX, IY, IZ 00002100 1520 KD = 2 00002110 GO TO (1521, 1527), KK1 00002120 1521 KKB = 5 00002130 GO TO (1600, 1526, 1523, 1524, 1525), KK 00002140 1523 KKB = 7 00002150 KKA = 6 00002160 GO TO 1527 00002170 1524 KKB = 6 00002180 KKA = 4 00002190 GO TO 1527 00002200 1525 KKB = 8 00002210 1526 KKA = 5 00002220 1527 KKC = 1 00002230 1600 K = 4 00002240 KB = 1 00002250 K2 = 1 00002260 GO TO 1202 00002270 1602 KB = 2 00002280 K1 = K4MAX + 1 00002290 1603 GO TO 2010 00002300 1610 IX = T4 + PX 00002310 IY = T5 + PY 00002320 IZ = T6 + PZ 00002330 K7 = WANT * DSQRT(T9) + .5 00002340 IF (WANT) 1611, 1612, 1612 00002350 1611 K7 = K7 - 1 00002360 1612 WRITE (IPR,55) ID1, ID(K1), IX, IY, IZ, K7 00002370 IF (NOCRD) 1614, 1613, 1614 00002380 1613 WRITE (IWCD,5) ID1, ID(K1), IX, IY, IZ, K7 00002390 1614 GO TO (1620, 1603), KB 00002400 1620 GO TO (1621, 1630), KK1 00002410 1621 GO TO (1204, 1622, 1622, 1622, 1622), KK 00002420 1622 IF (KKA-K1) 1623, 1623, 1204 00002430 1623 IF (K1-KKB) 1624, 1624, 1204 00002440 1624 SLIST(KKC) = AR7*(T4-BX) + AR8*(T5-BY) + AR9*(T6-BZ) 00002450 KKC = KKC + 1 00002460 GO TO 1204 00002470 1630 IF (KKC-10) 1631, 1631, 1204 00002480 1631 IF (LIST(KKC)-ID(K1)) 1204, 1624, 1204 00002490 1700 K8 = 5 00002500 GO TO (2904, 1701), KA 00002510 1701 ID1 = ID2 00002520 IF (ID1) 1000, 2999, 1702 00002530 1702 DO 1703 I=1,4 00002540 1703 PP(I) = W(I) 00002550 GO TO (1704, 1706), KK1 00002560 1704 GO TO (1018, 1706, 1705, 1706, 1705), KK 00002570 1705 K8 = 6 00002580 IF (KKB - K4) 1706, 1706, 2904 00002590 1706 K4 = K6 00002600 DO 1707 I = 1,10 00002610 1707 LYST(I) = LIST(I) 00002620 KK2 = KK1 00002630 PX = PX + BX 00002640 PY = PY + BY 00002650 PZ = PZ + BZ 00002660 DO 1708 I = 1,9 00002670 1708 AL(I) = AR(I) 00002680 MOD = 2 00002690 GO TO 1100 00002700 2010 READ (IRCD,3) ID2, ID(K1), (W(I),I=1,4), K6 00002710 IF (ID2-ID1) 1700, 2011, 1700 00002720 2011 U1 = (U1 - PP1) * CX 00002730 V1 = (V1 - PP2) * CY 00002740 U2 = (U2 - PP3) * CX 00002750 V2 = (V2 - PP4) * CY 00002760 KC = 1 00002770 T1 = U1 * U1 + V1 * V1 00002780 GO TO (2012, 2021), K3 00002790 2012 K3 = 2 00002800 FF = F 00002810 IF (U2-U1) 2013, 2013, 2021 00002820 2013 FF = -F 00002830 2021 T2 = DSQRT(T1) + .0000001 00002860 J = T2 / DELR + 1. 00002870 T3 = J 00002880 T3 = DELR * T3 - T2 00002890 T4 = (T3*CLIST(J)+(DELR-T3)*CLIST(J+1))/T2 +CR +CE*T1 00002900 K7 = KC + 1 00002910 DO 2022 I = KC,K7 00002920 2022 W(I) = (W(I) + W(I) * T4) / FF 00002930 GO TO (2024, 2999, 2025), KC 00002940 2024 KC = 3 00002950 T1 = U2 * U2 + V2 * V2 00002960 50 FORMAT('0 DLT ERROR_NODE',I4,'NOT IN LIST',I4) 00002970 GO TO 2021 00002980 2025 GO TO (1111, 2999, 2040, 2040), K 00002990 2030 DO 2031 M1 = 1,20 00003000 2031 S(M1) = 0. 00003010 GO TO (1110, 1202, 1202), K 00003020 2040 X2 = AR1 * U2 + AR4 * V2 + AR7 00003030 Y2 = AR2 * U2 + AR5 * V2 + AR8 00003040 Z2 = AR3 * U2 + AR6 * V2 + AR9 00003050 GO TO (2050, 2230), KA 00003060 2050 T7 = BY - V1 * BZ 00003070 T8 = U1 * BZ - 1. 00003080 T9 = V1 - U1 * BY 00003090 T1 = Y2 * T9 - Z2 * T8 00003100 T2 = Z2 * T7 - X2 * T9 00003110 T3 = X2 * T8 - Y2 * T7 00003120 T4 = X2 - U1 * Z2 00003130 T5 = U1 * Y2 - V1 * X2 00003140 T6 = Y2 - V1 * Z2 - T4 * BY - T5 * BZ 00003150 GO TO (2090, 2090, 2051), K 00003160 2051 IF (KKK) 2090, 2090, 2052 00003170 2052 T7 = 1./((.14+U1*U1+V1*V1)**2+(.14+U2*U2+V2*V2)**2) 00003180 T7 = DSQRT(T7) 00003190 DO 2053 I = 1,6 00003200 2053 T(I) = T7 * T(I) 00003210 2090 M3 = 1 00003220 DO 2091 M1 = 1,5 00003230 DO 2091 M2 = M1,6 00003240 S(M3) = S(M3) + T(M1) * T(M2) 00003250 2091 M3 = M3+1 00003260 GO TO (1121, 1204, 1301), K 00003270 2100 K7 = 0 00003280 DO 101 L1 = 2,5 00003290 M1 = K7 + 1 00003300 K7 = M1 + 7 - L1 00003310 M3 = K7 00003320 M4 = M1 + 1 00003330 M5 = K7 - 1 00003340 DO 102 L2 = M4,M5 00003350 T1 = S(L2) / S(M1) 00003360 DO 103 M2 = L2,K7 00003370 M3 = M3 + 1 00003380 103 S(M3) = S(M3) - T1 * S(M2) 00003390 102 S(L2) = T1 00003400 101 S(K7) = S(K7) / S(M1) 00003410 DBZ = DBZ / S(19) 00003420 M1 = 20 00003430 DO 105 L1 = 2,5 00003440 M1 = M1 - 1 00003450 M2 = 20 00003460 M3 = M3 - 2 00003470 S(M1) = S(M3) 00003480 DO 105 L2 = 2,L1 00003490 M3 = M3 - 1 00003500 S(M1) = S(M1) - S(M2) * S(M3) 00003510 105 M2 = M2 - 1 00003520 WRITE (IPR,4) K, (S(I),I=16,20) 00003530 T1 = .5 * A1 00003540 T2 = .5 * A2 00003550 T3 = .5 * A3 00003560 R1 = T2 * A2 + T3 * A3 00003570 R5 = T3 * A3 + T1 * A1 00003580 R9 = T1 * A1 + T2 * A2 00003590 R2 = T1 * A2 + A3 00003600 R4 = T1 * A2 - A3 00003610 R3 = T3 * A1 - A2 00003620 R7 = T3 * A1 + A2 00003630 R6 = T2 * A3 + A1 00003640 R8 = T2 * A3 - A1 00003650 T8 = 1. + .5 * R9 + T3 * T3 00003660 DO 2201 I=1,9 00003670 2201 R(I) = R(I) / T8 00003680 R1 = 1. - R1 00003690 R5 = 1. - R5 00003700 R9 = 1. - R9 00003710 J = 9 00003720 GO TO (2220, 2210, 2210), K 00003730 2210 DO 2211 I = 1,J,3 00003740 T1 = AR(I) 00003750 T2 = AR(I+1) 00003760 T3 = AR(I+2) 00003770 AR(I) = R1 * T1 + R4 * T2 + R7 * T3 00003780 AR(I+1) = R2 * T1 + R5 * T2 + R8 * T3 00003790 2211 AR(I+2) = R3 * T1 + R6 * T2 + R9 * T3 00003800 IF (J-9) 2999, 2220, 1510 00003810 2220 BY = BY + DBY 00003820 BZ = BZ + DBZ 00003830 GO TO (1150, 1300, 1400), K 00003840 2230 GO TO (2231, 2232), MOD 00003850 2231 X1 = U1 00003860 Y1 = V1 00003870 Z1 = 1. 00003880 GO TO 2234 00003890 2232 GO TO (2231, 2233), KD 00003900 2233 X1 = AL1 * U1 + AL4 * V1 + AL7 00003910 Y1 = AL2 * U1 + AL5 * V1 + AL8 00003920 Z1 = AL3 * U1 + AL6 * V1 + AL9 00003930 2234 T1 = Y1 * Z2 - Z1 * Y2 00003940 T2 = Z1 * X2 - X1 * Z2 00003950 T3 = X1 * Y2 - Y1 * X2 00003960 T9 = T1 * T1 + T2 * T2 + T3 * T3 00003970 T4 = BY * Z2 - BZ * Y2 00003980 T5 = BZ * X2 - BX * Z2 00003990 T6 = BX * Y2 - BY * X2 00004000 T7 = (T4 * T1 + T5 * T2 + T6 * T3) / T9 00004010 WANT = (BX * T1 + BY * T2 + BZ * T3) / T9 00004020 T8 = 0.5 * WANT 00004030 T4 = T7 * X1 + T8 * T1 00004040 T5 = T7 * Y1 + T8 * T2 00004050 T6 = T7 * Z1 + T8 * T3 00004060 GO TO (1423, 1610), KD 00004070 2901 K5 = K7 00004080 2902 WRITE (IPR,6) K8, K5 00004090 2910 READ (IRCD,3) ID2 00004100 2911 IF (ID2) 1000, 2999, 2910 00004110 2903 WRITE (IPR,6) K8, ID1, ID2 00004120 GO TO 2910 00004130 2904 WRITE (IPR,6) K8, ID2, ID(K1) 00004140 GO TO 2911 00004150 2999 CONTINUE 00004160 C2999 ITEM1=KLOCK(J) 00004170 C SEC = (ITEM1 - ITEM2) / 100. 00004180 C WRITE (IPR,94) SEC 00004190 CALL EXIT 00004200 RETURN END 00004210