C C C ************************************************************* C * * C * PROGRAM *GEBAT* * C * * C * GENERAL BUNDLE ADJUSTMENT TRIANGULATION * C * * C * BY : SABRY EL-HAKIM..1978/80 * C * ---------------------------- * C * * C * * C * PROGRAM FEATURES * C * ---------------- * C * -BUNDLE ADJUSTMENT WITH SELF CALIBRATION USING * C * HARMONIC FUNCTION * C * -BLOCK INVARIANT APPROACH * C * -COMBINED ADJUSTMENT WITH 3-D GEODESY MATH. MODEL * C * ACCEPTING : * C * -SLOPE DISTANCES * C * -HORIZONTAL DIRECTIONS * C * -VERTICAL ANGLES * C * -AST. AZIMUTHS * C * -AST. LATITUDES * C * -AST. LONGITUDES * C * * C * * C * LIMITATIONS : * C * ------------- * C * -MAX. NO. OF PHOTOS....................25 * C * -MAX. NO. OF IMAGE PTS................750 * C * -MAX. NO. OF OBJECT PTS...............250 * C * -MAX. NO. OF CONTROL PTS(V OR H).......99 * C * -MAX. NO. OF GEODETIC PTS..............23 * C * -UNLIMITED ; PTS./PHOTO * C * & GEODETIC OBSERVATIONS * C * -CASE OF DISTANCES ONLY : * C * MAX. NO. OF GEODETIC PTS. ............33 * C * MAX. NO. OF DISTANCES...............1089 * C * * C ************************************************************* C C IMPLICIT REAL*8(A-H,O-Z) MAIN 1 INTEGER Q 2 DIMENSION L(751),XP(751),YP(751),LG(751),XG(751),YG(751),ZG(751) 3 *,NF(99),PX2(3,750),LC(250),NV(99),XF(751),YF(751),LQ(250) 4 *,FF(10),PF(10) 5 REAL*8 AP1(2,150),AP1U(2,11),MP(2,2),MPA(2,150),W(2),UP1(161), 6 *NP11(6,150),NP11U(161,11),X(161),M(3,3),BP(2,2),AP2(2,750) 7 *,NP21(750,006),NP22(3,750),UP2(750),X1(33,33),X2(33,33),X3(33,33), 8 *X4(33,33),X5(33),X6(33),XN(33),XW(33),H(33),AG(1,161),P(275), 9 *NG(161,161),AN(33),UG(161),WG(275),XO(161),NE(161,1),NOG(99,99), 10 *VEC(750),R(1,161),B(161),Y(161),PX1(6),PX1U(11),XP21(750,6), 11 *XC(250),YC(250),ZC(250),XQ(250),YQ(250),ZQ(250),NP2U(750,11) 12 *,AL(751),ALG(751) 13 1 FORMAT(10A8) 14 12 FORMAT(I10) 15 30 FORMAT(5X,I6,19X,2F10.3) 16 1234 FORMAT(/10X,'PHOTO NO.',I6/9X,17('-')/) 17 39 FORMAT(I15,2F15.3) 18 47 FORMAT(I10,3F20.9) 19 409 FORMAT(3F10.4) 20 489 FORMAT(3F10.5) 21 201 FORMAT(6F10.4/4F20.9) 22 410 FORMAT(4I10) 23 419 FORMAT(3F20.4) 24 3030 FORMAT(//5X,'TOTAL NO. OF PHOTO PTS. =',I4/5X, 25 * 'TOTAL NO. OF PHOTOS =',I4////5X, 26 *'THE CONTROL POINTS COORDS. :-'/5X,29('-')/) 27 3103 FORMAT(/5X,'****ERROR****POINT',I6,'HAS NO APPROX. COORDS.****') 28 IRC=5 IWP=6 READ(IRC,410)IL1,IL2,IL3,IJKL 31 READ(IRC,410)IG,IQ 32 READ(IRC,410)NITR,Q 33 IS1=0 34 IS2=0 35 IS3=0 36 IGN=99 37 CCC=0.00009 38 2345 FORMAT(5X,16('*')/5X,'*',14X,'*'/5X,'*',' READ-IN DATA ','*'/5X, 39 *'*',14X,'*'/5X,16('*')/) 40 WRITE(IWP,2345) 41 READ(IRC,1)(PF(I),I=1,10) 42 READ(IRC,1)(FF(I),I=1,10) 43 IB=IG 44 IK=0 45 J=1 46 I=1 47 16 K=1 48 C 49 C READING THE PHOTO NO. NF AND POINT NO. L AND PHOTO COORD. XP & YP 50 C 51 READ(IRC,PF)NFIII,XPP,YPP 52 WRITE(IWP,1234)NFIII 53 NFIII=J 54 14 READ(IRC,FF)L(I),XP(I),YP(I) 55 IF(L(I).EQ.666666)GO TO 13 56 IF(L(I).EQ.-66)GO TO 15 57 WRITE(IWP,39)L(I),XP(I),YP(I) 58 XP(I)=XP(I)-XPP 59 YP(I)=YP(I)-YPP 60 XP(I)=XP(I)/1000. 61 YP(I)=YP(I)/1000. 62 LOIO=L(I) 63 IF(IG.EQ.0)GO TO 987 64 ICHEK=L(I)-L(I)/10*10 65 IF(ICHEK.EQ.IQ)GO TO 985 66 987 L(I)=L(I)+NFIII*1000000 67 III=0 68 IF(I.EQ.1)GO TO 190 69 I2=I-1 70 DO 913 II=1,I2 71 LL=(AL(II)+0.01)/1000000000 72 A1L=LL 73 AR=AL(II)-A1L*1000000000.0+0.01 74 IR=AR 75 LOIIO=IR/1000 76 IR=IR-IR/1000*1000 77 IF(LOIO.EQ.LOIIO)IM=IR 78 IF(LOIO.EQ.LOIIO)III=1 79 913 CONTINUE 80 190 IF(III.EQ.0)IB=IB+1 81 IF(III.EQ.0)IM=IB 82 AL(I)=L(I) 83 AL(I)=AL(I)*1000+IM 84 GO TO 986 85 985 CONTINUE 86 L(I)=L(I)+NFIII*1000000 87 III=0 88 IF(I.EQ.1)GO TO 19 89 I2=I-1 90 DO 113 II=1,I2 91 LL=(AL(II)+0.01)/1000000000 92 A1L=LL 93 AR=AL(II)-A1L*1000000000.0+0.01 94 IR=AR 95 LOIIO=IR/1000 96 IR=IR-IR/1000*1000 97 IF(LOIO.EQ.LOIIO)IM=IR 98 IF(LOIO.EQ.LOIIO)III=1 99 113 CONTINUE 100 19 IF(III.EQ.0)IK=IK+1 101 IF(III.EQ.0)IM=IK 102 AL(I)=L(I) 103 AL(I)=AL(I)*1000+IM 104 986 CONTINUE 105 I=I+1 106 K=K+1 107 GO TO 14 108 13 J=J+1 109 GO TO 16 110 15 NP=I-1 111 C 112 C NP : THE TOTAL NO. OF PHOTO POINTS 113 C 114 NFI=J 115 C 116 C NFI : THE NUMBER OF PHOTOS 117 C 118 WRITE(IWP,3030)NP,NFI 119 READ(IRC,1)(PF(I),I=1,10) 120 I=0 121 55 I=I+1 122 READ(IRC,PF)LC(I),XC(I),YC(I),ZC(I) 123 IF(LC(I).NE.-66)GO TO 55 124 IC=I-1 125 56 FORMAT(I10) 126 DO 762J=1,NP 127 LL=(AL(J)+0.01)/1000000000 128 A1L=LL 129 AR=AL(J)-A1L*1000000000.0+0.01 130 IR=AR 131 IR=IR/1000 132 DO 31 I=1,IC 133 IF(IR.EQ.LC(I))GO TO 32 134 GO TO 3102 135 32 XG(J)=XC(I) 136 YG(J)=YC(I) 137 ZG(J)=ZC(I) 138 ALG(J)=AL(J) 139 GO TO 762 140 3102 IF(I.EQ.IC)WRITE(IWP,3103)IR 141 31 CONTINUE 142 762 CONTINUE 143 C 144 C READING THE WT. MATRIX OF THE PHOTO POINTS & PX1, PX1U&PX2 145 C 146 READ(IRC,24)P1,P2,P3 147 24 FORMAT(3F10.3) 148 PX=P2/(P1*P2-P3**2) 149 PY=P1/(P1*P2-P3**2) 150 PXY=P3/(P1*P2-P3**2) 151 JC=IC 152 NU=6 153 N1=11 154 READ(IRC,201)(PX1(I),I=1,NU),PX1U(1),PX1U(2),PX1U(3),P8 155 DO 359 I=4,N1 156 359 PX1U(I)=P8 157 READ(IRC,489)PXC,PYC,PZC 158 READ(IRC,419)PXCO,PYCO,PZCO 159 NUN=NU*NFI 160 MBW=NUN 161 NAN=NUN+N1 162 DO 21 I=1,NAN 163 B(I)=0.0 164 21 X(I)=0.0 165 DO 98 I=1,NFI 166 READ(IRC,47)ND,XX,YY,ZZ 167 READ(IRC,489)AOM,PHI,AKP 168 NUII=NU*I 169 X(NUII-2)=AOM 170 X(NUII-1)=PHI 171 X(NUII)=AKP 172 X(NUII-5)=XX 173 X(NUII-4)=YY 174 98 X(NUII-3)=ZZ 175 READ(IRC,198)XO2,YO2,F2 176 XO2=XO2/1000. 177 YO2=YO2/1000. 178 F2=F2/1000. 179 198 FORMAT(3F10.3) 180 X(NUN+1)=XO2 181 X(NUN+2)=YO2 182 X(NUN+3)=F2 183 ITR=0 184 JCS=JC*3 185 DO 3698 I=1,3 186 DO 3698 J=1,JCS 187 VEC(J)=0.0 188 3698 PX2(I,J)=0.0 189 READ(IRC,56)ICON 190 READ(IRC,2121)(NF(I),I=1,ICON) 191 2121 FORMAT(8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10/8I10) 192 DO 6789 I=1,ICON 193 INF=NF(I) 194 DO 982 J=1,IC 195 982 IF(LC(J).EQ.NF(I))INP=J 196 6789 WRITE(IWP,39)INF,XC(INP),YC(INP) 197 READ(IRC,56)IVER 198 READ(IRC,2121)(NV(I),I=1,IVER) 199 DO 6790 I=1,IVER 200 INV=NV(I) 201 DO 981 J=1,IC 202 981 IF(LC(J).EQ.NV(I))INF=J 203 6790 WRITE(IWP,4747)INV,ZC(INF) 204 4747 FORMAT(I15,30X,F15.3) 205 READ(IRC,12)ICC 206 DO 879 I=1,ICC 207 879 READ(IRC,47)LQ(I),XQ(I),YQ(I),ZQ(I) 208 IDOR=0 209 501 CONTINUE 210 ICHT=0 211 IF(ITR.LE.2)GO TO 1529 212 IF(IS1.EQ.0)GO TO 1529 213 DO 1527 I=1,3 214 DO 1527 J=1,JCS 215 1527 PX2(I,J)=NP22(I,J) 216 GO TO 1531 217 1529 DO 524 I=1,JC 218 PX2(1,3*I-2)=PXC 219 PX2(2,3*I-1)=PYC 220 524 PX2(3,3*I)=PZC 221 1531 CONTINUE 222 DO 1122 I=1,ICON 223 DO 980 J=1,NP 224 LL=(AL(J)+0.01)/1000000000 225 A1L=LL 226 AR=AL(J)-A1L*1000000000.0+0.01 227 JR=AR 228 IR=JR/1000 229 980 IF(IR.EQ.NF(I))INF=JR-IR*1000 230 PX2(1,3*INF-2)=PXCO 231 1122 PX2(2,3*INF-1)=PYCO 232 DO 3344 I=1,IVER 233 DO 983 J=1,NP 234 LL=(AL(J)+0.01)/1000000000 235 A1L=LL 236 AR=AL(J)-A1L*1000000000.0+0.01 237 JR=AR 238 IR=JR/1000 239 983 IF(IR.EQ.NV(I))INV=JR-IR*1000 240 3344 PX2(3,3*INV)=PZCO 241 DO 142 I=1,NU 242 DO 142 J=1,NUN 243 142 NP11(I,J)=0.0 244 DO 143 I=1,NAN 245 UP1(I)=0.0 246 DO 143 J=1,N1 247 143 NP11U(I,J)=0.0 248 9182 CONTINUE 249 ITR=ITR+1 250 INB=NP 251 NPS=INB 252 NP=751 253 INU=NUN 254 INC=JCS 255 JCS=750 256 NUN=150 257 NAN=161 258 C 259 C CALLING PNORM1 TO FORM NP11, NP21, NP22, UP1 AND UP2 260 C 261 CALL PNORM1(L,XP,YP,XG,YG,ZG,X,NU,NP,NFI,AP1,AP1U,MP,W,NP11U, 262 *NP11,INB,MPA,N1,NUN,NAN,UP1,LG,M,BP,PX,PY,PXY,AP2,NP21,JCS, 263 *NP22,UP2,MBW,B,VEC,ITR,XF,YF,IDOR,IS2,INC,INU,PX1,NP2U,AL,ALG) 264 NP=INB 265 IF(ITR.GE.Q)GO TO 629 266 DO 634 I=1,N1 267 DO 635 J=1,NAN 268 635 NP11U(J,I)=0.0 269 DO 636 J=1,JCS 270 636 NP2U(J,I)=0.0 271 K=INU+I 272 634 UP1(K)=0.0 273 629 CONTINUE 274 IF(IDOR.EQ.1)GO TO 1000 275 IF(IG.EQ.0)GO TO 228 276 IF(IJKL.EQ.1)GO TO 227 277 IF(ITR.GT.1)GO TO 229 278 NM=161 279 NMO=161 280 NN=275 281 NP=751 282 CALL GNORM(X1,X2,X3,X4,X5,X6,XN,XW,H,AG,P,NG,AN,UG,WG,NM,NN,NO1, 283 *NO2,NO3,NO4,NO5,NO6,33,XO,NMO,LC,ITR,NNO,LG,INB,NP,ALG) 284 NM=INB 285 INB=NPS 286 NP=INB 287 GO TO 226 288 227 IGM=33 289 NP=751 290 CALL SNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,ITR,X1,X2,INB,IGM,L,XO,IGN, 291 *AG,WG,AL) 292 NP=INB 293 GO TO 228 294 226 CONTINUE 295 REWIND8 296 DO 328 J=1,NM 297 328 WRITE(8)(NG(I,J),I=1,NM) 298 ENDFILE8 299 GO TO 228 300 229 REWIND8 301 DO 329 J=1,NM 302 329 READ(8)(NG(I,J),I=1,NM) 303 228 CONTINUE 304 JCS=750 305 NUN=150 306 NAN=161 307 IF(IJKL.EQ.1.AND.IG.NE.0)CALL NORM2(NP22,NOG,PX2,NM,JCS,JC,IG,Y,M, 308 *MGW,IGN,3,INC,INU) 309 IF(IJKL.EQ.1.AND.IG.NE.0)GO TO 7723 310 CALL NORM2(NP22,NG,PX2,NM,JCS,JC,IG,Y,M,MGW,NAN,2,INC,INU) 311 7723 CONTINUE 312 JCS=750 313 NUN=150 314 NAN=161 315 IF(IJKL.EQ.1.AND.IG.NE.0)CALL CVEC(UP1,NOG,NP21,UG,UP2,NP22,P,JC, 316 *IGN,JCS,IG,NNO,VEC,3,INC,INU,NP2U) 317 IF(IJKL.EQ.1.AND.IG.NE.0)GO TO 7724 318 CALL CVEC(UP1,NG,NP21,UG,UP2,NP22,P,JC,NAN,JCS,IG,NNO,VEC,1, 319 *INC,INU,NP2U) 320 7724 CONTINUE 321 ICR=0 322 DO 7714 I=1,NFI 323 DO 7714 J=1,NU 324 ICR=ICR+1 325 7714 UP1(ICR)=UP1(ICR)-B(ICR)*PX1(J) 326 IF(ITR.LT.Q)GO TO 6715 DO 6714 I=1,N1 ICR=ICR+1 6714 UP1(ICR)=UP1(ICR)-B(ICR)*PX1U(I) 6715 CONTINUE JCS=750 327 NUN=150 328 NAN=161 329 CALL NORMAL(JC,JCS,NUN,NAN,NG,NP22,NP11,NP11U,PX1,PX1U,NE,IG,VEC, 330 *NP21,INC,INU,NP2U,XP21,NOG,IGN,IJKL) 331 JCS=INC 332 NUN=INU 333 NAN=NUN+N1 334 DO 1420 I=1,NAN 335 DO 1420 J=I,NAN 336 1420 NG(I,J)=NG(J,I) 337 NAN=161 338 CALL CHLSK(NAN,NG,R,INDEF,ITR,N1,INU,Q) 339 NAN=INU+11 340 DO 207 I=1,NAN 341 DO 207 J=1,NAN 342 207 NG(I,J)=0.0 343 REWINDIL1 344 IF(ITR.LT.Q)NAN=NUN 345 DO 205 I=1,NAN 346 READ(IL1)(R(1,II),II=1,NAN) 347 DO 206 J=I,NAN 348 206 NG(I,J)=R(1,J) 349 205 CONTINUE 350 NAN=161 351 CALL CHLSKS(NAN,NG,UP1,B,Y,ITR,N1,INU,Q) 352 NAN=NUN+11 353 DO 427 I=1,NAN 354 IF(B(I).GT.CCC)ICHT=1 355 IF(ITR.LT.Q.AND.I.GT.NUN)GO TO 22 356 427 X(I)=X(I)+B(I) 357 22 WRITE(IWP,330) 358 DO 320 I=1,NFI 359 M1=I*6-5 360 M2=M1-1+2 361 M3=M2-1+2 362 M4=M3+1 363 M5=M4+1 364 M6=M5+1 365 WRITE(IWP,340)I,X(M1),X(M2),X(M3),X(M4),X(M5),X(M6) 366 320 WRITE(IWP,341)B(M1),B(M2),B(M3),B(M4),B(M5),B(M6) 367 340 FORMAT(I3,'-ADJ.PAR.=',3F12.3,3F10.4) 368 341 FORMAT(4X,'RESIDUAL=',3F12.3,3F10.4//) 369 330 FORMAT(20X,'XC',10X,'YC',10X,'ZC',7X,'OMEGA',7X,'PHI',5X,'KAPPA'/) 370 I6=NUN+1 371 WRITE(IWP,342) 372 DO 345 I=I6,NAN 373 345 WRITE(IWP,350)X(I),B(I) 374 350 FORMAT(5X,2F15.9) 375 342 FORMAT(5X,'THE BASIC PAR. + HARMONIC PAR.'/5X,32('-')) 376 REWINDIL2 377 DO 40 I=1,JCS 378 40 P(I)=0.0 379 IK=0 380 DO 402 J=1,NFI 381 DO 402 I=1,NU 382 IK=IK+1 383 READ(IL2)(NP21(KK,I),KK=1,JCS) 384 DO 402 K=1,JCS 385 402 P(K)=P(K)+NP21(K,I)*B(IK) 386 IK=IK+1 387 DO 45 I=1,JCS 388 IO=0 389 DO 45 J=IK,NAN 390 IO=IO+1 391 45 P(I)=P(I)+NP2U(I,IO)*B(J) 392 IF(IG.EQ.0.OR.IJKL.EQ.1)GO TO 1000 393 REWIND8 394 DO 707 J=1,NM 395 707 READ(8)(NG(I,J),I=1,NM) 396 1000 CONTINUE 397 DO 708 J=1,3 398 DO 708 I=1,JCS 399 708 PX2(J,I)=0.0 400 DO 709 I=1,JCS 401 709 UP2(I)=UP2(I)+P(I) 402 IF(IDOR.EQ.1)GO TO 1002 403 JCS=750 404 NUN=150 405 NAN=161 406 IF(IJKL.EQ.1.AND.IG.NE.0)GO TO 1071 407 CALL NORM2(NP22,NG,PX2,NM,JCS,JC,IG,Y,M,MGW,NAN,1,INC,INU) 408 1071 CONTINUE 409 JCS=750 410 NUN=150 411 NAN=161 412 1002 CONTINUE 413 IF(IJKL.EQ.1.AND.IG.NE.0)CALL CVEC(UP1,NOG,NP21,UG,UP2,NP22,P,JC, 414 *IGN,JCS,IG,NNO,VEC,4,INC,INU,NP2U) 415 IF(IJKL.EQ.1.AND.IG.NE.0)GO TO 1072 416 CALL CVEC(UP1,NG,NP21,UG,UP2,NP22,P,JC,NAN,JCS,IG,NNO,VEC,2, 417 *INC,INU,NP2U) 418 1072 CONTINUE 419 JCS=INC 420 NUN=INU 421 NAN=NUN+N1 422 WRITE(IWP,20) 423 20 FORMAT(///5X,'THE COORDS. RESIDUALS'/5X,22('*')) 424 3330 FORMAT(/I10,'-',3F12.3/11X,3F12.3) 425 DO 780 I=1,NP 426 LL=(AL(I)+0.01)/1000000000 427 A1L=LL 428 AR=AL(I)-A1L*1000000000.0+0.01 429 KN=AR 430 KN=KN-KN/1000*1000 431 LGOIO=AR/1000 432 XG(I)=XG(I)-VEC(3*KN-2) 433 YG(I)=YG(I)-VEC(3*KN-1) 434 ZG(I)=ZG(I)-VEC(3*KN) 435 IF(I.EQ.1)GO TO 1003 436 I2=I-1 437 DO 1004 JJ=1,I2 438 LL=(ALG(I)+0.01)/1000000000 439 A1L=LL 440 AR=ALG(I)-A1L*1000000000.0+0.01 441 LGOIO=AR 442 LGOIO=LGOIO/1000 443 LL=(ALG(JJ)+0.01)/1000000000 444 A1L=LL 445 AR=ALG(JJ)-A1L*1000000000.0+0.01 446 LGOJJO=AR 447 LGOJJO=LGOJJO/1000 448 IF(LGOIO.EQ.LGOJJO)GO TO 780 449 1004 CONTINUE 450 1003 WRITE(IWP,3330)LGOIO,XG(I),YG(I),ZG(I),VEC(3*KN-2), 451 1VEC(3*KN-1),VEC(3*KN) 452 780 CONTINUE 453 IF(IG.EQ.0.OR.IJKL.EQ.0)GO TO 881 454 DO 882 I=1,IG 455 I3=I*3 456 XO(I3)=XO(I3)-VEC(I3) 457 XO(I3-1)=XO(I3-1)-VEC(I3-1) 458 882 XO(I3-2)=XO(I3-2)-VEC(I3-2) 459 881 CONTINUE 460 NP=751 461 CALL CHECK(LG,XG,YG,ZG,LQ,XQ,YQ,ZQ,NQ,NP,INB,ICC,ALG) 462 NP=INB 463 IF(IS3.EQ.0)GO TO 5011 464 IF(ICHT.EQ.0)IDOR=1 465 IF(ICHT.EQ.0)GO TO 9182 466 5011 IF(ITR.LT.NITR)GO TO 501 467 IF(IG.EQ.0)STOP NP=751 CALL SNORM(NOG,UG,NAN,NP,XG,YG,ZG,IG,0,X1,X2,INB,IGM,L,XO,IGN,AG, *WG,AL) STOP 468 END 469 SUBROUTINE CHECK(L,X,Y,Z,LC,XC,YC,ZC,IC,NB,INB,IN,AL) 470 IMPLICIT REAL*8 (A-H,O-Z) 471 REAL*8 X(NB),Y(NB),Z(NB),XC(250),YC(250),ZC(250) 472 *,AL(NB) 473 DIMENSION L(NB),LC(250) 474 IWP=6 NP=INB 476 SX=0. 477 SY=0.0 478 SZ=0.0 479 WRITE(IWP,17) 480 17 FORMAT(//5X,'CHECK PT. DIFFS.'/5X,16('-')/) 481 DO 30 I=1,IN 482 J=LC(I) 483 A=XC(I) 484 B=YC(I) 485 C=ZC(I) 486 DO 980 II=1,NP 487 LL=(AL(II)+0.01)/1000000000 488 A1L=LL 489 AR=AL(II)-A1L*1000000000.0+0.01 490 JR=AR 491 IR=JR/1000 492 980 IF(LC(I).EQ.IR)K=II 493 AA=X(K)-A 494 AD=AA**2 495 SX=SX+AD 496 BB=Y(K)-B 497 BD=BB**2 498 SY=SY+BD 499 CC=Z(K)-C 500 WRITE(IWP,22)J,AA,BB,CC 501 CC=CC**2 502 SZ=SZ+CC 503 30 CONTINUE 504 ZN=IN 505 RMX=DSQRT(SX/ZN) 506 RMY=DSQRT(SY/ZN) 507 RMZ=DSQRT(SZ/ZN) 508 WRITE(IWP,21)RMX,RMY,RMZ 509 22 FORMAT(I10,3F10.3) 510 21 FORMAT(///5X,'RMS =',3F10.4) 511 RETURN 512 END 513 SUBROUTINE CVEC(UP1,NG,NP21,UG,UP2,NP22,Y,JC,NAM,JJS,IG,NNO,VEC,L, 516 *INC,INU,NP2U) 517 IMPLICIT REAL*8(A-H,O-Z) 518 REAL*8 Y(275),UP1(NAM),UP2(JJS),NG(NAM,NAM),UG(NAM),NP22(3,JJS) 519 *,NP21(JJS,006),VEC(JJS),NP2U(JJS,11) 520 NU=6 521 IL2=22 522 JCS=INC 523 NUN=INU 524 NAN=NUN+11 525 IF(IG.EQ.0)GO TO 222 526 II=7*IG-NNO 527 IF(L.GE.3)II=3*IG 528 DO 15 I=1,II 529 15 Y(I)=UG(I) 530 DO 10 K=1,IG 531 DO 10 I=1,3 532 IP=3*K-3+I 533 JG=5*K-5+I 534 IF(L.GE.3)JG=IP 535 10 Y(JG)=Y(JG)+UP2(IP) 536 DO 20 I=1,IG 537 DO 20 K=1,3 538 JG=5*I-5+K 539 JP=3*I-3+K 540 IF(L.GE.3)JG=JP 541 VEC(JP)=0.0 542 DO 25 J=1,II 543 25 VEC(JP)=VEC(JP)+NG(JG,J)*Y(J) 544 20 CONTINUE 545 222 IX=IG+1 546 DO 30 I=IX,JC 547 DO 30 J=1,3 548 JP=3*I-3+J 549 VEC(JP)=0.0 550 DO 31 K=1,3 551 KP=3*I-3+K 552 31 VEC(JP)=VEC(JP)+NP22(J,KP)*UP2(KP) 553 30 CONTINUE 554 IF(L.EQ.2.OR.L.EQ.4)RETURN 555 NF=NUN/6 556 REWINDIL2 557 II=0 558 DO 40 J=1,NF 559 DO 40 I=1,NU 560 READ(IL2)(NP21(KK,I),KK=1,JCS) 561 II=II+1 562 DO 50 K=1,JCS 563 50 UP1(II)=UP1(II)-NP21(K,I)*VEC(K) 564 40 CONTINUE 565 II=II+1 566 IO=0 567 DO 60 I=II,NAN 568 IO=IO+1 569 DO 60 J=1,JCS 570 60 UP1(I)=UP1(I)-NP2U(J,IO)*VEC(J) 571 RETURN 572 END 573 SUBROUTINE CHLSK(M,A,R,INDEF,ITR,N1,INU,Q) 574 C 575 C ************************************************************************ 576 C *SUBROUTINE CHLSK * 577 C * * 578 C * THIS SUBROUTINE COMPUTES THE UPPER TRIANGULAR MATRIX R * 579 C * WHICH, WHEN MULTIPLIED BY ITS OWN TRANSPOSE, YIELDS * 580 C * THE GIVEN SYMMETRIC MATRIX A, BY THE CHOLESKY METHOD * 581 C * * 582 C * THE MATRIX R IS STORED ( ROW BY ROW ) IN AN ON LINE DISK * 583 C * * 584 C * MODIFIED BY S. EL-HAKIM AUG9,77 * 585 C ************************************************************************ 586 C 587 IMPLICIT REAL*8(A-H,O-Z) 588 INTEGER Q 589 REAL*8 A(M,M),R(1,M) 590 IWP=6 IL1=21 592 N=INU+N1 593 IF(ITR.LT.Q)N=N-N1 594 REWINDIL1 595 INDEF=0 596 DO 30 IP=1,N 597 IF(A(IP,IP))05,05,10 598 5 WRITE(IWP,6)IP 599 6 FORMAT(/5X,'ARRAY NOT POSITIVE DEFINITE,COMP#',I4,'IS ADDED TO 1.0 600 *D-04'/5X,56('$')/) 601 A(IP,IP)=1.0D-08 602 10 R(01,IP)=DSQRT(A(IP,IP)) 603 IF(IP-N)12,40,40 604 12 IPPLS1=IP+1 605 DO 20 K=IPPLS1,N 606 R(01,K)=A(IP,K)/R(01,IP) 607 20 CONTINUE 608 DO 31 I=IPPLS1,N 609 DO 31 K=I,N 610 31 A(I,K)=A(I,K)-R(01,I)*R(01,K) 611 WRITE(IL1)(R(1,II),II=1,N) 612 30 CONTINUE 613 40 WRITE(IL1)(R(1,II),II=1,N) 614 ENDFILEIL1 615 RETURN 616 END 617 SUBROUTINE CHLSKS(M,R,B,X,Y,ITR,N1,INU,Q) 618 C ************************************************************************ 619 C *SUBROUTINE CHLSKS * 620 C * * 621 C * THIS SUBROUTINE SOLVES VECTOR EQUATION AX+B=0 FOR * 622 C * VECTOR X, WHERE VECTOR B IS GIVEN AND WHERE R, THE * 623 C * CHOLESKY DECOMPOSITION OF THE POSITIVE DEFINITE MATRIX * 624 C * A, IS KNOWN.. * 625 C * * 626 C ************************************************************************ 627 C 628 IMPLICIT REAL*8(A-H,O-Z) 629 INTEGER Q 630 REAL*8 R(M,M),B(M),X(M),Y(M) 631 N=INU+N1 632 IF(ITR.LT.Q)N=N-N1 633 DO 30 K=1,N 634 S=B(K) 635 IF(K-1)25,25,10 636 10 KLSS1=K-1 637 DO 20 I=1,KLSS1 638 S=S+R(I,K)*Y(I) 639 20 CONTINUE 640 25 Y(K)=-S/R(K,K) 641 30 CONTINUE 642 DO 50 IDUM=1,N 643 I=N+1-IDUM 644 S=Y(I) 645 IF(N-I)45,45,34 646 34 IPLS1=I+1 647 DO 40 K=IPLS1,N 648 S=S-R(I,K)*X(K) 649 40 CONTINUE 650 45 X(I)=S/R(I,I) 651 50 CONTINUE 652 RETURN 653 END 654 SUBROUTINE GNORM(X1,X2,X3,X4,X5,X6,XN,XW,H,AG,P,Q,AN,UG,WG,NI,NJ, 655 *NO1,NO2,NO3,NO4,NO5,NO6,MG,XO,MMO,LC,ITR,NNO,LK,INB,NQ,ALK) 656 C 657 C ************************************************************************ 658 C *SUBROUTINE GNORM * 659 C * * 660 C * THIS SUBROUTINE FORMS THE GEODETIC NORMAL EQUATIONS * 661 C * (( NG )) AND THE VECTOR UG * 662 C * * 663 C * WRITTEN BY S. EL-HAKIM JULY20,77 * 664 C ************************************************************************ 665 C 666 C 667 IMPLICIT REAL*8(A-H,O-Z) 668 DIMENSION X1(MG,MG),X2(MG,MG),X3(MG,MG),X4(MG,MG),X5(MG),X6(MG), 669 1XN(MG),XW(MG),H(MG),AG(1,NI),P(NJ),Q(NI,NI),AN(MG),UG(MMO),WG(NJ) 670 2,XO(MMO),LC(250),LK(NQ),ALK(NQ) 671 IRC=5 NP=INB 673 IF(ITR.GT.1)GO TO 555 674 C 675 C READING : IG...NO. OF GEODETIC POINTS 676 C NO1..NO. OF OBSERVED DISTANCES 677 C NO2..NO. OF OBSERVED VERTICAL ANGLEA 678 C NO3..NO OF OBSERVED DIRECTIONS 679 C NO4..NO. OF OBSERVED AZIMUTHS 680 C NO5..NO. OF OBSERVED LATITUDES 681 C NO6..NO. OF OBSERVED LONGITUDES 682 C 683 READ(IRC,10)NG,NO1,NO2,NO3,NO4,NO5,NO6 684 10 FORMAT(7I10) 685 C 686 C NM : NO. OF UNKNOWNS (INCLUDING INITIAL AZ. AND 687 C VERTICAL ANGLE REFRACTION COOF.) 688 C NN : NO. OF GEODETIC OBSERVATIONS 689 C 690 NMO=7*NG 691 C 692 C NNO : NO. OF POINTS WITHOUT DIRECTIONS OBSERVED 693 C 694 READ(IRC,12)NNO 695 12 FORMAT(I11) 696 DO 329 I=1,NMO 697 329 XO(I)=0.0 698 DO 20 I=1,NG 699 DO 20 J=1,NG 700 X1(I,J)=0.0 701 X2(I,J)=0.0 702 X3(I,J)=0.0 703 X4(I,J)=0.0 704 X5(J)=0.0 705 X6(J)=0.0 706 20 CONTINUE 707 C 708 C READING THE GEODETIC OBSERVATIONS 709 C X1 : DISTANCE < <................O1 710 C X2 : VERTICAL ANGLE < <..........O2 711 C X3 : DIRECTION < <...............O3 712 C X4 : AST. AZIMUTH < <............O4 713 C 714 C ** OBSERVATIONS X1,X2,X3&X4 ARE FROM I TO J ** 715 C 716 C X5 : AST. LATITUDE 717 C X6 : AST. LONGITUDE 718 C ------------------------------------- 719 C O5 : INITIAL AZIMUTH 720 C I5 : SERIAL NO. FOR INITIAL AZIMUTH 721 C 722 600 CONTINUE 723 READ(IRC,9)K,M,O1,O2,O3,O4,O5,I5 724 IF(K.EQ.-9)GO TO 808 725 DO 980 II=1,NP 726 LL=(ALK(II)+0.01)/1000000000 727 A1L=LL 728 AR=ALK(II)-A1L*1000000000.0+0.01 729 JR=AR 730 IR=JR/1000 731 IF(K.EQ.IR)I=JR-IR*1000 732 980 IF(M.EQ.IR)J=JR-IR*1000 733 X1(I,J)=O1 734 X2(I,J)=O2 735 X3(I,J)=O3 736 X4(I,J)=O4 737 IUI=6*NG+I5 738 XO(IUI)=O5 739 GO TO 600 740 808 CONTINUE 741 DO 29 K=1,NO5 742 C 743 C READING THE GEODETIC LATITUDE ( XN ) 744 C AND THE GEODETIC LONGITUDE( XW ) 745 C 746 DO 981 II=1,NP 747 LL=(ALK(II)+0.01)/1000000000 748 A1L=LL 749 AR=ALK(II)-A1L*1000000000.0+0.01 750 JR=AR 751 IR=JR/1000 752 981 IF(K.EQ.IR)I=JR-IR*1000 753 29 READ(IRC,33)X5(I),X6(I) 754 33 FORMAT(10X,2F20.9) 755 READ(IRC,30)A,E2 772 DO 32 K=1,NG 756 DO 982 II=1,NP 757 LL=(ALK(II)+0.01)/1000000000 758 A1L=LL 759 AR=ALK(II)-A1L*1000000000.0+0.01 760 JR=AR 761 IR=JR/1000 762 982 IF(K.EQ.IR)I=JR-IR*1000 763 READ(IRC,39)XN(I),XW(I),H(I) 764 AN(I)=A/DSQRT(DCOS(XN(I))**2+(1.0-E2)*DSIN(XN(I))**2) 765 XO(5*I-4)=(AN(I)+H(I))*DCOS(XN(I))*DCOS(XW(I)) 766 XO(5*I-3)=(AN(I)+H(I))*DCOS(XN(I))*DSIN(XW(I)) 767 XO(5*I-2)=(AN(I)*(1.0-E2)+H(I))*DSIN(XN(I)) 768 XO(5*I-1)=XN(I) 769 XO(5*I)=XW(I) 770 32 CONTINUE 771 39 FORMAT(10X,3F20.9) 773 555 CONTINUE 774 NM=NMO-NNO 775 INB=NM 776 D=1.0 777 RAD=4.0*DATAN(D) 778 NN=NO1+NO2+NO3+NO4+NO5+NO6 779 G=206265.0 780 DO 163 I=1,NN 781 163 WG(I)=0.0 782 DO 705 I=1,NMO 783 UG(I)=0.0 784 DO 705 J=1,NM 785 705 Q(I,J)=0.0 786 JH=1 787 J=NO1+1 788 NOO=NO1+NO2 789 DO 105 I=J,NOO 790 105 P(I)=0.25 791 J=NO2+1+NO1 792 MK=NO3+NO4+NO1+NO2 793 DO 106 I=J,MK 794 106 P(I)=1.0 795 J=MK+1 796 NK=NO5+NO6+MK 797 DO 107 I=J,NK 798 107 P(I)=2.0 799 999 CONTINUE 800 II=0 801 JJ=0 802 DO 100 K=1,NO1 803 DO 94 I=1,NM 804 94 AG(1,I)=0.0 805 DO 101 I=1,NG 806 DO 101 J=1,NG 807 IF(I.LT.II)GO TO 101 808 IF(I.LE.II.AND.J.LE.JJ)GO TO 101 809 IF(I.EQ.J)GO TO 101 810 IF(X1(I,J).EQ.0.0)GO TO 101 811 SC=DSQRT((XO(5*I-4)-XO(5*J-4))**2+(XO(5*I-3)-XO(5*J-3))**2 812 *+(XO(5*I-2)-XO(5*J-2))**2) 813 S=SC 814 P(K)=1.0/((0.003*SC+0.03)**2) 815 C1=(XO(5*J-4)-XO(5*I-4))/S 816 C2=(XO(5*J-3)-XO(5*I-3))/S 817 C3=(XO(5*J-2)-XO(5*I-2))/S 818 AG(1,5*I-4)=-C1*P(K) 819 AG(1,5*I-3)=-C2*P(K) 820 AG(1,5*I-2)=-C3*P(K) 821 AG(1,5*J-4)=C1*P(K) 822 AG(1,5*J-3)=C2*P(K) 823 AG(1,5*J-2)=C3*P(K) 824 WG(K)=(SC-X1(I,J))*P(K) 825 II=I 826 JJ=J 827 GO TO 701 828 101 CONTINUE 829 701 DO 601 IH=1,NM 830 DO 602 KU=1,NM 831 602 Q(IH,KU)=Q(IH,KU)+AG(JH,IH)*AG(JH,KU)/P(K) 832 601 CONTINUE 833 DO 620 IH=1,NM 834 KU=K 835 620 UG(IH)=UG(IH)+AG(1,IH)*WG(KU)/P(K) 836 100 CONTINUE 837 N=NO1+1 838 M=NO1+NO2 839 II=0 840 JJ=0 841 DO 110 K=N,M 842 DO 631 I=1,NM 843 631 AG(1,I)=0.0 844 L=0 845 NAD=0 846 DO 111 I=1,NG 847 L=L+1 848 DO 111 J=1,NG 849 IF(I.LT.II)GO TO 111 850 IF(I.LE.II.AND.J.LE.JJ)GO TO 111 851 IF(I.EQ.J)GO TO 111 852 IF(X2(I,J).EQ.0.0)GO TO 111 853 S=DSQRT((XO(5*I-4)-XO(5*J-4))**2+(XO(5*I-3)-XO(5*J-3))**2) 854 S=DSQRT(S**2+(XO(5*J-2)-XO(5*I-2))**2) 855 BC=DARSIN(((XO(5*J-4)-XO(5*I-4))*DCOS(XN(I))*DCOS(XW(I)) 856 *+(XO(5*J-3)-XO(5*I-3))*DCOS(XN(I))*DSIN(XW(I))+(XO(5*J-2) 857 *-XO(5*I-2))*DSIN(XN(I)))/S) 858 VS=(XO(5*J-3)-XO(5*I-3))*DCOS(XW(I))-(XO(5*J-4)-XO(5*I-4))* 859 *DSIN(XW(I)) 860 VC=(XO(5*J-2)-XO(5*I-2))*DCOS(XN(I))-(XO(5*J-4)-XO(5*I-4))*DSIN(XN 861 *(I))*DCOS(XW(I))-(XO(5*J-3)-XO(5*I-3))*DSIN(XN(I))*DSIN(XW(I)) 862 AZ=DATAN(VS/VC) 863 IF(VS.LT.0.0.AND.VC.LT.0.0)AZ=AZ+RAD 864 IF(VS.LT.0.0.AND.VC.GT.0.0)AZ=2.0*RAD+AZ 865 IF(VS.GT.0.0.AND.VC.LT.0.0)AZ=RAD+AZ 866 BB=XN(I) 867 BO=XW(I) 868 XN(I)=X5(I) 869 XW(I)=X6(I) 870 IF(X6(I).EQ.0.0)XW(I)=BO 871 B1=(S*DCOS(XN(I))*DCOS(XW(I))+(XO(5*I-4)-XO(5*J-4))*DSIN(X2(I,J))) 872 */(S**2*DCOS(X2(I,J))) 873 B2=(S*DCOS(XN(I))*DSIN(XW(I))+(XO(5*I-3)-XO(5*J-3))*DSIN(X2(I,J))) 874 */(S**2*DCOS(X2(I,J))) 875 B3=(S*DSIN(XN(I))+(XO(5*I-2)-XO(5*J-2))*DSIN(X2(I,J)))/(S**2* 876 *DCOS(X2(I,J))) 877 B4=DCOS(AZ) 878 B5=DCOS(XN(I))*DSIN(AZ) 879 AG(1,5*I-4)=-B1*P(K) 880 AG(1,5*I-3)=-B2*P(K) 881 AG(1,5*I-2)=-B3*P(K) 882 AG(1,5*J-4)=B1*P(K) 883 AG(1,5*J-3)=B2*P(K) 884 AG(1,5*J-2)=B3*P(K) 885 AG(1,5*I-1)=B4*P(K) 886 AG(1,5*I)=B5*P(K) 887 IINGIL=5*NG+L 888 AG(1,IINGIL)=-S/1000.0*P(K) 889 WG(K)=(BC-X2(I,J)-S*XO(IINGIL)/1000.0)*P(K) 890 XN(I)=BB 891 XW(I)=BO 892 II=I 893 JJ=J 894 GO TO 711 895 111 CONTINUE 896 711 DO 603 IH=1,NM 897 DO 604 KU=1,NM 898 604 Q(IH,KU)=Q(IH,KU)+AG(JH,IH)*AG(JH,KU)/P(K) 899 603 CONTINUE 900 DO 630 IH=1,NM 901 KU=K 902 630 UG(IH)=UG(IH)+AG(1,IH)*WG(KU)/P(K) 903 110 CONTINUE 904 N=M+1 905 M=M+NO3 906 II=0 907 JJ=0 908 L=NG 909 DO 120 K=N,M 910 DO 632 I=1,NM 911 632 AG(1,I)=0.0 912 NAD=0 913 DO 121 I=1,NG 914 DO 121 J=1,NG 915 IF(I.LT.II)GO TO 121 916 IF(I.LE.II.AND.J.LE.JJ)GO TO 121 917 IF(I.EQ.J)GO TO 121 918 IF(X3(I,J).GE.7.0)L=L+1 919 IF(X3(I,J).EQ.0.0)GO TO 121 920 IF(X3(I,J).GE.7.0)X3(I,J)=0.0 921 S=DSQRT((XO(5*I-4)-XO(5*J-4))**2+(XO(5*I-3)-XO(5*J-3))**2) 922 S=DSQRT(S**2+(XO(5*J-2)-XO(5*I-2))**2) 923 VS=(XO(5*J-3)-XO(5*I-3))*DCOS(XW(I))-(XO(5*J-4)-XO(5*I-4))* 924 *DSIN(XW(I)) 925 VC=(XO(5*J-2)-XO(5*I-2))*DCOS(XN(I))-(XO(5*J-4)-XO(5*I-4))*DSIN(XN 926 *(I))*DCOS(XW(I))-(XO(5*J-3)-XO(5*I-3))*DSIN(XN(I))*DSIN(XW(I)) 927 AZ=DATAN(VS/VC) 928 IF(VS.LT.0.0.AND.VC.LT.0.0)AZ=AZ+RAD 929 IF(VS.LT.0.0.AND.VC.GT.0.0)AZ=2.0*RAD+AZ 930 IF(VS.GT.0.0.AND.VC.LT.0.0)AZ=RAD+AZ 931 BC=DARSIN(((XO(5*J-4)-XO(5*I-4))*DCOS(XN(I))*DCOS(XW(I)) 932 *+(XO(5*J-3)-XO(5*I-3))*DCOS(XN(I))*DSIN(XW(I))+(XO(5*J-2) 933 *-XO(5*I-2))*DSIN(XN(I)))/S) 934 BB=XN(I) 935 BO=XW(I) 936 XN(I)=X5(I) 937 XW(I)=X6(I) 938 IF(X6(I).EQ.0.0)XW(I)=BO 939 A1=(DSIN(AZ)*DSIN(XN(I))*DCOS(XW(I))-DCOS(AZ)*DSIN(XW(I)))/ 940 *(S*DCOS(BC)) 941 A2=(DSIN(AZ)*DSIN(XN(I))*DSIN(XW(I))+DCOS(AZ)*DCOS(XW(I)))/ 942 *(S*DCOS(BC)) 943 A3=-DSIN(AZ)*DCOS(XN(I))/(S*DCOS(BC)) 944 A4=DSIN(AZ)*DTAN(BC) 945 A5=DSIN(XN(I))-DCOS(AZ)*DCOS(XN(I))*DTAN(BC) 946 AG(1,5*I-4)=-A1*P(K) 947 AG(1,5*I-3)=-A2*P(K) 948 AG(1,5*I-2)=-A3*P(K) 949 AG(1,5*J-4)=A1*P(K) 950 AG(1,5*J-3)=A2*P(K) 951 AG(1,5*J-2)=A3*P(K) 952 AG(1,5*I-1)=A4*P(K) 953 AG(1,5*I)=A5*P(K) 954 IINGIL=5*NG+L 955 AG(1,IINGIL)=-1.0*P(K) 956 SH=XO(IINGIL)+X3(I,J)-RAD/180.0 957 WG(K)=(AZ-SH)*P(K) 958 XN(I)=BB 959 XW(I)=BO 960 II=I 961 JJ=J 962 GO TO 721 963 121 CONTINUE 964 721 DO 605 IH=1,NM 965 DO 606 KU=1,NM 966 606 Q(IH,KU)=Q(IH,KU)+AG(JH,IH)*AG(JH,KU)/P(K) 967 605 CONTINUE 968 DO 640 IH=1,NM 969 KU=K 970 640 UG(IH)=UG(IH)+AG(1,IH)*WG(KU)/P(K) 971 120 CONTINUE 972 N=M+1 973 M=M+NO4 974 II=0 975 JJ=0 976 DO 130 K=N,M 977 DO 633 I=1,NM 978 633 AG(1,I)=0.0 979 DO 131 I=1,NG 980 DO 131 J=1,NG 981 IF(I.LT.II)GO TO 131 982 IF(I.LE.II.AND.J.LE.JJ)GO TO 131 983 IF(I.EQ.J)GO TO 131 984 IF(X4(I,J).EQ.0.0)GO TO 131 985 S=DSQRT((XO(5*I-4)-XO(5*J-4))**2+(XO(5*I-3)-XO(5*J-3))**2) 986 S=DSQRT(S**2+(XO(5*J-2)-XO(5*I-2))**2) 987 VS=(XO(5*J-3)-XO(5*I-3))*DCOS(XW(I))-(XO(5*J-4)-XO(5*I-4))* 988 *DSIN(XW(I)) 989 VC=(XO(5*J-2)-XO(5*I-2))*DCOS(XN(I))-(XO(5*J-4)-XO(5*I-4))*DSIN(XN 990 *(I))*DCOS(XW(I))-(XO(5*J-3)-XO(5*I-3))*DSIN(XN(I))*DSIN(XW(I)) 991 ZZ=DATAN(VS/VC) 992 IF(VS.LT.0.0.AND.VC.LT.0.0)ZZ=ZZ+RAD 993 IF(VS.LT.0.0.AND.VC.GT.0.0)ZZ=2.0*RAD+ZZ 994 IF(VS.GT.0.0.AND.VC.LT.0.0)ZZ=RAD+ZZ 995 AZ=X4(I,J) 996 BC=DARSIN(((XO(5*J-4)-XO(5*I-4))*DCOS(XN(I))*DCOS(XW(I)) 997 *+(XO(5*J-3)-XO(5*I-3))*DCOS(XN(I))*DSIN(XW(I))+(XO(5*J-2) 998 *-XO(5*I-2))*DSIN(XN(I)))/S) 999 BB=XN(I) 1000 BO=XW(I) 1001 XN(I)=X5(I) 1002 XW(I)=X6(I) 1003 IF(X6(I).EQ.0.0)XW(I)=BO 1004 A1=(DSIN(AZ)*DSIN(XN(I))*DCOS(XW(I))-DCOS(AZ)*DSIN(XW(I)))/ 1005 *(S*DCOS(BC)) 1006 A2=(DSIN(AZ)*DSIN(XN(I))*DSIN(XW(I))+DCOS(AZ)*DCOS(XW(I)))/ 1007 *(S*DCOS(BC)) 1008 A3=-DSIN(AZ)*DCOS(XN(I))/(S*DCOS(BC)) 1009 A4=DSIN(AZ)*DTAN(BC) 1010 A5=DSIN(XN(I))-DCOS(AZ)*DCOS(XN(I))*DTAN(BC) 1011 AG(1,5*I-4)=-A1*P(K) 1012 AG(1,5*I-3)=-A2*P(K) 1013 AG(1,5*I-2)=-A3*P(K) 1014 AG(1,5*J-4)=A1*P(K) 1015 AG(1,5*J-3)=A2*P(K) 1016 AG(1,5*J-2)=A3*P(K) 1017 AG(1,5*I-1)=A4*P(K) 1018 AG(1,5*I)=A5*P(K) 1019 WG(K)=(ZZ-X4(I,J))*P(K) 1020 XN(I)=BB 1021 XW(I)=BO 1022 II=I 1023 JJ=J 1024 GO TO 731 1025 131 CONTINUE 1026 731 DO 607 IH=1,NM 1027 DO 608 KU=1,NM 1028 608 Q(IH,KU)=Q(IH,KU)+AG(JH,IH)*AG(JH,KU)/P(K) 1029 607 CONTINUE 1030 DO 650 IH=1,NM 1031 KU=K 1032 650 UG(IH)=UG(IH)+AG(1,IH)*WG(KU)/P(K) 1033 130 CONTINUE 1034 N=M+1 1035 M=M+NO5 1036 II=0 1037 DO 140 K=N,M 1038 DO 634 I=1,NM 1039 634 AG(1,I)=0.0 1040 DO 141 I=1,NG 1041 IF(I.LE.II)GO TO 141 1042 IF(X5(I).EQ.0.0)GO TO 141 1043 AG(1,5*I-1)=1.0*P(K) 1044 WG(K)=(XO(5*I-1)-X5(I))*P(K) 1045 II=I 1046 GO TO 741 1047 141 CONTINUE 1048 741 DO 609 IH=1,NM 1049 DO 610 KU=1,NM 1050 610 Q(IH,KU)=Q(IH,KU)+AG(JH,IH)*AG(JH,KU)/P(K) 1051 609 CONTINUE 1052 DO 660 IH=1,NM 1053 KU=K 1054 660 UG(IH)=UG(IH)+AG(1,IH)*WG(KU)/P(K) 1055 140 CONTINUE 1056 N=M+1 1057 M=M+NO6 1058 II=0 1059 DO 150 K=N,M 1060 DO 635 I=1,NM 1061 635 AG(1,I)=0.0 1062 DO 151 I=1,NG 1063 IF(I.LE.II)GO TO 151 1064 IF(NAD.EQ.1)GO TO 151 1065 IF(X6(I).EQ.0.0)X6(I)=XW(I) 1066 AG(1,5*I)=1.0*P(K) 1067 WG(K)=(XO(5*I)-X6(I))*P(K) 1068 II=I 1069 GO TO 751 1070 151 CONTINUE 1071 751 DO 611 IH=1,NM 1072 DO 612 KU=1,NM 1073 612 Q(IH,KU)=Q(IH,KU)+AG(JH,IH)*AG(JH,KU)/P(K) 1074 611 CONTINUE 1075 DO 670 IH=1,NM 1076 KU=K 1077 670 UG(IH)=UG(IH)+AG(1,IH)*WG(KU)/P(K) 1078 150 CONTINUE 1079 30 FORMAT(2F20.9) 1080 9 FORMAT(2I9,F10.3,4F12.6,I3) 1081 11 FORMAT(3I10,2F20.9) 1082 RETURN 1083 END 1084 SUBROUTINE PNORM1(L,XP,YP,XG,YG,ZG,X,NU,NQ,NFI,AP1,AP1U,MP,W,NP11U 1085 *,NP11,INB,MPA,N1,NUM,NAM,UP1,LG,M,BP,PX,PY,PXY,AP2,NP21,JJS, 1086 *NP22,UP2,MBW,B,VEC,ITR,XF,YF,ICHT,IS2,INC,INU,PX1,NP2U,AL,ALG) 1087 C 1088 C ************************************************************************ 1089 C *SUBROUTINE PNORM1 * 1090 C * * 1091 C * THIS SUBROUTINE FORMS THE PHOTOGRAMMETRIC NORMAL * 1092 C * EQUATIONS ((NP11))-((NP21))-((NP22)) * 1093 C * || || || * 1094 C * $$$$$$$$$$$$$$$$$$$$$$$$$ * 1095 C * AND THE VECTORS UP1 AND UP2 * 1096 C * $$$$$$$$$$$ $ 1097 C * WHERE THE INPUT JCS=3*NO. OF OBJECT POINTS * 1098 C * **CPU = 24 SEC. FOR 30 PHOTOS & 200 POINTS * 1099 C * * 1100 C * WRITTEN BY S. EL-HAKIM JULY16,77 * 1101 C ************************************************************************ 1102 C 1103 IMPLICIT REAL*8(A-H,O-Z) 1104 DIMENSION L(NQ),XP(NQ),YP(NQ),XG(NQ),YG(NQ),ZG(NQ),LG(NQ) 1105 REAL*8 AP1(2,NUM),AP1U(2,N1),MP(2,2),MPA(2,NUM),W(2),UP1(NAM), 1106 *NP11(06,NUM),NP11U(NAM,N1),X(NAM),M(3,3),BP(2,2),NP21(JJS,006), 1107 *AP2(2,JJS),NP22(3,JJS),UP2(JJS),B(NAM),VEC(JJS),XF(NQ),YF(NQ) 1108 *,PX1(6),NP2U(JJS,11) 1109 *,AL(NQ),ALG(NQ) 1110 IWP=6 IL2=22 1112 IS2=1 1113 IZO=0 1114 WRITE(IWP,50) 50 FORMAT(// 5X,'RESIDUALS AT IMAGE COORDS. :'/5X,28('-')/17X,'PT.', *10X,'VX',13X,'VY'/) JCS=INC 1115 NUN=INU 1116 NAN=NUN+11 1117 NP=INB 1118 NBW=NUN 1119 REWINDIL2 1120 INDX=1 1121 SEGX=0.0 1122 SEGY=0.0 1123 DO 364 I=1,JCS 1124 DO 364 J=1,11 1125 364 NP2U(I,J)=0.0 1126 DO 365 I=1,JCS 1127 DO 365 J=1,6 1128 365 NP21(I,J)=0.0 1129 IF(ITR.GT.1)GO TO 3073 1130 DO 3074 I=1,NP 1131 XF(I)=XP(I) 1132 3074 YF(I)=YP(I) 1133 3073 CONTINUE 1134 DO1262 J=1,JCS 1135 1262 UP2(J)=0.0 1136 IF(ICHT.EQ.1)GO TO 200 1137 DO 262 J=1,JCS 1138 DO 262 I=1,3 1139 262 NP22(I,J)=0.0 1140 IDO=NBW+N1 1141 200 J=1 1142 100 IX=(AL(J)+0.01)/1000000000 1143 NUONFI=NUN 1144 NUOIX=NU*IX 1145 IIO=0 1146 IF(IX.GT.INDX)IIO=1 1147 IF(IX.GT.INDX)INDX=IX 1148 XR=XP(J)-X(NUONFI+1) 1149 YR=YP(J)-X(NUONFI+2) 1150 R=DSQRT(XR**2+YR**2) 1151 KO=IX*NU 1152 SO=DSIN(X(NUOIX-2)) 1153 CO=DCOS(X(NUOIX-2)) 1154 CF=DCOS(X(NUOIX-1)) 1155 SF=DSIN(X(NUOIX-1)) 1156 CK=DCOS(X(NUOIX)) 1157 SK=DSIN(X(NUOIX)) 1158 M(1,1)=CF*CK 1159 M(1,2)=CO*SK+SO*SF*CK 1160 M(1,3)=SO*SK-SF*CK*CO 1161 M(2,1)=-CF*SK 1162 M(2,2)=CO*CK-SO*SF*SK 1163 M(2,3)=SO*CK+SF*SK*CO 1164 M(3,1)=SF 1165 M(3,2)=-SO*CF 1166 M(3,3)=CO*CF 1167 T1=M(1,1)*(XG(J)-X(NUOIX-5))+M(1,2)*(YG(J)-X(NUOIX-4))+M(1,3)* 1168 *(ZG(J)-X(NUOIX-3)) 1169 T2=M(2,1)*(XG(J)-X(NUOIX-5))+M(2,2)*(YG(J)-X(NUOIX-4))+M(2,3)* 1170 *(ZG(J)-X(NUOIX-3)) 1171 T3=M(3,1)*(XG(J)-X(NUOIX-5))+M(3,2)*(YG(J)-X(NUOIX-4))+M(3,3)* 1172 *(ZG(J)-X(NUOIX-3)) 1173 SL=YR/R 1174 CL=XR/R 1175 S2L=2.*SL*CL 1176 C2L=2.*CL**2-1.0 1177 NU=NUN 1178 T=X(NU+4)+X(NU+5)*CL+X(NU+6)*SL+X(NU+7)*R+X(NU+8)*R*C2L+X(NU+9)*R* 1179 *S2L+X(NU+10)*R**2*CL+X(NU+11)*R**2*SL 1180 NU=6 1181 DVX=T*XR 1182 DVY=T*YR 1183 IF(ICHT.EQ.1)GO TO 300 1184 AP1(1,KO-5)=-(XR+DVX)*M(3,1)-X(NUONFI+3)*M(1,1) 1185 AP1(2,KO-5)=-(YR+DVY)*M(3,1)-X(NUONFI+3)*M(2,1) 1186 AP1(1,KO-4)=-(XR+DVX)*M(3,2)-X(NUONFI+3)*M(1,2) 1187 AP1(2,KO-4)=-(YR+DVY)*M(3,2)-X(NUONFI+3)*M(2,2) 1188 AP1(1,KO-3)=-(XR+DVX)*M(3,3)-X(NUONFI+3)*M(1,3) 1189 AP1(2,KO-3)=-(YR+DVY)*M(3,3)-X(NUONFI+3)*M(2,3) 1190 AP1(1,KO-2)=(XR+DVX)*(-M(3,3)*(YG(J)-X(NUOIX-4))+M(3,2)*(ZG(J)- 1191 *X(NUOIX-3)))+X(NUONFI+3)*(-M(1,3)*(YG(J)-X(NUOIX-4))+M(1,2)*(ZG(J) 1192 *-X(NUOIX-3))) 1193 AP1(2,KO-2)=(YR+DVY)*(-M(3,3)*(YG(J)-X(NUOIX-4))+M(3,2)*(ZG(J)- 1194 *X(NUOIX-3)))+X(NUONFI+3)*(-M(2,3)*(YG(J)-X(NUOIX-4))+M(2,2)*(ZG(J) 1195 *-X(NUOIX-3))) 1196 AP1(1,KO-1)=(XR+DVX)*(CF*(XG(J)-X(NUOIX-5))+SO*SF*(YG(J)-X(NUOIX-4 1197 *))-CO*SF*(ZG(J)-X(NUOIX-3)))+X(NUONFI+3)*(-SF*CK*(XG(J)-X(NUOIX-5) 1198 *)+SO*CF*(YG(J)-X(NUOIX-4))*CK-CF*CK*CO*(ZG(J)-X(NUOIX-3))) 1199 AP1(2,KO-1)=(YR+DVY)*(CF*(XG(J)-X(NUOIX-5))+SO*SF*(YG(J)-X(NUOIX-4 1200 *))-CO*SF*(ZG(J)-X(NUOIX-3)))+X(NUONFI+3)*(SF*SK*(XG(J)-X(NUOIX-5)) 1201 *-SO*CF*SK*(YG(J)-X(NUOIX-4))+CF*SK*CO*(ZG(J)-X(NUOIX-3))) 1202 AP1(1,KO)=X(NUONFI+3)*(M(2,1)*(XG(J)-X(NUOIX-5))+M(2,2)*(YG(J)- 1203 *X(NUOIX-4))+M(2,3)*(ZG(J)-X(NUOIX-3))) 1204 AP1(2,KO)=X(NUONFI+3)*(-M(1,1)*(XG(J)-X(NUOIX-5))-M(1,2)*(YG(J)- 1205 *X(NUOIX-4))-M(1,3)*(ZG(J)-X(NUOIX-3))) 1206 AP1U(1,1)=-T3-T*T3 1207 AP1U(2,1)=0.0 1208 AP1U(1,2)=0.0 1209 AP1U(2,2)=AP1U(1,1) 1210 AP1U(1,3)=T1 1211 AP1U(2,3)=T2 1212 AP1U(1,4)=XR*T3 1213 AP1U(2,4)=YR*T3 1214 AP1U(1,5)=XR*T3*CL 1215 AP1U(2,5)=YR*T3*CL 1216 AP1U(1,6)=XR*T3*SL 1217 AP1U(2,6)=YR*T3*SL 1218 AP1U(1,7)=XR*T3*R 1219 AP1U(2,7)=YR*T3*R 1220 AP1U(1,8)=XR*T3*R*C2L 1221 AP1U(2,8)=YR*T3*R*C2L 1222 AP1U(1,9)=XR*T3*R*S2L 1223 AP1U(2,9)=YR*T3*R*S2L 1224 AP1U(1,10)=XR*T3*R**2*CL 1225 AP1U(2,10)=YR*T3*R**2*CL 1226 AP1U(1,11)=XR*T3*R**2*SL 1227 AP1U(2,11)=YR*T3*R**2*SL 1228 300 LL=(AL(J)+0.01)/1000000000 1229 A1L=LL 1230 AR=AL(J)-A1L*1000000000.0+0.01 1231 ISO=AR 1232 KN=ISO 1233 ISO=ISO/1000 1234 KN=KN-KN/1000*1000 1235 AP2(1,3*KN-2)=(XR+DVX)*M(3,1)+X(NUN+3)*M(1,1) 1236 AP2(2,3*KN-2)=(YR+DVY)*M(3,1)+X(NUN+3)*M(2,1) 1237 AP2(1,3*KN-1)=(XR+DVX)*M(3,2)+X(NUN+3)*M(1,2) 1238 AP2(2,3*KN-1)=(YR+DVY)*M(3,2)+X(NUN+3)*M(2,2) 1239 AP2(1,3*KN)=(XR+DVX)*M(3,3)+X(NUN+3)*M(1,3) 1240 AP2(2,3*KN)=(YR+DVY)*M(3,3)+X(NUN+3)*M(2,3) 1241 T5=X(NUN+5)*SL**2/R-X(NUN+6)*CL*SL/R+X(NUN+7)*CL+X(NUN+8)*(CL*C2L+ 1242 *2.*S2L*SL)+X(NUN+9)*(CL*S2L-2.*SL*C2L)+X(NUN+10)*(XP(J)*CL+ 1243 *R)+X(NUN+11)*XP(J)*SL 1244 T6=-X(NUN+5)*SL*CL/R+X(NUN+6)*CL**2/R+X(NUN+7)*SL+X(NUN+8)*(SL*C2L 1245 *-2.*S2L*CL)+X(NUN+9)*(SL*S2L+2.*C2L*CL)+X(NUN+10)*XP(J)*SL+ 1246 *X(NUN+11)*(R+YP(J)*SL) 1247 BP(1,1)=T3+T3*(T+XR*T5) 1248 BP(1,2)=T3*(XR*T6) 1249 BP(2,1)=T3*YR*T5 1250 BP(2,2)=T3+T3*(T+YR*T6) 1251 AMP11=BP(1,1)*(PX*BP(1,1)+PXY*BP(1,2))+BP(1,2)*(PXY*BP(1,1)+PY* 1252 *BP(1,2)) 1253 AMP12=BP(1,1)*(PX*BP(2,1)+PXY*BP(2,2))+BP(1,2)*(PXY*BP(2,1)+PY* 1254 *BP(2,2)) 1255 AMP21=BP(2,1)*(PX*BP(1,1)+PXY*BP(1,2))+BP(2,2)*(PXY*BP(1,1)+PY* 1256 *BP(1,2)) 1257 AMP22=BP(2,1)*(PX*BP(2,1)+PXY*BP(2,2))+BP(2,2)*(PXY*BP(2,1)+PY*BP 1258 *(2,2)) 1259 DET=AMP11*AMP22-AMP12*AMP21 1260 MP(1,1)=AMP22/DET 1261 MP(1,2)=-AMP12/DET 1262 MP(2,1)=-AMP21/DET 1263 MP(2,2)=AMP11/DET 1264 W1=(XR+DVX)*T3+X(NUN+3)*T1 1265 W2=(YR+DVY)*T3+X(NUN+3)*T2 1266 IF(IS2.EQ.0)GO TO 600 1267 WX=W1+AP1(1,KO-5)*B(KO-5)+AP1(1,KO-4)*B(KO-4)+AP1(1,KO-3)*B(KO-3) 1268 *+AP1(1,KO-2)*B(KO-2)+AP1(1,KO-1)*B(KO-1)+AP1(1,KO)*B(KO) 1269 WY=W2+AP1(2,KO-5)*B(KO-5)+AP1(2,KO-4)*B(KO-4)+AP1(2,KO-3)*B(KO-3) 1270 *+AP1(2,KO-2)*B(KO-2)+AP1(2,KO-1)*B(KO-1)+AP1(2,KO)*B(KO) 1271 WX=WX-AP2(1,3*KN-2)*VEC(3*KN-2)-AP2(1,3*KN-1)*VEC(3*KN-1) 1272 *-AP2(1,3*KN)*VEC(3*KN) 1273 WY=WY-AP2(2,3*KN-2)*VEC(3*KN-2)-AP2(2,3*KN-1)*VEC(3*KN-1) 1274 *-AP2(2,3*KN)*VEC(3*KN) 1275 VX=MP(1,1)*WX+MP(1,2)*WY 1276 VY=MP(2,1)*WX+MP(2,2)*WY 1277 CX=-PX*(BP(1,1)*VX+BP(2,1)*VY) 1278 CY=-PY*(BP(1,2)*VX+BP(2,2)*VY) 1279 2973 FORMAT(I20,2F15.6) 1280 WRITE(IWP,2973)ISO,CX,CY 1281 SEGX=SEGX+CX**2 1282 SEGY=SEGY+CY**2 1283 600 W(1)=W1*MP(1,1)+W2*MP(1,2) 1284 W(2)=W1*MP(2,1)+W2*MP(2,2) 1285 IF(ICHT.EQ.1)GO TO 400 1286 DO 110 II=1,2 1287 DO 120 K=1,NU 1288 I=IX*NU-NU+K 1289 MPA(II,I)=0.0 1290 DO 120 JJ=1,2 1291 120 MPA(II,I)=MPA(II,I)+MP(II,JJ)*AP1(JJ,I) 1292 110 CONTINUE 1293 IF(IIO.EQ.0)GO TO 5150 1294 1005 DO 6150 KL=1,6 1295 WRITE(IL2)(NP21(JJ,KL),JJ=1,JCS) 1296 6150 CONTINUE 1297 IF(IZO.EQ.1)GO TO 1000 1298 DO 363 I=1,JCS 1299 DO 363 K=1,6 1300 363 NP21(I,K)=0.0 1301 5150 CONTINUE 1302 DO 150 II=1,3 1303 I2=KN*3-3+II 1304 DO 160 K=1,NU 1305 KK=IX*NU-NU+K 1306 KD=KK 1307 DO 161 JJ=1,2 1308 161 NP21(I2,K)=NP21(I2,K)+AP2(JJ,I2)*MPA(JJ,KK) 1309 160 CONTINUE 1310 150 CONTINUE 1311 DO 130 II=1,NU 1312 I2=IX*NU-NU+II 1313 DO 140 K=1,NU 1314 KK=IX*NU-NU+K 1315 DO 140 JJ=1,2 1316 140 NP11(II,KK)=NP11(II,KK)+AP1(JJ,I2)*MPA(JJ,KK) 1317 130 CONTINUE 1318 DO 115 II=1,2 1319 DO 125 K=1,N1 1320 MPA(II,K)=0.0 1321 DO 125 JJ=1,2 1322 125 MPA(II,K)=MPA(II,K)+MP(II,JJ)*AP1U(JJ,K) 1323 115 CONTINUE 1324 I=NBW+1 1325 N2=NBW+N1 1326 IK=0 1327 DO 170 II=1,3 1328 I2=KN*3-3+II 1329 IO=0 1330 DO 180 K=I,N2 1331 IO=IO+1 1332 DO 180 JJ=1,2 1333 180 NP2U(I2,IO)=NP2U(I2,IO)+AP2(JJ,I2)*MPA(JJ,IO) 1334 170 CONTINUE 1335 I=NUN+1 1336 N2=NUN+N1 1337 IK=0 1338 DO 135 II=I,N2 1339 IK=IK+1 1340 IO=0 1341 DO 145 K=I,N2 1342 IO=IO+1 1343 DO 145 JJ=1,2 1344 145 NP11U(II,IO)=NP11U(II,IO)+AP1U(JJ,IK)*MPA(JJ,IO) 1345 135 CONTINUE 1346 DO 136 II=1,NU 1347 IK=IX*NU-NU+II 1348 DO 137 K=1,N1 1349 DO 137 JJ=1,2 1350 137 NP11U(IK,K)=NP11U(IK,K)+AP1(JJ,IK)*MPA(JJ,K) 1351 136 CONTINUE 1352 DO 210 II=1,NU 1353 IK=IX*NU-NU+II 1354 DO 210 KK=1,2 1355 210 UP1(IK)=UP1(IK)+AP1(KK,IK)*W(KK) 1356 IFO=NU*NFI+1 1357 IFE=IFO+N1-1 1358 IK=0 1359 DO 310 II=IFO,IFE 1360 IK=IK+1 1361 DO 310 KK=1,2 1362 310 UP1(II)=UP1(II)+AP1U(KK,IK)*W(KK) 1363 400 CONTINUE 1364 DO 111 II=1,2 1365 DO 121 K=1,3 1366 MPA(II,K)=0.0 1367 KK=KN*3-3+K 1368 DO 121 JJ=1,2 1369 121 MPA(II,K)=MPA(II,K)+MP(II,JJ)*AP2(JJ,KK) 1370 111 CONTINUE 1371 IF(ICHT.EQ.1)GO TO 500 1372 DO 131 II=1,3 1373 I2=KN*3-3+II 1374 DO 141 K=1,3 1375 KK=KN*3-3+K 1376 DO 141 JJ=1,2 1377 141 NP22(II,KK)=NP22(II,KK)+AP2(JJ,I2)*MPA(JJ,K) 1378 131 CONTINUE 1379 500 DO 211 II=1,3 1380 IK=KN*3-3+II 1381 DO 211 KK=1,2 1382 211 UP2(IK)=UP2(IK)+AP2(KK,IK)*W(KK) 1383 IF(J.EQ.NP)IZO=1 1384 IF(J.EQ.NP)GO TO 1005 1385 J=J+1 1386 GO TO 100 1387 1000 DF=NP-NUN 1388 SEGX=DSQRT(SEGX/DF) 1389 SEGY=DSQRT(SEGY/DF) 1390 WRITE(IWP,1001)SEGX,SEGY 1391 1001 FORMAT(//5X,'ST. ERRORS'/2F20.6//) 1392 SEG=SEGX+SEGY 1393 DO 1002 I=1,NUN 1394 K=I-(I-1)/6*6 1395 1002 SEG=SEG+B(I)**2*PX1(K) 1396 DF=NP 1397 ENDFILEIL2 1398 RETURN 1399 END 1400 SUBROUTINE NORM2(NP22,NG,PX2,NM,JJS,JC,IG,SE,M,MGW,NAM,IRM,INC 1401 *,INU) 1402 IMPLICIT REAL*8(A-H,O-Z) 1403 REAL*8 NP22(3,JJS),NG(NAM,NAM),PX2(3,JJS),SE(NAM),M(3,3) 1404 IF(IRM.GE.3)NM=3*IG 1405 JCS=INC 1406 NUN=INU 1407 NAN=NUN+11 1408 DO 100 I=1,JC 1409 III=3*I 1410 NP22(1,III-2)=NP22(1,III-2)+PX2(1,III-2) 1411 NP22(2,III-1)=NP22(2,III-1)+PX2(2,III-1) 1412 100 NP22(3,III-0)=NP22(3,III-0)+PX2(3,III) 1413 IF(IG.EQ.0)GO TO 1000 1414 DO 200 I=1,IG 1415 KII=5*I 1416 III=3*I 1417 IF(IRM.GE.3)KII=III+2 1418 NG(KII-4,KII-4)=NG(KII-4,KII-4)+NP22(1,III-2) 1419 NG(KII-4,KII-3)=NG(KII-4,KII-3)+NP22(1,III-1) 1420 NG(KII-4,KII-2)=NG(KII-4,KII-2)+NP22(1,III) 1421 NG(KII-3,KII-4)=NG(KII-3,KII-4)+NP22(2,III-2) 1422 NG(KII-3,KII-3)=NG(KII-3,KII-3)+NP22(2,III-1) 1423 NG(KII-3,KII-2)=NG(KII-3,KII-2)+NP22(2,III) 1424 NG(KII-2,KII-4)=NG(KII-2,KII-4)+NP22(3,III-2) 1425 NG(KII-2,KII-3)=NG(KII-2,KII-3)+NP22(3,III-1) 1426 NG(KII-2,KII-2)=NG(KII-2,KII-2)+NP22(3,III) 1427 200 CONTINUE 1428 CALL SYINV(NG,NAM,NM,SE) 1429 1000 CONTINUE 1430 IF(IRM.EQ.1.OR.IRM.EQ.4)RETURN 1431 K=IG+1 1432 DO 300 I=K,JC 1433 III=3*I 1434 DO 310 J=1,3 1435 M(J,1)=NP22(J,III-2) 1436 M(J,2)=NP22(J,III-1) 1437 310 M(J,3)=NP22(J,III) 1438 CALL SYINV3(M) 1439 DO 320 L=1,3 1440 NP22(L,III-2)=M(L,1) 1441 NP22(L,III-1)=M(L,2) 1442 320 NP22(L,III)=M(L,3) 1443 300 CONTINUE 1444 RETURN 1445 END 1446 SUBROUTINE NORMAL(JC,JJS,NUM,NAM,NG,NP22,NP11,NP11U,PX1,PX1U,NE 1447 *,IG,VEC,NP21,INC,INU,NP2U,XP21,NO,IG1,IGL) 1448 IMPLICIT REAL*8(A-H,O-Z) 1449 REAL*8 NG(NAM,NAM),NP22(3,JJS),NP11(6,NUM),NP11U(NAM,11),VEC(JJS) 1450 *,PX1(6),PX1U(11),NE(NAM,1),NP21(JJS,006),NP2U(JJS,11),XP21(JJS,6) 1451 *,NO(IG1,IG1) 1452 IL2=22 1453 REWINDIL2 1454 JCS=INC 1455 NUN=INU 1456 NAN=NUN+11 1457 IF(IG.EQ.0.OR.IGL.EQ.1)GO TO 101 1458 DO 100 I=1,IG 1459 I3=3*I 1460 I5=5*I 1461 DO 100 J=1,IG 1462 J3=3*J 1463 J5=5*J 1464 NO(I3-2,J3-2)=NG(I5-4,J5-4) 1465 NO(I3-2,J3-1)=NG(I5-4,J5-3) 1466 NO(I3-2,J3-0)=NG(I5-4,J5-2) 1467 NO(I3-1,J3-2)=NG(I5-3,J5-4) 1468 NO(I3-1,J3-1)=NG(I5-3,J5-3) 1469 NO(I3-1,J3-0)=NG(I5-3,J5-2) 1470 NO(I3-0,J3-2)=NG(I5-2,J5-4) 1471 NO(I3-0,J3-1)=NG(I5-2,J5-3) 1472 NO(I3-0,J3-0)=NG(I5-2,J5-2) 1473 100 CONTINUE 1474 101 CONTINUE 1475 IT=IG+1 1476 NFI=NUN/6 1477 IXI=NUN+1 1478 L=1 1479 READ(IL2)(XP21(KK,L),KK=1,JCS) 1480 DO 10 K=1,NAN 1481 IF(K.GT.NUN)GO TO 200 1482 L=K-(K-1)/6*6 1483 DO 83 IQ=1,JCS 1484 83 NP21(IQ,L)=XP21(IQ,L) 1485 200 CONTINUE 1486 IF(IG.EQ.0)GO TO 222 1487 DO 20 J=1,IG 1488 DO 20 II=1,3 1489 JJ=3*J-3+II 1490 VEC(JJ)=0.0 1491 IF(K.GT.NUN)GO TO 250 1492 DO 15 I=1,IG 1493 15 VEC(JJ)=VEC(JJ)+NP21(3*I-2,L)*NO(JJ,3*I-2)+NP21(3*I-1,L)* 1494 *NO(JJ,3*I-1)+NP21(3*I,L)*NO(JJ,3*I-0) 1495 GO TO 20 1496 250 L=K-NUN 1497 DO 16 I=1,IG 1498 16 VEC(JJ)=VEC(JJ)+NP2U(3*I-2,L)*NO(JJ,3*I-2)+NP2U(3*I-1,L)* 1499 *NO(JJ,3*I-1)+NP2U(3*I,L)*NO(JJ,3*I-0) 1500 20 CONTINUE 1501 222 CONTINUE 1502 DO 25 J=IT,JC 1503 DO 25 II=1,3 1504 JJ=3*J-3+II 1505 IF(K.GT.NUN)GO TO 350 1506 VEC(JJ)=NP21(3*J-2,L)*NP22(II,3*J-2)+NP21(3*J-1,L)*NP22(II,3*J-1) 1507 *+NP21(3*J,L)*NP22(II,3*J) 1508 GO TO 25 1509 350 L=K-NUN 1510 VEC(JJ)=NP2U(3*J-2,L)*NP22(II,3*J-2)+NP2U(3*J-1,L)*NP22(II,3*J-1) 1511 *+NP2U(3*J,L)*NP22(II,3*J) 1512 25 CONTINUE 1513 REWINDIL2 1514 II=0 1515 IL=K+1 1516 DO 40 I=1,NFI 1517 DO 40 J=1,6 1518 II=II+1 1519 NE(II,1)=0.0 1520 IF(II.GE.K)GO TO 401 1521 READ(IL2) 1522 GO TO 40 1523 401 CONTINUE 1524 READ(IL2)(NP21(M,J),M=1,JCS) 1525 IF(K.GT.NUN)GO TO 81 1526 IF(II.NE.IL)GO TO 81 1527 DO 82 IQ=1,JCS 1528 82 XP21(IQ,J)=NP21(IQ,J) 1529 81 CONTINUE 1530 DO 41 KK=1,JCS 1531 41 NE(II,1)=NE(II,1)-NP21(KK,J)*VEC(KK) 1532 40 CONTINUE 1533 II=II+1 1534 IO=0 1535 DO 42 I=II,NAN 1536 IO=IO+1 1537 NE(I,1)=0.0 1538 DO 42 J=1,JCS 1539 42 NE(I,1)=NE(I,1)-NP2U(J,IO)*VEC(J) 1540 IK=K-(K-1)/6*6 1541 KF=(K-1)/6+1 1542 IF(K.GT.NUN)GO TO 50 1543 DO 43 III=1,6 1544 JJ=KF*6-6+III 1545 43 NE(JJ,1)=NP11(III,K)+NE(JJ,1) 1546 NE(K,1)=NE(K,1)+PX1(IK) 1547 IO=0 1548 DO 44 I=IXI,NAN 1549 IO=IO+1 1550 44 NE(I,1)=NP11U(K,IO)+NE(I,1) 1551 GO TO 60 1552 50 IO=0 1553 IK=K-NUN 1554 DO 70 I=IXI,NAN 1555 IO=IO+1 1556 70 NE(I,1)=NP11U(K,IO)+NE(I,1) 1557 NE(K,1)=NE(K,1)+PX1U(IK) 1558 60 DO 120 I=1,NAN 1559 120 NG(I,K)=NE(I,1) 1560 10 CONTINUE 1561 RETURN 1562 END 1563 SUBROUTINE SYINV3(A) 1564 IMPLICIT REAL*8(A-H,O-Z) 1565 REAL*8 A(3,3) 1566 DET=A(1,1)*A(2,2)-A(1,2)**2 1567 PX=A(2,2)/DET 1568 PY=A(1,1)/DET 1569 PXY=-A(1,2)/DET 1570 A1=A(1,3)*PX+A(2,3)*PXY 1571 A2=A(1,3)*PXY+A(2,3)*PY 1572 EP=A(3,3)-A1*A(3,1)-A2*A(3,2) 1573 A(1,1)=PX+A1/EP*A1 1574 A(1,2)=PXY+A1/EP*A2 1575 A(2,1)=PXY+A2/EP*A1 1576 A(2,2)=PY+A2/EP*A2 1577 A(1,3)=-A1/EP 1578 A(2,3)=-A2/EP 1579 A(3,1)=A(1,3) 1580 A(3,2)=A(2,3) 1581 A(3,3)=1.0/EP 1582 RETURN 1583 END 1584 SUBROUTINE SYINV(A,N1,NN,SE) 1585 IMPLICIT REAL*8(A-H,O-Z) 1586 REAL*8 A(N1,N1),SE(N1) 1587 DET=A(1,1)*A(2,2)-A(1,2)**2 1588 PX=A(2,2)/DET 1589 PY=A(1,1)/DET 1590 PXY=-A(1,2)/DET 1591 A1=A(1,3)*PX+A(2,3)*PXY 1592 A2=A(1,3)*PXY+A(2,3)*PY 1593 EP=A(3,3)-A1*A(3,1)-A2*A(3,2) 1594 A(1,1)=PX+A1/EP*A1 1595 A(1,2)=PXY+A1/EP*A2 1596 A(2,1)=PXY+A2/EP*A1 1597 A(2,2)=PY+A2/EP*A2 1598 A(1,3)=-A1/EP 1599 A(2,3)=-A2/EP 1600 A(3,1)=A(1,3) 1601 A(3,2)=A(2,3) 1602 A(3,3)=1.0D00/EP 1603 N=3 1604 200 N=N+1 1605 M=N-1 1606 DO 11 I=1,M 1607 SE(I)=0.0D00 1608 DO 10 J=1,M 1609 10 SE(I)=SE(I)+A(I,J)*A(J,N) 1610 11 CONTINUE 1611 EP=0.0D00 1612 DO 12 J=1,M 1613 12 EP=EP+A(N,J)*SE(J) 1614 EB=A(N,N)-EP 1615 EI=1.0D00/EB 1616 DO 100 I=1,M 1617 DO 100 J=1,M 1618 100 A(I,J)=A(I,J)+SE(I)*SE(J)*EI 1619 DO 150 I=1,M 1620 A(I,N)=-SE(I)*EI 1621 150 A(N,I)=A(I,N) 1622 A(N,N)=EI 1623 IF(NN.GT.N)GO TO 200 1624 RETURN 1625 END 1626 SUBROUTINE SNORM(NG,UG,NAN,NP,XG,YG,ZG,IG,ITR,S,WT,INB,IGM,L,XO, *IGN,AG,WG,AL) IMPLICIT REAL*8(A-H,O-Z) REAL*8 NG(IGN,IGN),UG(NAN),XG(NP),YG(NP),ZG(NP),S(IGM,IGM), *WT(IGM,IGM),XO(NAN),AG(1,NAN),AL(NP) DIMENSION L(NP) 10 FORMAT(6X,2I5,4X,2F10.3) IRC=5 IWP=6 IF(ITR.EQ.1)WRITE(IWP,11) 11 FORMAT(//5X,'READ IN DIST. OBSERVATIONS'/11X,'FROM',8X,'TO',10X, *'ST.E.',10X,'DIST.'/) 111 FORMAT(//5X,'DISTANCES FROM ADJ. COORDS'/11X,'FROM',8X,'TO',10X, *'DIFF.',10X,'DIST.'/) 12 FORMAT(5X,2I10,5X,F10.3,5X,F10.3) IGS=IG*3 DO 20 II=1,IGS UG(II)=0.0 DO 20 JJ=1,IGS 20 NG(II,JJ)=0.0 IF(ITR.NE.1)GO TO 40 DO 33 K=1,IG DO 33 M=1,IG 33 S(K,M)=0.0 READ(IRC,10)NO DO 77 KM=1,NO READ(IRC,10)I,J,P,SS WRITE(IWP,12)I,J,P,SS IF(I.EQ.-999)GO TO 40 I1=0 I2=0 DO 105 II=1,INB IR=L(II)-L(II)/1000000*1000000 IF(IR.EQ.I)GO TO 104 IF(IR.EQ.J)GO TO 106 GO TO 105 104 LL=(AL(II)+0.01)/1000000000 A1L=LL AR=AL(II)-A1L*1000000000.+.01 IR=AR K=IR-IR/1000*1000 I1=1 KK=3*K XO(KK-2)=XG(II) XO(KK-1)=YG(II) XO(KK)=ZG(II) GO TO 107 106 LL=(AL(II)+0.01)/1000000000 A1L=LL AR=AL(II)-A1L*1000000000.+.01 IR=AR M=IR-IR/1000*1000 I2=1 MM=3*M XO(MM-2)=XG(II) XO(MM-1)=YG(II) XO(MM)=ZG(II) 107 IF(I1.EQ.1.AND.I2.EQ.1)GO TO 108 105 CONTINUE 108 CONTINUE S(K,M)=SS WT(K,M)=1/P*2.0 77 CONTINUE WRITE(IWP,13)IG,NO 13 FORMAT(/5X,'NO. OF GEODETIC PTS. =',I4/5X,'NO. OF GEODETIC OBS. =' *,I4) 40 CONTINUE IF(ITR.EQ.1)GO TO 45 DO 48 II=1,IG III=3*II DO 505 I=1,INB LM=AL(I)/1000. AJ=LM LL=AL(I)-AJ*1000+0.01 IF(LL.EQ.II)GO TO 506 505 CONTINUE 506 CONTINUE XO(III-2)=XG(I) XO(III-1)=YG(I) XO(III-0)=ZG(I) 48 CONTINUE 45 CONTINUE II=0 JJ=0 WRITE(IWP,111) DO 100 KK=1,NO DO 94 I=1,IGS 94 AG(1,I)=0.0 DO 101 I=1,IG DO 101 J=1,IG IF(I.LT.II)GO TO 101 IF(I.LE.II.AND.J.LE.JJ)GO TO 101 IF(I.EQ.J)GO TO 101 IF(S(I,J).EQ.0.0)GO TO 101 I3=3*I J3=3*J SC=DSQRT((XO(I3-2)-XO(J3-2))**2+(XO(I3-1)-XO(J3-1))**2+(XO(I3)- *XO(J3))**2) D=SC-S(I,J) WGOKKO=D*WT(I,J) DO 97 IN=1,INB LM=AL(IN)/1000 LL=AL(IN)-LM*1000+0.01 IF(LL.EQ.I)LOIO=L(IN)-L(IN)/1000000*1000000 IF(LL.EQ.J)LOJO=L(IN)-L(IN)/1000000*1000000 97 CONTINUE WRITE(IWP,12)LOIO,LOJO,D,SC II=I JJ=J IF(ITR.EQ.0)GO TO 100 C1=(XO(J3-2)-XO(I3-2))/SC C2=(XO(J3-1)-XO(I3-1))/SC C3=(XO(J3)-XO(I3))/SC AG(1,I3-2)=-C1*WT(I,J) AG(1,I3-1)=-C2*WT(I,J) AG(1,I3)=-C3*WT(I,J) AG(1,J3-2)=-AG(1,I3-2) AG(1,J3-1)=-AG(1,I3-1) AG(1,J3)=-AG(1,I3) GO TO 701 101 CONTINUE IF(I.EQ.IG.AND.J.EQ.IG)RETURN 701 CONTINUE DO 301 IB=1,3 IA=IB-1 DO 301 JB=1,3 JA=JB-1 I3OJA=I3-JA I3OIA=I3-IA 301 NG(I3OIA,I3OJA)=NG(I3OIA,I3OJA)+AG(1,I3OIA)*AG(1,I3OJA)/WT(I,J) DO 302 IB=1,3 IA=IB-1 DO 302 JB=1,3 JA=JB-1 J3OJA=J3-JA I3OIA=I3-IA 302 NG(I3OIA,J3OJA)=NG(I3OIA,J3OJA)+AG(1,I3OIA)*AG(1,J3OJA)/WT(I,J) DO 303 IB=1,3 IA=IB-1 DO 303 JB=1,3 JA=JB-1 J3OIA=J3-IA I3OJA=I3-JA 303 NG(J3OIA,I3OJA)=NG(I3OJA,J3OIA) DO 304 IB=1,3 DO 304 JB=1,3 IA=IB-1 JA=JB-1 J3OJA=J3-JA J3OIA=J3-IA 304 NG(J3OIA,J3OJA)=NG(J3OIA,J3OJA)+AG(1,J3OIA)*AG(1,J3OJA)/WT(I,J) DO 305 IB=1,3 IA=IB-1 305 UG(I3-IA)=UG(I3-IA)+AG(1,I3-IA)*WGOKKO/WT(I,J) DO 306 IB=1,3 IA=IB-1 306 UG(J3-IA)=UG(J3-IA)+AG(1,J3-IA)*WGOKKO/WT(I,J) 100 CONTINUE RETURN END